<s> from math import atan , cos , sin , pi <EOL> from tvtk . tools . visual import vector , MVector , Box , Helix , Frame , Cylinder , curve , color , iterate , show <EOL> top = vector ( <NUM_LIT:0> , <NUM_LIT:1.> , <NUM_LIT:0> ) <EOL> ks = <NUM_LIT> <EOL> Lspring = <NUM_LIT:1.> <EOL> Rspring = <NUM_LIT> <EOL> Dspring = <NUM_LIT> <EOL> Lshaft = <NUM_LIT:1.> <EOL> Rshaft = <NUM_LIT> <EOL> M = <NUM_LIT:1.> <EOL> Rrotor = <NUM_LIT> <EOL> Drotor = <NUM_LIT:0.1> <EOL> Dsquare = <NUM_LIT> * Drotor <EOL> I = <NUM_LIT:0.5> * M * Rrotor ** <NUM_LIT> <EOL> omega = <NUM_LIT> <EOL> g = <NUM_LIT> <EOL> Fgrav = MVector ( <NUM_LIT:0> , - M * g , <NUM_LIT:0> ) <EOL> precession = M * g * ( Lshaft / <NUM_LIT> ) / ( I * abs ( omega ) ) <EOL> phi = atan ( precession ** <NUM_LIT:2> * ( Lshaft / <NUM_LIT> ) / g ) <EOL> s = M * g / ( ks * cos ( phi ) ) <EOL> phi = <NUM_LIT:1.> / ( ( ( I * abs ( omega ) ) / ( M * Lshaft / <NUM_LIT> ) ) ** <NUM_LIT:2> / ( g * Lshaft / <NUM_LIT> ) - ( Lspring + s ) / ( Lshaft / <NUM_LIT> ) ) <EOL> s = M * g / ( ks * cos ( phi ) ) <EOL> phi = <NUM_LIT:1.> / ( ( ( I * abs ( omega ) ) / ( M * Lshaft / <NUM_LIT> ) ) ** <NUM_LIT:2> / ( g * Lshaft / <NUM_LIT> ) - ( Lspring + s ) / ( Lshaft / <NUM_LIT> ) ) <EOL> pprecess = vector ( <NUM_LIT:0> , <NUM_LIT:0> , M * precession * ( Lshaft / <NUM_LIT> + ( Lspring + s ) * sin ( phi ) ) ) <EOL> if omega < <NUM_LIT:0> : <EOL> pprecess = - pprecess <EOL> support = Box ( pos = top + MVector ( <NUM_LIT:0> , <NUM_LIT> , <NUM_LIT:0> ) , size = ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) , color = ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:0> ) ) <EOL> spring = Helix ( pos = top , axis = vector ( - <NUM_LIT> , - <NUM_LIT> , <NUM_LIT:0> ) , radius = Rspring , color = ( <NUM_LIT:1> , <NUM_LIT> , <NUM_LIT> ) ) <EOL> gyro1 = Frame ( pos = top + spring . axis ) <EOL> gyro1 . axis = MVector ( <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> shaft = Cylinder ( pos = gyro1 . pos , axis = Lshaft * gyro1 . axis , radius = <EOL> Rshaft , color = ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) , length = <NUM_LIT:1.0> ) <EOL> rotor = Cylinder ( pos = <NUM_LIT:0.5> * gyro1 . axis * ( Lshaft - Drotor ) , axis = <EOL> gyro1 . axis * Drotor , radius = Rrotor , color = ( <NUM_LIT:0.5> , <NUM_LIT:0.5> , <NUM_LIT:0.5> ) , length = <NUM_LIT:0.1> ) <EOL> stripe1 = curve ( color = color . green , <EOL> points = [ rotor . pos + <NUM_LIT> * rotor . axis + vector ( <NUM_LIT:0> , Rrotor , <NUM_LIT:0> ) , <EOL> rotor . pos + <NUM_LIT> * rotor . axis - vector ( <NUM_LIT:0> , Rrotor , <NUM_LIT:0> ) ] ) <EOL> stripe2 = curve ( color = color . green , <EOL> points = [ rotor . pos - <NUM_LIT> * rotor . axis + vector ( <NUM_LIT:0> , Rrotor , <NUM_LIT:0> ) , <EOL> rotor . pos - <NUM_LIT> * rotor . axis - vector ( <NUM_LIT:0> , Rrotor , <NUM_LIT:0> ) ] ) <EOL> gyro = Frame ( stripe1 , stripe2 ) <EOL> gyro . pos = top + spring . axis <EOL> gyro . axis = vector ( <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> gyro . rotate ( axis = ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:0> ) , angle = <NUM_LIT> , origin = gyro . pos ) <EOL> cm = gyro . pos + <NUM_LIT:0.5> * Lshaft * gyro . axis <EOL> Lrot = I * omega * gyro . axis <EOL> p = pprecess <EOL> dt = <NUM_LIT> <EOL> t = <NUM_LIT:0.> <EOL> c = <NUM_LIT:0> <EOL> def anim ( ) : <EOL> global Lrot , p , cm , t <EOL> Fspring = - ks * vector . norm ( spring . axis ) * ( vector . mag ( spring . axis ) - Lspring ) <EOL> torque = vector . cross ( - <NUM_LIT:0.5> * Lshaft * gyro . axis , Fspring ) <EOL> Lrot = Lrot + torque * dt <EOL> p = p + ( Fgrav + Fspring ) * dt <EOL> cm = cm + ( p / M ) * dt <EOL> if omega > <NUM_LIT:0> : <EOL> gyro . axis = vector . norm ( Lrot ) <EOL> else : <EOL> gyro . axis = - vector . norm ( Lrot ) <EOL> gyro . pos = cm - <NUM_LIT:0.5> * Lshaft * gyro . axis <EOL> spring . axis = gyro . pos - top <EOL> spring . length = vector . mag ( spring . axis ) <EOL> gyro . rotate ( <NUM_LIT> * omega * dt / ( <NUM_LIT> * pi ) , gyro . axis , gyro . pos ) <EOL> shaft . pos = gyro . pos <EOL> shaft . axis = Lshaft * gyro . axis <EOL> rotor . pos = gyro . pos + <NUM_LIT:0.5> * gyro . axis * ( Lshaft - Drotor ) <EOL> rotor . axis = gyro . axis * Drotor <EOL> t = t + dt <EOL> a = iterate ( <NUM_LIT> , anim ) <EOL> a . edit_traits ( ) <EOL> show ( ) </s>
