MGHumanOnTurntableWithGyroFBD.html  (MotionGenesis input/output).
   (1) % MotionGenesis file:  MGHumanOnTurntableWithGyroFBD.txt
   (2) %------------------------------------------------------------
   (3) NewtonianFrame  N   
   (4) RigidBody       A   % Turntable, human legs, torso, and head.
   (5) RigidFrame      B   % Human shoulder, arms, and wheel's axle.
   (6) RigidBody       C   % Bicycle wheel (rotor).
   (7) %------------------------------------------------------------
   (8) Variable    qA''        % Az> measure of angle from Nx> to Ax>
   (9) Specified   qB''        % Ax> measure of angle from Ay> to By> 
   (10) Variable    wC'         % By> measure of C's angular velocity in B
   (11) Constant    Lz = 1.2 m  % Az> measure of Bo from No 
   (12) Constant    Lx = 0.5 m  % Ax> measure of Ccm from Bo
   (13) C.SetMass( mC = 2 kg )
   (14) A.SetInertia( Acm,  IAxx,  IAyy,  IAzz = 0.64 kg*m^2 )
   (15) C.SetInertia( Ccm, B,  IC = 0.12 kg*m^2,  JC = 0.24 kg*m^2,  IC  )
   (16) %------------------------------------------------------------
   (17) %       Rotational kinematics.
   (18) A.RotateZ( N, qA )
-> (19) A_N = [cos(qA), sin(qA), 0;  -sin(qA), cos(qA), 0;  0, 0, 1]
-> (20) w_A_N> = qA'*Az>
-> (21) alf_A_N> = qA''*Az>

   (22) B.RotateX( A, qB )
-> (23) B_A = [1, 0, 0;  0, cos(qB), sin(qB);  0, -sin(qB), cos(qB)]
-> (24) w_B_A> = qB'*Bx>
-> (25) w_B_N> = qA'*Az> + qB'*Bx>
-> (26) alf_B_A> = qB''*Bx>
-> (27) alf_B_N> = qB'*qA'*Ay> + qA''*Az> + qB''*Bx>

   (28) C.SetAngularVelocityAcceleration( B, wC*By> )
-> (29) w_C_B> = wC*By>
-> (30) w_C_N> = qB'*Bx> + (wC+sin(qB)*qA')*By> + cos(qB)*qA'*Bz>
-> (31) alf_C_B> = wC'*By>
-> (32) alf_C_N> = (qB''-cos(qB)*wC*qA')*Bx> + (wC'+cos(qB)*qB'*qA'+sin(qB)*qA
        '')*By> + (qB'*wC+cos(qB)*qA''-sin(qB)*qB'*qA')*Bz>

   (33) %------------------------------------------------------------
   (34) %       Translational kinematics.
   (35) Acm.SetVelocity( N, 0> )
-> (36) v_Acm_N> = 0>

   (37) Bo.Translate(  No, Lz*Az> )
-> (38) p_No_Bo> = Lz*Az>
-> (39) v_Bo_N> = 0>
-> (40) a_Bo_N> = 0>

   (41) Ccm.Translate( Bo, Lx*Ax> )
-> (42) p_Bo_Ccm> = Lx*Ax>
-> (43) v_Ccm_N> = Lx*qA'*Ay>
-> (44) a_Ccm_N> = -Lx*qA'^2*Ax> + Lx*qA''*Ay>

   (45) %------------------------------------------------------------
   (46) %       Form equations of motion (angular momentum principle).
   (47) Dynamics[1] = Dot( Az>,  System(A,C).GetDynamics(Bo)  )    
-> (48) Dynamics[1] = (IAzz+mC*Lx^2)*qA'' + JC*sin(qB)*(wC'+cos(qB)*qB'*qA'+sin
        (qB)*qA'') - cos(qB)*((IC-JC)*qB'*(wC+sin(qB)*qA')+IC*(sin(qB)*qB'*qA'-
        qB'*wC-cos(qB)*qA''))

   (49) Dynamics[2] = Dot( By>,            C.GetDynamics(Ccm) )    
-> (50) Dynamics[2] = JC*(wC'+cos(qB)*qB'*qA'+sin(qB)*qA'')

   (51) FactorQuadratic( Dynamics,  qA', qB', wC )
-> (52) Dynamics[1] = JC*cos(qB)*qB'*wC + JC*sin(qB)*wC' + IC*cos(qB)^2*qA'' +
        JC*sin(qB)^2*qA'' + (IAzz+mC*Lx^2)*qA'' - 2*(IC-JC)*sin(qB)*cos(qB)*qB'
        *qA'

-> (53) Dynamics[2] = JC*(wC'+cos(qB)*qB'*qA'+sin(qB)*qA'')

   (54) %------------------------------------------------------------
   (55) %       System angular momentum about vertical axis passing through Ao.
   (56) %       Note:  Its magnitude is not constant whereas Az> measure is constant.
   (57) SystemAngularMomentumAboutBo> = System.GetAngularMomentum( Bo )
-> (58) SystemAngularMomentumAboutBo> = (IAzz+mC*Lx^2)*qA'*Az> + IC*qB'*Bx> + 
        JC*(wC+sin(qB)*qA')*By> + IC*cos(qB)*qA'*Bz>

   (59) SystemAngularMomentumAboutBoZ = Dot( SystemAngularMomentumAboutBo>,  Az> )
-> (60) SystemAngularMomentumAboutBoZ = IC*cos(qB)^2*qA' + (IAzz+mC*Lx^2)*qA'
         + JC*sin(qB)*(wC+sin(qB)*qA')

   (61) SystemAngularMomentumMagnitudeAboutBo = GetMagnitude( SystemAngularMomentumAboutBo> )
-> (62) SystemAngularMomentumMagnitudeAboutBo = sqrt(IC^2*qB'^2+IC^2*cos(qB)^2*
        qA'^2+(IAzz+mC*Lx^2)^2*qA'^2+JC^2*(wC+sin(qB)*qA')^2+2*IC*(IAzz+mC*Lx^2)
        *cos(qB)^2*qA'^2+2*JC*(IAzz+mC*Lx^2)*sin(qB)*qA'*(wC+sin(qB)*qA'))

   (63) SystemAngularMomentumMagnitudeAboutNo = GetMagnitude( System.GetAngularMomentum(No) )
-> (64) SystemAngularMomentumMagnitudeAboutNo = sqrt(IC^2*cos(qB)^2*qA'^2+(IAzz
        +mC*Lx^2)^2*qA'^2+JC^2*(wC+sin(qB)*qA')^2+(IC*qB'-mC*Lx*Lz*qA')^2+2*IC*
        (IAzz+mC*Lx^2)*cos(qB)^2*qA'^2+2*JC*(IAzz+mC*Lx^2)*sin(qB)*qA'*(wC+sin(
        qB)*qA'))

   (65) %------------------------------------------------------------
   (66) %       Bicycle C's angular momentum about vertical axis passing through Ccm.
   (67) %       Note:  Its magnitude is not constant whereas By measure is constant.
   (68) CAngularMomentumAboutCcm> = C.GetAngularMomentum( Ccm )
-> (69) CAngularMomentumAboutCcm> = IC*qB'*Bx> + JC*(wC+sin(qB)*qA')*By> + IC*
        cos(qB)*qA'*Bz>

   (70) CAngularMomentumMagnitudeAboutCcm = GetMagnitude( CAngularMomentumAboutCcm> )
-> (71) CAngularMomentumMagnitudeAboutCcm = sqrt(IC^2*qB'^2+IC^2*cos(qB)^2*qA'^
        2+JC^2*(wC+sin(qB)*qA')^2)

   (72) CAngularMomentumAboutCcmY = Dot( CAngularMomentumAboutCcm>, By> )
-> (73) CAngularMomentumAboutCcmY = JC*(wC+sin(qB)*qA')

   (74) %------------------------------------------------------------
   (75) %       Sum of kinetic and potential energy (not constant).
   (76) SystemKineticEnergy = System.GetKineticEnergy()      
-> (77) SystemKineticEnergy = 0.5*IC*qB'^2 + 0.5*IAzz*qA'^2 + 0.5*mC*Lx^2*qA'^2
        + 0.5*IC*cos(qB)^2*qA'^2 + 0.5*JC*(wC+sin(qB)*qA')^2

   (78) FactorQuadratic( SystemKineticEnergy, qA', qB', wC )
-> (79) SystemKineticEnergy = 0.5*IC*qB'^2 + 0.5*JC*wC^2 + JC*sin(qB)*wC*qA' +
        0.5*(IAzz+IC+mC*Lx^2-(IC-JC)*sin(qB)^2)*qA'^2

   (80) %------------------------------------------------------------
   (81) %       Integration parameters and initial values.
   (82) Input  tFinal = 8 sec,  tStep = 0.02 sec,  absError = 1.0E-08
   (83) Input  qA = 0 deg,  qA' = 0 rad/sec,  wC = 40 rad/sec
   (84) %------------------------------------------------------------
   (85) %       Specified expressions 
   (86) SetDt( qB = pi/6 * sin( pi/2*t ) )
-> (87) qB = 0.5235988*sin(1.570796*t)
-> (88) qB' = 0.822467*cos(1.570796*t)
-> (89) qB'' = -1.291928*sin(1.570796*t)

   (90) %------------------------------------------------------------
   (91) %       Output quantities by ODE command.
   (92) Output t sec,  qA deg,  qA' deg/sec,  qB deg,  wC rad/sec,  &
       SystemAngularMomentumMagnitudeAboutNo kg*m^2/sec,    &
       SystemAngularMomentumMagnitudeAboutBo kg*m^2/sec,    &
       SystemAngularMomentumAboutBoZ         kg*m^2/sec,    &
       CAngularMomentumMagnitudeAboutCcm     kg*m^2/sec,    &
       CAngularMomentumAboutCcmY             kg*m^2/sec,    &
       SystemKineticEnergy                   Joules
   (93) %------------------------------------------------------------
   (94) %       Form and solve ODEs.
   (95) Solve( Dynamics = 0,   qA'', wC' )  
-> (96) qA'' = cos(qB)*qB'*(2*IC*sin(qB)*qA'-JC*wC-JC*sin(qB)*qA')/(IAzz+mC*Lx^
        2+IC*cos(qB)^2)
-> (97) wC' = -cos(qB)*qB'*(qA'-sin(qB)*(JC*wC-(2*IC-JC)*sin(qB)*qA')/(IAzz+mC*
        Lx^2+IC*cos(qB)^2))

   (98) ODE()  MGHumanOnTurntableWithGyroFBD

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