% Modeling the velocity of a skydiver from free-fall through deployment
%    of the parachute.  We take
%    v' = -32.2 + 0.00104 v^2 (before parachute)
%    v' = -32.2 - 2.01 v (after)
% and v(0) = 0.  The solution is graphed and the Euler approximation
% shown.
%
% by Justin Dunmyre
% excerpted and modified by Gavin LaRose, 9/5/12

% (c)2012 Regents of the University of Michigan
% This work is licensed under the
% Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons,
% 444 Castro Street, Suite 900,
% Mountain View, CA, 94041, USA.

function SkydiverEulerBasic()

%-------------------------------------------------------------------------------
% configuration variables

% number of Euler steps (good values give a stepsize less than 1; for smaller
%    values the error is large enough that the Euler values don't appear
%    in the graph window) % '
    numEulerSteps = 50;
% time at which parachute is opened
    timeOpenParachute = 15;
% total time to which to graph the solution
    totalTime = 30;

% the time at which we show the local error (best is shortly after the
%    parachute deploys)
    localErrorTime = 16;
    
% colors
    closedColor='c';
    openColor='g';
    EulerColor='w';
    annotationColor='k';
    EulerErrorColor='k';
    
%-------------------------------------------------------------------------------
% the solutions to the differential equations; the following is all taken
%    from the non-Euler demonstration

% the solution for the velocity pre-deployment
    function v=ClosedExactSol(t)
        v = 176*(1 - exp(.366*t))./(1 + exp(.366*t));
    end

% the solution for the velocity post-deployment
    function v=OpenExactSol(t)
       % Initial Velocity based on closed parachute solution, adjusted
       initVelocity=ClosedExactSol(timeOpenParachute);
       v = -16 + (initVelocity+16)*exp(-2.01*(t-timeOpenParachute));
    end

% the differential equations, before deployment
    function vprime=skydiverClosedODE(t,v)
        vprime = -32.2 + 0.00104*v.^2;
    end
% and after deployment
    function vprime=skydiverOpenODE(t,v)
        vprime = -32.2 - 2.01*v;
    end

% the times at which we graph the velocities and corresponding
%    values of the velocities
    tClosed=linspace(0,timeOpenParachute,100);
    tOpen=linspace(timeOpenParachute,totalTime,100);
    vClosed=ClosedExactSol(tClosed);
    vOpen=OpenExactSol(tOpen);

% graph these
    scrsz = get(0,'ScreenSize'); % Get the screen size
    figh  = figure('Position',[1 scrsz(4)/4 scrsz(3)/4 scrsz(4)/2]);
    axes('Parent', figh, 'Color', [0.6 0.6 0.6]);
    hold on;
    axis([0 totalTime min(vOpen)-20 20]);
    xlabel('t','FontSize',20);
    ylabel('v','FontSize',20);
    title('Skydiver velocity vs t','FontSize',20);

% a line at zero
    plot( [0, totalTime],[0, 0],'--k' );
    
    eqntext = annotation( figh, 'textbox', [0.21 0.66 0.33 0.092], ...
                          'String', ...
                          {'v\prime=-32.2 + 0.00104 v^2', ...
                           '(before parachute)'}, ...
                          'FontSize', 20, 'Color',closedColor, ...
                          'FitBoxToText', 'on' );
    eqntext = annotation( figh, 'textbox', [0.57 0.32 0.28 0.092], ...
                          'String', ...
                          {'v\prime=-32.2 - 2.01 v', ...
                           '(after parachute)'}, ...
                          'FontSize', 20, 'Color', openColor,...
                          'FitBoxToText','on' );
    
    plot(tClosed, vClosed, [closedColor '-'], 'LineWidth',2);
    plot(tOpen, vOpen, [openColor '-'], 'LineWidth',2);

   
%-------------------------------------------------------------------------------
% the Euler approximation

% stepsize
    h = totalTime/numEulerSteps;
% initialize Euler t and v vectors
    tEuler(1)=0;
    vEuler(1)=0;

    deployIndex = ceil(timeOpenParachute*numEulerSteps/totalTime);
    deployIndex = ceil(numEulerSteps/2);
% calculate the approximation to the point when the parachute opens
    for j=2:deployIndex
       vEuler(j) = vEuler(j-1) + h*skydiverClosedODE(tEuler(j-1),vEuler(j-1));
       tEuler(j) = tEuler(j-1) + h;
    end
% and thereafter
    for j=(deployIndex+1):(numEulerSteps+1)
        vEuler(j) = vEuler(j-1) + h*skydiverOpenODE(tEuler(j-1),vEuler(j-1));
        tEuler(j) = tEuler(j-1) + h;
    end

% plot the Euler approximation
    plot( tEuler, vEuler, [EulerColor 'x'], 'MarkerSize',20, 'LineWidth',2 );
    annotation( figh, 'textbox', [0.55 0.13 0.26 0.072],...
                'String',{['Euler approximation, n=', ...
	                   num2str(numEulerSteps)]},...
                'FontSize',20, 'Color',EulerColor, 'FitBoxToText','on');

end

% all done
%-------------------------------------------------------------------------------
