MGGyroPrecessionNutationSpinWithWxyz.html  (MotionGenesis input/output).
   (1) % MotionGenesis file:  MGGyroPrecessionNutationSpinWithWxyz.txt
   (2) % Copyright (c) 2009 Motion Genesis LLC.  All rights reserved.
   (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) Zero> = wx*Bx> + wy*By> + wz*Bz> - C.GetAngularVelocity(N)
-> (35) Zero> = (wx+phi')*Bx> + (wy-sin(phi)*theta')*By> + (wz+cos(phi)*theta'-
        wC)*Bz>

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

   (38) Solve( ChangeVariablesEqn = 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).
   (60) CCm.AddForce( -m*g*Nz> )
-> (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.