(1) %    File: MGSkyDiverParachuteDynamics.txt
   (2) % Purpose: Simulation of a sky-diver with free-fall then parachute deployment.
   (3) %    Note: Air-resistance forces change suddenly due to opening parachute.
   (4) %--------------------------------------------------------------------
   (5) NewtonianFrame N        % Earth with Ny> vertically upward.
   (6) Particle       Q        % Parachutist.
   (7) %--------------------------------------------------------------------
   (8) Variable v'                  % Downward measure of Q's velocity relative to Earth.
   (9) Constant vTerminal = 8 m/s   % Terminal velocity with open parachute (28.8 km/hr, 17.9 mph).
   (10) Constant g = 9.8 m/s^2       % Earth's gravitational acceleration.
   (11) Constant b                   % Air-resistance drag constant.
   (12) Q.SetMass( m = 100 kg )
   (13) %--------------------------------------------------------------------
   (14) %   Set Q's velocity and acceleration relative to Earth.
   (15) Q.SetVelocityAcceleration( N,  -v*Ny> )
-> (16) v_Q_N> = -v*Ny>
-> (17) a_Q_N> = -v'*Ny>

   (18) %--------------------------------------------------------------------
   (19) %   Add air-resistance and gravity forces on parachutist.
   (20) Fair> = -b * Q.GetVelocity(N)
-> (21) Fair> = b*v*Ny>

   (22) Q.AddForce( -m*g*Ny> + Fair> )
-> (23) Force_Q> = (b*v-m*g)*Ny>

   (24) %--------------------------------------------------------------------
   (25) %   Form statics equation which is relevant for terminal velocity.
   (26) Statics = Evaluate( Dot( Q.GetStatics(), Ny> ),  v = vTerminal )
-> (27) Statics = b*vTerminal - m*g

   (28) Solve( Statics = 0,  b )
-> (29) b = m*g/vTerminal

   (30) %--------------------------------------------------------------------
   (31) %   Create an expression for b that makes b = 0 for the first 4 seconds (free-fall)
   (32) %   and thereafter ensure b is the value determined from terminal velocity.
   (33) isTimeGreaterThan4 = IsPositive( t - 4 )
-> (34) isTimeGreaterThan4 = IsPositive(-4+t)

   (35) b *= isTimeGreaterThan4
-> (36) b = m*g*isTimeGreaterThan4/vTerminal

   (37) %--------------------------------------------------------------------
   (38) %   Form dynamic equations with F = m a.
   (39) Dynamics = Dot( Q.GetDynamics(),  Ny> )
-> (40) Dynamics = m*g - b*v - m*v'

   (41) Solve( Dynamics = 0,   v' )
-> (42) v' = g - b*v/m

   (43) %--------------------------------------------------------------------
   (44) %   Height y is measured positive upward.
   (45) Variable y' = -v
-> (46) y' = -v

   (47) %--------------------------------------------------------------------
   (48) %   Initial values for variables -- for subsequent ODE() command.
   (49) Input   y = 200 m,  v = 0 m/s
   (50) %--------------------------------------------------------------------
   (51) %   List output quantities from the ODE() command.
   (52) Output  t sec,  y m,  v m/s,  v' m/s,  Magnitude( Fair> ) Newtons
   (53) %--------------------------------------------------------------------
   (54) %   Set numerical integration parameters and solve the ODEs for y(t), v(t).
   (55) Input  tFinal = 9 sec,  tStep = 0.02 sec
   (56) ODE()  MGSkyDiverParachuteDynamics

   (57) %--------------------------------------------------------------------
   (58) %   Optional: Plot results for v(t) vs. t.
   (59) % Plot MGSkyDiverParachuteDynamics.1[ 1, 3 ]
   (60) %--------------------------------------------------------------------
   (61) %   Optional: Verify dynamics with Kane's method.
   (62) SetGeneralizedSpeed( v )
   (63) KaneDynamics = System.GetDynamicsKane()
-> (64) KaneDynamics = [b*v + m*v' - m*g]

   (65) isSameDynamics = IsSimplifyEqual( KaneDynamics, Dynamics )
-> (66) isSameDynamics = true

   (67) %--------------------------------------------------------------------
   (68) %   Record input together with responses.