MGWrapParticleAroundHorizontalCylinder2D.html  (MotionGenesis input/output).
   (1) % MotionGenesis file:  MGWrapParticleAroundHorizontalCylinder2D.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      B                % Taut string with By> directed towards Q.
   (6) Particle        Q                % Particle at end of string.
   (7) Point           BN( B )          % Point of string B in contact with cylinder N.
   (8) %--------------------------------------------------------------------
   (9) Variable    theta''              % Angle from Ny> to By> with sense +Nz>.
   (10) Variable    Tension              % Tension in string.
   (11) Variable    L''                  % Length of deployed string (between BN and Q).
   (12) Constant    R = 3 cm             % Radius of cylinder N.
   (13) Constant    g = 9.8 m/s^2        % Earth's gravitational acceleration.
   (14) Q.SetMass( m = 100 g )
   (15) %--------------------------------------------------------------------
   (16) %   Rotational and translational kinematics.
   (17) B.RotateZ( N, theta )
-> (18) B_N = [cos(theta), sin(theta), 0;  -sin(theta), cos(theta), 0;  0, 0, 1]
-> (19) w_B_N> = theta'*Bz>
-> (20) alf_B_N> = theta''*Bz>

   (21) Q.Translate( No, R*Bx> + L*By> ) 
-> (22) p_No_Q> = R*Bx> + L*By>
-> (23) v_Q_N> = -L*theta'*Bx> + (L'+R*theta')*By>
-> (24) a_Q_N> = (-2*L'*theta'-R*theta'^2-L*theta'')*Bx> + (L''+R*theta''-L*theta'^2)*By>

   (25) BN.SetPositionVelocity( Q,  -L*By>, B )
-> (26) p_Q_BN> = -L*By>
-> (27) v_BN_N> = (L'+R*theta')*By>

   (28) %--------------------------------------------------------------------
   (29) %   Constraint: Total length of string does not change.
   (30) %   Alternate:  Dot( v_BN_N>, By> ) = 0.
   (31) %   Let  s = Additional length of string wrapped around cylinder after t = 0.
   (32) %        s = R * (theta - theta0) =  L0 - L   
   (33) Constant    L0 = 70 cm           % Initial value of L.
   (34) Input       theta = 0 deg        % Initial value of theta.
   (35) s = R * (theta - Input(theta))
-> (36) s = R*theta

   (37) SolveDt( s = L0 - L,   L )       % Find L', L'' in terms of theta', theta''.
-> (38) L = L0 - s
-> (39) L' = -R*theta'
-> (40) L'' = -R*theta''

   (41) %--------------------------------------------------------------------
   (42) %   Q's velocity/acceleration simplifies by substituting L', L''.
   (43) Explicit( v_Q_N>,  theta' )
-> (44) v_Q_N> = -L*theta'*Bx>

   (45) Explicit( a_Q_N>,  theta', theta'' )
-> (46) a_Q_N> = (R*theta'^2-L*theta'')*Bx> - L*theta'^2*By>

   (47) %--------------------------------------------------------------------
   (48) %   Add relevant contact/distance forces.
   (49) Q.AddForce(  -Tension*By> - m*g*Ny> )
-> (50) Force_Q> = -Tension*By> - m*g*Ny>

   (51) %--------------------------------------------------------------------
   (52) %   Form equations of motion via F = m*a and solve for theta'', Tension.
   (53) Dynamics[1] = Dot(  Q.GetDynamics(),  Bx>  )
-> (54) Dynamics[1] = m*(g*sin(theta)+R*theta'^2-L*theta'')

   (55) Dynamics[2] = Dot(  Q.GetDynamics(),  By>  )
-> (56) Dynamics[2] = Tension + m*g*cos(theta) - m*L*theta'^2

   (57) Solve( Dynamics = 0,   theta'',  Tension )
-> (58) theta'' = (g*sin(theta)+R*theta'^2)/L
-> (59) Tension = -m*(g*cos(theta)-L*theta'^2)

   (60) %--------------------------------------------------------------------
   (61) %   Use initial velocity to solve for for initial value of theta'.
   (62) Constant  v0 =  -4 m/s          % Initial Bx> measure of Q's velocity.
   (63) vx = Dot( Q.GetVelocity(N), Bx> )
-> (64) vx = -L*theta'

   (65) SolveSetInput( vx = v0,   theta' = 0 rad/sec )

->   %  INPUT has been assigned as follows:
->   %   theta'                    5.714285714285714       rad/sec

   (66) %--------------------------------------------------------------------
   (67) %   Calculations for plotting Q's trajectory in N.
   (68) x = Dot(  Q.GetPosition(No),  Nx>  )
-> (69) x = R*cos(theta) - L*sin(theta)

   (70) y = Dot(  Q.GetPosition(No),  Ny>  )
-> (71) y = R*sin(theta) + L*cos(theta)

   (72) %--------------------------------------------------------------------
   (73) %   Calculations for plotting angular momentum and kinetic + potential energy.
   (74) KineticEnergy = Q.GetKineticEnergy()
-> (75) KineticEnergy = 0.5*m*L^2*theta'^2

   (76) GravityPotentialEnergy = Q.GetForceGravityPotentialEnergy( -g*Ny>, No )
-> (77) GravityPotentialEnergy = m*g*(R*sin(theta)+L*cos(theta))

   (78) MechanicalEnergy = KineticEnergy + GravityPotentialEnergy
-> (79) MechanicalEnergy = GravityPotentialEnergy + KineticEnergy

   (80) HNo = Dot(  Q.GetAngularMomentum( No ),  Nz>  )
-> (81) HNo = m*L^2*theta'

   (82) HBN = Dot(  Q.GetAngularMomentum( BN ),  Nz>  )
-> (83) HBN = m*L^2*theta'

   (84) vSpeed = Q.GetSpeed( N )
-> (85) vSpeed = abs(L*theta')

   (86) %--------------------------------------------------------------------
   (87) %   List output quantities, provide integration parameters, and solve ODEs.
   (88) Output  t sec,  x cm,  y cm,  L cm, theta deg, Tension N, MechanicalEnergy Joules, vSpeed cm/s, HNo kg*m^2/s^2, HBN kg*m^2/s^2
   (89) Input  tFinal = 1.536 sec,  tStep = 0.002 sec,  absError = 1.0E-07
   (90) ODE()  MGWrapParticleAroundHorizontalCylinder2D

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