MGWrapParticleAroundVerticalCylinder3D.html  (MotionGenesis input/output).
   (1) % MotionGenesis file:  MGWrapParticleAroundVerticalCylinder3D.txt
   (2) % Copyright (c) 2016 Motion Genesis LLC.  Only for use with MotionGenesis.
   (3) %--------------------------------------------------------------------
   (4) NewtonianFrame  N                % Cylinder with Nz> along symmetry axis.
   (5) RigidFrame      A                % Associated with theta rotation about +Nz>
   (6) RigidFrame      B                % Taut string with By> directed towards Q.
   (7) Particle        Q                % Particle at end of string.
   (8) Point           BN( B )          % Point of string B in contact with cylinder N.
   (9) Point           P                % Path point:  Traces out wrap on cylinder.
   (10) %--------------------------------------------------------------------
   (11) Variable    qA''                 % Angle from Nx> to Ax> with sense +Nz>.
   (12) Variable    qB''                 % Angle from Ay> to By> with sense +Ax>.
   (13) Variable     z''                 % Nz> measure of P's position from No.
   (14) Variable    L''                  % Length of deployed string (between P and Q).
   (15) Variable    Tension              % Tension in string.
   (16) Constant    R = 3 cm             % Radius of cylinder N.
   (17) Constant    g = 9.8 m/s^2        % Earth's gravitational acceleration.
   (18) Q.SetMass( m = 100 g )
   (19) %--------------------------------------------------------------------
   (20) %   Rotational and translational kinematics.
   (21) A.RotateZ( N, qA )
-> (22) A_N = [cos(qA), sin(qA), 0;  -sin(qA), cos(qA), 0;  0, 0, 1]
-> (23) w_A_N> = qA'*Az>
-> (24) alf_A_N> = qA''*Az>

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

   (31) P.Translate( No, R*Ax> + z*Az> ) 
-> (32) p_No_P> = R*Ax> + z*Az>
-> (33) v_P_N> = R*qA'*Ay> + z'*Az>
-> (34) a_P_N> = -R*qA'^2*Ax> + R*qA''*Ay> + z''*Az>

   (35) Q.Translate( P,  L*By> ) 
-> (36) p_P_Q> = L*By>
-> (37) v_Q_N> = R*qA'*Ay> + z'*Az> - L*cos(qB)*qA'*Bx> + L'*By> + L*qB'*Bz>
-> (38) a_Q_N> = (2*L*sin(qB)*qA'*qB'-R*qA'^2-2*cos(qB)*L'*qA'-L*cos(qB)*qA'')*Ax>
        + (R*qA''-L*cos(qB)*qA'^2)*Ay> + z''*Az> + (L''-L*qB'^2)*By> + (2*L'*
        qB'+L*qB'')*Bz>

   (39) BN.SetPositionVelocity( Q,  -L*By>, B )
-> (40) p_Q_BN> = -L*By>
-> (41) v_BN_N> = R*qA'*Ay> + z'*Az> + L'*By>

   (42) %--------------------------------------------------------------------
   (43) %   Constraint 1: Total length of string does not change.
   (44) %   Constraint 2: Path point P only moves in By> direction (no sliding in Bx> direction).
   (45) %   Alternate:    v_BN_N> = 0>.
   (46) %   Let  s = Additional length of string wrapped around cylinder after t = 0.
   (47) %        s = FindByIntegration =  L0 - L  
   (48) Variable s' = Dot( P.GetVelocity(N), By> )
-> (49) s' = sin(qB)*z' + R*cos(qB)*qA'

   (50) Input    s  = 0 cm                     % Initial value of s.
   (51) Constant    L0 = 70 cm                 % Initial value  of L.
   (52) SolveDt( s = L0 - L,   L )             % Find L', L'' in terms of theta', theta''.
-> (53) L = L0 - s
-> (54) L' = -s'
-> (55) L'' = R*sin(qB)*qA'*qB' - cos(qB)*qB'*z' - sin(qB)*z'' - R*cos(qB)*qA''

   (56) SolveDt( Dot( P.GetVelocity(N), Bz> ),  z' )
-> (57) z' = R*tan(qB)*qA'
-> (58) z'' = (sin(qB)*qB'*z'+R*(cos(qB)*qA'*qB'+sin(qB)*qA''))/cos(qB)

   (59) %--------------------------------------------------------------------
   (60) %   Add relevant contact/distance forces.
   (61) Q.AddForce(  -Tension*By> - m*g*Nz> )
-> (62) Force_Q> = -Tension*By> - m*g*Nz>

   (63) %--------------------------------------------------------------------
   (64) %   Form equations of motion via F = m*a and solve for qA'', qB'', Tension.
   (65) Dynamics[1] = Dot(  Q.GetDynamics(),  Bx>  )
-> (66) Dynamics[1] = m*(2*L*sin(qB)*qA'*qB'-R*qA'^2-2*cos(qB)*L'*qA'-L*cos(qB)
        *qA'')

   (67) Dynamics[2] = Dot(  Q.GetDynamics(),  By>  )
-> (68) Dynamics[2] = Tension + m*sin(qB)*(g+z'') - m*(L*qB'^2-L'') - m*cos(qB)
        *(L*cos(qB)*qA'^2-R*qA'')

   (69) Dynamics[3] = Dot(  Q.GetDynamics(),  Bz>  )
-> (70) Dynamics[3] = m*(2*L'*qB'+L*qB''+cos(qB)*(g+z'')+sin(qB)*(L*cos(qB)*qA'^2
        -R*qA''))

   (71) Solve( Dynamics = 0,   qA'',  qB'', Tension )
-> (72) qA'' = -qA'*(R*qA'+2*cos(qB)*L'-2*L*sin(qB)*qB')/(L*cos(qB))
-> (73) qB'' = -sin(qB)*cos(qB)*qA'^2 - (g*cos(qB)+2*L'*qB'+sin(qB)*qB'*z'+R*
        cos(qB)*qA'*qB')/L
-> (74) Tension = m*(L*cos(qB)^2*qA'^2+qB'*(L*qB'+z'/cos(qB))-sin(qB)*(g+qB'*(R
        *qA'+tan(qB)*z')))

   (75) %--------------------------------------------------------------------
   (76) %   Calculations for plotting Q's trajectory in N and solving initial values.
   (77) xQ = Dot(  Q.GetPosition(No),  Nx>  )
-> (78) xQ = R*cos(qA) - L*sin(qA)*cos(qB)

   (79) yQ = Dot(  Q.GetPosition(No),  Ny>  )
-> (80) yQ = R*sin(qA) + L*cos(qA)*cos(qB)

   (81) zQ = Dot(  Q.GetPosition(No),  Nz>  )
-> (82) zQ = z + L*sin(qB)

   (83) %--------------------------------------------------------------------
   (84) %   Use Q's initial position to solve for initial values of qA, qB, z.
   (85) SolveSetInput(  [xQ;  yQ;  zQ] = [R;  L0;  0],  qA =  0 deg,  qB =  0 deg,  z =  0 cm )

->   %  INPUT has been assigned as follows:
->   %   qA                        0                       deg
->   %   qB                        0                       deg
->   %   z                         0                       cm

   (86) %--------------------------------------------------------------------
   (87) %   Use Q's initial velocity to solve for initial value of qA', qB'.
   (88) %   Note: Already used Bz> measure of P's velocity in N to solve for z'.
   (89) vx = Explicit( Dot( Q.GetVelocity(N), Ax> ),  qA', qB' )
-> (90) vx = -L*cos(qB)*qA'

   (91) vy = Explicit( Dot( Q.GetVelocity(N), Ay> ),  qA', qB'  )
-> (92) vy = -L*sin(qB)*qB'

   (93) vz = Explicit( Dot( Q.GetVelocity(N), Az> ),  qA', qB'  )
-> (94) vz = L*cos(qB)*qB'

   (95) SolveSetInput( [vx; vz] = [-4;  0],   qA' = 4 rad/sec,  qB' = 0 rad/sec )

->   %  INPUT has been assigned as follows:
->   %   qA'                       5.714285714285714       rad/sec
->   %   qB'                       0                       rad/sec

   (96) %--------------------------------------------------------------------
   (97) %   Calculations for plotting angular momentum and kinetic + potential energy.
   (98) KineticEnergy = Q.GetKineticEnergy()
-> (99) KineticEnergy = -0.5*m*(2*R*L*sin(qB)*qA'*qB'-L'^2-z'^2-R^2*qA'^2-L^2*qB'^2
        -2*sin(qB)*L'*z'-2*R*cos(qB)*L'*qA'-2*L*cos(qB)*qB'*z'-L^2*cos(qB)^2*qA'^2)

   (100) GravityPotentialEnergy = Q.GetForceGravityPotentialEnergy( -g*Nz>, No )
-> (101) GravityPotentialEnergy = m*g*(z+L*sin(qB))

   (102) MechanicalEnergy = KineticEnergy + GravityPotentialEnergy
-> (103) MechanicalEnergy = GravityPotentialEnergy + KineticEnergy

   (104) HNo = Dot(  Q.GetAngularMomentum( No ),  Nz>  )
-> (105) HNo = -m*(R*L*sin(qB)*qB'-R^2*qA'-cos(qB)*(R*L'+L^2*cos(qB)*qA'))

   (106) HBN = Dot(  Q.GetAngularMomentum( BN ),  Nz>  )
-> (107) HBN = m*L^2*cos(qB)^2*qA'

   (108) vSpeed = Q.GetSpeed( N )
-> (109) vSpeed = sqrt(L'^2+z'^2+R^2*qA'^2+L^2*qB'^2+2*sin(qB)*L'*z'+2*R*cos(
         qB)*L'*qA'+2*L*cos(qB)*qB'*z'+L^2*cos(qB)^2*qA'^2-2*R*L*sin(qB)*qA'*
         qB')

   (110) %--------------------------------------------------------------------
   (111) %   List output quantities, provide integration parameters, and solve ODEs.
   (112) Output  t sec,  xQ cm,  yQ cm, zQ cm,  L cm, qA deg, qB deg, z m, Tension N, MechanicalEnergy Joules, vSpeed cm/s, HNo kg*m^2/s^2, HBN kg*m^2/s^2
   (113) Input  tFinal = 1.7259 sec,  tStep = 0.002 sec,  absError = 1.0E-07
   (114) ODE()  MGWrapParticleAroundVerticalCylinder3D

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