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

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

   (29) C.SetAngularVelocityAcceleration( B, wC*By> )
-> (30) w_C_B> = wC*By>
-> (31) w_C_N> = qB'*Bx> + (wC+sin(qB)*qA')*By> + cos(qB)*qA'*Bz>
-> (32) alf_C_B> = wC'*By>
-> (33) 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>

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

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

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

   (46) %------------------------------------------------------------
   (47) %   Equations of motion with angular momentum principle (MG road-maps).
   (48) Dynamics[1] = Dot( Az>,  System(A,C).GetDynamics(Bo)  )
-> (49) 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''))

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

   (52) FactorQuadratic( Dynamics,  qA', qB', wC )
-> (53) 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'

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

   (55) %------------------------------------------------------------
   (56) %   Optional: Equations of motion with Kane's method.
   (57) SetGeneralizedSpeed( qA', wC )
   (58) KaneDynamics = System.GetDynamicsKane()
-> (59) KaneDynamics[1] = cos(qB)*qB'*(JC*sin(qB)*qA'+IC*(wC-sin(qB)*qA')-(IC-
        JC)*(wC+sin(qB)*qA')) + JC*sin(qB)*wC' + (IAzz+mC*Lx^2+IC*cos(qB)^2+JC*
        sin(qB)^2)*qA''

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

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

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

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

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

   (71) SystemAngularMomentumMagnitudeAboutBo = GetMagnitude( SystemAngularMomentumAboutBo> )
-> (72) 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'))

   (73) SystemAngularMomentumMagnitudeAboutNo = GetMagnitude( System.GetAngularMomentum(No) )
-> (74) 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'))

   (75) %------------------------------------------------------------
   (76) %   Wheel C's angular momentum about vertical axis passing through Ccm.
   (77) %   Note:  Its magnitude is not constant whereas By measure is constant.
   (78) CAngularMomentumAboutCcm> = C.GetAngularMomentum( Ccm )
-> (79) CAngularMomentumAboutCcm> = IC*qB'*Bx> + JC*(wC+sin(qB)*qA')*By> + IC*
        cos(qB)*qA'*Bz>

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

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

   (84) %------------------------------------------------------------
   (85) %   Sum of kinetic and potential energy (not constant).
   (86) SystemKineticEnergy = System.GetKineticEnergy()
-> (87) 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

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

   (90) %------------------------------------------------------------
   (91) %   Integration parameters and initial values.
   (92) Input  tFinal = 8 sec,  tStep = 0.02 sec,  absError = 1.0E-08
   (93) Input  qA = 0 deg,  qA' = 0 rad/sec,  wC = 40 rad/sec
   (94) %------------------------------------------------------------
   (95) %   Specified expressions
   (96) SetDt( qB = pi/6 * sin( pi/2*t ) )
-> (97) qB = 0.5235988*sin(1.570796*t)
-> (98) qB' = 0.822467*cos(1.570796*t)
-> (99) qB'' = -1.291928*sin(1.570796*t)

   (100) %--------------------------------------------------------------------
   (101) %   Solve equations of motion for qA'', wC'.
   (102) Solve( Dynamics = 0,   qA'', wC' )
-> (103) qA'' = cos(qB)*qB'*(2*IC*sin(qB)*qA'-JC*wC-JC*sin(qB)*qA')/(IAzz+mC*Lx^2
         +IC*cos(qB)^2)
-> (104) wC' = -cos(qB)*qB'*(qA'-sin(qB)*(JC*wC-(2*IC-JC)*sin(qB)*qA')/(IAzz+
         mC*Lx^2+IC*cos(qB)^2))

   (105) %------------------------------------------------------------
   (106) %   List output quantities and solve ODEs.
   (107) 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
   (108) ODE() MGHumanOnTurntableWithGyroDynamics

   (109) %------------------------------------------------------------
Saved by Motion Genesis LLC.   Command names and syntax: Copyright (c) 2009-2021 Motion Genesis LLC. All rights reserved.