MGSpinStability3DRigidBodyKane.html  (MotionGenesis input/output).
```   (1) % MotionGenesis file:  MGSpinStability3DRigidBodyKane.txt
(3) % Purpose: Spin stability of rigid body (books, footballs, aircraft).
(4) %-------------------------------------------------------------------
(5) %       Physical objects
(6) NewtonianFrame  N
(7) RigidBody       B
(8) %-------------------------------------------------------------------
(9) Variable   wx',  wy',  wz'         % Angular velocity measures
(10) Constant   b = 0.0 N*m*s/rad       % Viscous damping constant
(11) B.SetMassInertia( m = 0.4 kg,  Ixx = 1 kg*m^2,  Iyy = 2 kg*m^2,  Izz = 3 kg*m^2 )
(12) %-------------------------------------------------------------------
(13) %       Rotational and translational kinematics.
(14) B.SetAngularVelocityAcceleration( N,  wx*Bx> + wy*By> + wz*Bz> )
-> (15) w_B_N> = wx*Bx> + wy*By> + wz*Bz>
-> (16) alf_B_N> = wx'*Bx> + wy'*By> + wz'*Bz>

(17) Bcm.SetVelocity( N, 0> )
-> (18) v_Bcm_N> = 0>

(19) %-------------------------------------------------------------------
(21) B.AddTorque( -b * B.GetAngularVelocity(N) )
-> (22) Torque_B> = -b*wx*Bx> - b*wy*By> - b*wz*Bz>

(23) %-------------------------------------------------------------------
(24) %       Form Kane's equations of motion.
(25) SetGeneralizedSpeed( wx,  wy,  wz )
(26) Dynamics = System.GetDynamicsKane()
-> (27) Dynamics[1] = b*wx + Ixx*wx' - (Iyy-Izz)*wy*wz
-> (28) Dynamics[2] = b*wy + (Ixx-Izz)*wx*wz + Iyy*wy'
-> (29) Dynamics[3] = b*wz + Izz*wz' - (Ixx-Iyy)*wx*wy

(30) Solve( Dynamics = 0,   wx',  wy',  wz'  )
-> (31) wx' = -(b*wx-(Iyy-Izz)*wy*wz)/Ixx
-> (32) wy' = -(b*wy+(Ixx-Izz)*wx*wz)/Iyy
-> (33) wz' = -(b*wz-(Ixx-Iyy)*wx*wy)/Izz

(34) %-------------------------------------------------------------------
(35) %       Numerical integration parameters and initial values.
(36) Input  tFinal = 4 sec,  tStep = 0.02 sec,  absError = 1.0E-08
(38) %-------------------------------------------------------------------
(39) %       Form output quantities and solve ODEs.
(40) H> = B.GetAngularMomentum( Bcm )
-> (41) H> = Ixx*wx*Bx> + Iyy*wy*By> + Izz*wz*Bz>

(42) Hmag = GetMagnitude( H> )
-> (43) Hmag = sqrt(Ixx^2*wx^2+Iyy^2*wy^2+Izz^2*wz^2)

(44) Wmag = B.GetAngularSpeed( N )
-> (45) Wmag = sqrt(wx^2+wy^2+wz^2)

(46) theta = acos( Dot( H>, B.GetAngularVelocity(N) ) / (Hmag * Wmag) )
-> (47) theta = acos((Ixx*wx^2+Iyy*wy^2+Izz*wz^2)/(Hmag*Wmag))