%    File: MGSkyDiverParachuteDynamics.txt
% Purpose: Simulation of a sky-diver with free-fall then parachute deployment.
%    Note: Air-resistance forces change suddenly due to opening parachute.
%--------------------------------------------------------------------
NewtonianFrame N        % Earth with Ny> vertically upward.
Particle       Q        % Parachutist.
%--------------------------------------------------------------------
Variable v'                  % Downward measure of Q's velocity relative to Earth.
Constant vTerminal = 8 m/s   % Terminal velocity with open parachute (28.8 km/hr, 17.9 mph).
Constant g = 9.8 m/s^2       % Earth's gravitational acceleration.
Constant b                   % Air-resistance drag constant.
Q.SetMass( m = 100 kg )
%--------------------------------------------------------------------
%   Set Q's velocity and acceleration relative to Earth.
Q.SetVelocityAcceleration( N,  -v*Ny> )
%--------------------------------------------------------------------
%   Add air-resistance and gravity forces on parachutist.
Fair> = -b * Q.GetVelocity(N)
Q.AddForce( -m*g*Ny> + Fair> )
%--------------------------------------------------------------------
%   Form statics equation which is relevant for terminal velocity.
Statics = Evaluate( Dot( Q.GetStatics(), Ny> ),  v = vTerminal )
Solve( Statics = 0,  b )
%--------------------------------------------------------------------
%   Create an expression for b that makes b = 0 for the first 4 seconds (free-fall)
%   and thereafter ensure b is the value determined from terminal velocity.
isTimeGreaterThan4 = IsPositive( t - 4 )
b *= isTimeGreaterThan4
%--------------------------------------------------------------------
%   Form dynamic equations with F = m a.
Dynamics = Dot( Q.GetDynamics(),  Ny> )
Solve( Dynamics = 0,   v' )
%--------------------------------------------------------------------
%   Height y is measured positive upward.
Variable y' = -v
%--------------------------------------------------------------------
%   Initial values for variables -- for subsequent ODE() command.
Input   y = 200 m,  v = 0 m/s
%--------------------------------------------------------------------
%   List output quantities from the ODE() command.
Output  t sec,  y m,  v m/s,  v' m/s,  Magnitude( Fair> ) Newtons
%--------------------------------------------------------------------
%   Set numerical integration parameters and solve the ODEs for y(t), v(t).
Input  tFinal = 9 sec,  tStep = 0.02 sec
ODE()  MGSkyDiverParachuteDynamics
%--------------------------------------------------------------------
%   Optional: Plot results for v(t) vs. t.
% Plot MGSkyDiverParachuteDynamics.1[ 1, 3 ]
%--------------------------------------------------------------------
%   Optional: Verify dynamics with Kane's method.
SetGeneralizedSpeed( v )
KaneDynamics = System.GetDynamicsKane()
isSameDynamics = IsSimplifyEqual( KaneDynamics, Dynamics )
%--------------------------------------------------------------------
%   Record input together with responses.
Save  MGSkyDiverParachuteDynamics.html
Quit