%filename: freefall.m
% ------> Written 11/9/05 by JE Hertel
%This program simulates jumping out of an airplane by integrating
%           acceleration twice. (Uses simple rectangular integration.)
%Note: would be nice to add a parachute phase to soften the landing.
%
clc,clear,clf

%Initialize variables:
t(1)=0;                     %(sec) time of jump
xvel(1)=125*5280/60/60;     %(ft/sec) airplane is flying straight and level at 125 mph
yvel(1)=0;
xpos(1)=0;                  %(ft) position at jump time
ypos(1)=2000;               %(ft) altitude at jump time
wt=180;                     %(lbf) weight of jumper
mass=wt/32.2;               %(slug)mass of jumper
drag_coef=1.2;
frt_area=6;                 %(ft^2)frontal area
air_den=.00242236;          %(slug/ft^3) air density

%Setup Integration
i=1;
dt=.1;                       %timestep (sec)

while ypos(i)>0              %loop til you hit the ground
    
    %calculate forces
    tot_vel=(xvel(i)^2+yvel(i)^2)^.5;   
    drag_force=.5*drag_coef*frt_area*air_den*tot_vel^2;
    xforce=-xvel(i)/tot_vel*drag_force;
    yforce=-yvel(i)/tot_vel*drag_force-wt;
    
    i=i+1;                      % use i for a counter/pointer
    
    xacc=xforce/mass;           %(ft/sec^2) x-dir acceleration
    yacc=yforce/mass;           %(ft/sec^2) y-dir acceleration
    
    xvel(i)=               %integrate x-acceleration to get x-velocity
    yvel(i)=                %integrate y-acceleration to get y-velocity
       
    xpos(i)=                %integrate x-velocity to get x-position
    ypos(i)=                %integrate y-velocity to get y-position
    
    t(i)=                   %increment time

end

%create plots
subplot(3,1,1)          %divide figure window into three rows --- plot top area
plot(t,xvel,t,yvel,'--')
title('Figure 1:  Velocity of a Skydiver in Freefall')
xlabel('Time (sec)')
ylabel('X and Y velocity (ft/sec)')
legend('X velocity', 'Y velocity')

subplot(3,1,2)          %plot middle area
plot(t,xpos,t,ypos,'--')
title('Figure 2:  Position of a Skydiver in Freefall')
xlabel('Time (sec)')
ylabel('X and Y position (ft)')
legend('X position', 'Y position')
yaxis_lim=ylim; 
ylim([0, yaxis_lim(2)]);  %set y axis minimum = 0

subplot(3,1,3)          %plot bottom area
plot(xpos,ypos)
title('Figure 3:  Trajectory of a Skydiver in Freefall')
xlabel('X position (ft)')
ylabel('Y position (ft)')
yaxis_lim=ylim; 
ylim([0, yaxis_lim(2)]);  %set y axis minimum = 0