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

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

(20) %-------------------------------------------------------------------
(21) %       Form equations of motion (angular momentum principle).
(22) Dynamics[1] = Dot(  Bx>,  B.GetDynamics(Bcm)  )
-> (23) Dynamics[1] = b*wx + Ixx*wx' - (Iyy-Izz)*wy*wz

(24) Dynamics[2] = Dot(  By>,  B.GetDynamics(Bcm)  )
-> (25) Dynamics[2] = b*wy + (Ixx-Izz)*wx*wz + Iyy*wy'

(26) Dynamics[3] = Dot(  Bz>,  B.GetDynamics(Bcm)  )
-> (27) Dynamics[3] = b*wz + Izz*wz' - (Ixx-Iyy)*wx*wy

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

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

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

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

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