MGGyroPrecessionNutationSpinWithWxyz.html  (MotionGenesis input/output).
```   (1) % MotionGenesis file:  MGGyroPrecessionNutationSpinWithWxyz.txt
(3) %--------------------------------------------------------------
(4) NewtonianFrame  N
(5) RigidFrame      A, B              % Intermediate frames.
(6) RigidBody       C                 % Rotor.
(7) %--------------------------------------------------------------
(8) Variable   theta', phi', wC       % Angles and spin rate.
(9) Constant   g = 9.81 m/s^2         % Earth's local gravity.
(10) Constant   L = 0.2  m             % Distance between No and Ccm.
(11) Constant   r = 0.2  m             % Rotor radius.
(12) C.SetMass( m = 0.1 kg )
(13) I = m*r^2/2                       % Moment of inertia about symmetry axis.
-> (14) I = 0.5*m*r^2

(15) J = m*(r^2/4 + L^2)               % Moment of inertia about transverse axis through Co.
-> (16) J = 0.25*m*(r^2+4*L^2)

(17) C.SetInertia( Co,  J*Bx>*Bx> + J*By>*By> + I*Bz>*Bz> )
(18) %--------------------------------------------------------------------
(19) %       Rotational kinematics.
(20) A.RotateNegativeZ( N, theta )
-> (21) A_N = [cos(theta), -sin(theta), 0;  sin(theta), cos(theta), 0;  0, 0, 1]
-> (22) w_A_N> = -theta'*Az>

(23) B.RotateNegativeX( A, phi )
-> (24) B_A = [1, 0, 0;  0, cos(phi), -sin(phi);  0, sin(phi), cos(phi)]
-> (25) w_B_A> = -phi'*Bx>
-> (26) w_B_N> = -theta'*Az> - phi'*Bx>

(27) Express( w_B_N>, B )
-> (28) w_B_N> = -phi'*Bx> + sin(phi)*theta'*By> - cos(phi)*theta'*Bz>

(29) C.SetAngularVelocity( B, wC*Bz> )
-> (30) w_C_B> = wC*Bz>

(31) %--------------------------------------------------------------------
(32) %       Efficient change of variables to wx, wy, wz.
(33) Variable    wx', wy', wz'
(34) ZeroAngularVelocity> = wx*Bx> + wy*By> + wz*Bz> - C.GetAngularVelocity(N)
-> (35) ZeroAngularVelocity> = (wx+phi')*Bx> + (wy-sin(phi)*theta')*By> + (wz+
cos(phi)*theta'-wC)*Bz>

(36) ChangeAngularVelocityVariablesEqn = Matrix( B, ZeroAngularVelocity> )
-> (37) ChangeAngularVelocityVariablesEqn = [wx + phi';  wy - sin(phi)*theta';  wz + cos(phi)*theta' - wC]

(38) Solve( ChangeAngularVelocityVariablesEqn = 0,  theta', phi', wC )
-> (39) theta' = wy/sin(phi)
-> (40) phi' = -wx
-> (41) wC = wz + wy*cos(phi)/sin(phi)

(42) C.SetAngularVelocity( N, wx*Bx> + wy*By> + wz*Bz> )
-> (43) w_C_N> = wx*Bx> + wy*By> + wz*Bz>

(44) alf_C_N> = Dt( w_C_N>, C )
-> (45) alf_C_N> = (wC*wy+wx')*Bx> + (wy'-wC*wx)*By> + wz'*Bz>

(46) Explicit( w_B_N>,  wx, wy, wz )
-> (47) w_B_N> = wx*Bx> + wy*By> - wy*cos(phi)/sin(phi)*Bz>

(48) %--------------------------------------------------------------------
(49) %       Translational kinematics.
(50) Co.Translate( No, 0> )
-> (51) p_No_Co> = 0>
-> (52) v_Co_N> = 0>
-> (53) a_Co_N> = 0>

(54) CCm.Translate( Co, L*Bz> )
-> (55) p_Co_Ccm> = L*Bz>
-> (56) v_Ccm_N> = L*wy*Bx> - L*wx*By>
-> (57) a_Ccm_N> = -L*(wC*wx-wx*wz-wy')*Bx> + L*(wy*wz-wC*wy-wx')*By> - L*(wx^2
+wy^2)*Bz>

(58) %--------------------------------------------------------------------
(59) %       Add relevant forces (gravity).
-> (61) Force_Ccm> = -m*g*Nz>

(62) %--------------------------------------------------------------------
(63) %       Form equations of motion (angular momentum principle).
(64) Dynamics[1] = Dot( Bx>,  System.GetDynamics(Co) )
-> (65) Dynamics[1] = m*g*L*sin(phi) + (I-J)*wy*wz + J*(wC*wy+wx')

(66) Dynamics[2] = Dot( By>,  System.GetDynamics(Co) )
-> (67) Dynamics[2] = -(I-J)*wx*wz - J*(wC*wx-wy')

(68) Dynamics[3] = Dot( Bz>,  System.GetDynamics(Co) )
-> (69) Dynamics[3] = I*wz'

(70) Solve( Dynamics = 0,  wx', wy', wz' )
-> (71) wx' = -wC*wy - (m*g*L*sin(phi)+(I-J)*wy*wz)/J
-> (72) wy' = -wx*(wz-wC-I*wz/J)
-> (73) wz' = 0

(74) %--------------------------------------------------------------------
(75) %       System's linear/angular momentum and kinetic/potential energy.
(76) L> = System.GetLinearMomentum()
-> (77) L> = m*L*wy*Bx> - m*L*wx*By>

(78) H> = System.GetAngularMomentum( Co )
-> (79) H> = J*wx*Bx> + J*wy*By> + I*wz*Bz>

(80) KE = System.GetKineticEnergy()
-> (81) KE = 0.5*I*wz^2 + 0.5*J*wx^2 + 0.5*J*wy^2

(82) PE = System.GetForceGravityPotentialEnergy( -g*Nz>, No )
-> (83) PE = m*g*L*cos(phi)

(84) KePe = KE + PE             % Total mechanical energy
-> (85) KePe = KE + PE

(86) %--------------------------------------------------------------------
(87) Hnz = Dot( Nz>,  H> )      % Nz> measure of angular momentum (conserved).
-> (88) Hnz = I*wz*cos(phi) - J*wy*sin(phi)

(89) Hbz = Dot( Bz>,  H> )      % Bz> measure of angular momentum (conserved).
-> (90) Hbz = I*wz

(91) %--------------------------------------------------------------------
(92) %       Integration parameters and initial values.
(93) Input  tFinal = 4 sec,  tStep = 0.05 sec,  absError = 1.0E-07
(94) Input  theta  = 0 deg,      phi  = 20 deg
(95) Input  wx = 0 rad/sec,       wy  = 0 rad/sec,     wz = 300 rpm
(96) %--------------------------------------------------------------------
(97) %       List output quantities and solve ODEs.
(98) Output  t sec,  theta deg,  phi deg,  KePe Joules,  Hnz kg*m^2/sec,  Hbz kg*m^2/sec
(99) ODE()  MGGyroPrecessionNutationSpinWithWxyz

(100) %--------------------------------------------------------------------
```
Saved by Motion Genesis LLC.   Portions copyright (c) 2009-2017 Motion Genesis LLC. Rights reserved. Only for use with MotionGenesis.