% 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 different approximations
% shown.
%
% by Gavin LaRose, 9/10/12
% parts excerpted from previous demonstrations by Justin Dunmyre

% (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 SkydiverNumericalPlot()

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

% number of 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) % '
    numSteps = 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 may be shortly before the
%    parachute deploys)
    localErrorTime = 14;
    
% should we include RK in the calculation?  set to zero to exclude
    includeRK = 1;
    
% should we pause between grsphs?  set to zero to graph all at once
    pauseForGraphs = 1;

% colors
    closedColor='c';
    openColor='g';
    EulerColor='w';
    ImpEulerColor='b';
    RKColor='r';
    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);

% we may insert a pause after each method
    if pauseForGraphs == 1
        pause on;
    end
    
% 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' );

% only put annotations for the first graph in if we're pausing
    if pauseForGraphs == 1
        eqntxt1 = 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' );
        eqntxt2 = annotation( figh, 'textbox', [0.62 0.45 0.28 0.092], ...
                              'String', ...
                              {'v\prime=-32.2 - 2.01 v', ...
                               '(after parachute)'}, ...
                              'FontSize', 20, 'Color', openColor,...
                              'FitBoxToText','on' );
    end
    
    plot(tClosed, vClosed, [closedColor '-'], 'LineWidth',2);
    plot(tOpen, vOpen, [openColor '-'], 'LineWidth',2);

    if pauseForGraphs == 1
        disp('Press Enter to continue');
        pause;
        delete(eqntxt1);
        delete(eqntxt2);
    end
  
%-------------------------------------------------------------------------------
% the Euler approximation

% stepsize
    h = totalTime/numSteps;

% initialize Euler t and v vectors
    tEuler(1)=0;
    vEuler(1)=0;

    deployIndex = ceil(timeOpenParachute*numSteps/totalTime);
% 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):(numSteps+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',10, 'LineWidth',2 );
    annotation( figh, 'textbox', [0.55 0.3 0.26 0.072],...
                'String',{['Euler approximation, n=', ...
	                   num2str(numSteps)]},...
                'FontSize',20, 'Color',EulerColor, 'FitBoxToText','on');

    if pauseForGraphs == 1
        disp('Press Enter to continue');
        pause;
    end

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

% initialize improved Euler t and v vectors
    tImpEuler(1)=0;
    vImpEuler(1)=0;

% calculate the approximation to the point when the parachute opens
    for j=2:deployIndex
        k = skydiverClosedODE(tImpEuler(j-1),vImpEuler(j-1));
        vImpEuler(j) = vImpEuler(j-1) + ...
            0.5*h*(k + skydiverClosedODE(tImpEuler(j-1)+h,vImpEuler(j-1) + h*k));
        tImpEuler(j) = tImpEuler(j-1) + h;
    end
% and thereafter
    for j=(deployIndex+1):(numSteps+1)
        k = skydiverOpenODE(tImpEuler(j-1),vImpEuler(j-1));
        vImpEuler(j) = vImpEuler(j-1) + ...
            0.5*h*(k + skydiverOpenODE(tImpEuler(j-1)+h,vImpEuler(j-1) + h*k));
        tImpEuler(j) = tImpEuler(j-1) + h;
    end

% plot the approximation
    plot( tImpEuler, vImpEuler, [ImpEulerColor 'x'], 'MarkerSize',10, ...
          'LineWidth',2 );
    annotation( figh, 'textbox', [0.55 0.2 0.26 0.072],...
                'String',{['ImpEuler approximation, n=', ...
	                   num2str(numSteps)]},...
                'FontSize',20, 'Color',ImpEulerColor, 'FitBoxToText','on');

    if pauseForGraphs == 1
        disp('Press Enter to continue');
        pause;
    end
    
%-------------------------------------------------------------------------------
% the Runge Kutta approximation

    if includeRK == 1
    
% initialize RK t and v vectors
        tRK(1)=0;
        vRK(1)=0;

        for j=2:deployIndex
            k1 = skydiverClosedODE(tRK(j-1),vRK(j-1));
            k2 = skydiverClosedODE(tRK(j-1)+.5*h,vRK(j-1)+.5*h*k1);
            k3 = skydiverClosedODE(tRK(j-1)+.5*h,vRK(j-1)+.5*h*k2);
            k4 = skydiverClosedODE(tRK(j-1)+h,vRK(j-1)+h*k3);
            vRK(j) = vRK(j-1) + (h/6)*(k1 + 2*(k2 + k3) + k4);
            tRK(j) = tRK(j-1) + h;
        end

        for j=(deployIndex+1):(numSteps+1)
            k1 = skydiverOpenODE(tRK(j-1),vRK(j-1));
            k2 = skydiverOpenODE(tRK(j-1)+.5*h,vRK(j-1)+.5*h*k1);
            k3 = skydiverOpenODE(tRK(j-1)+.5*h,vRK(j-1)+.5*h*k2);
            k4 = skydiverOpenODE(tRK(j-1)+h,vRK(j-1)+h*k3);
            vRK(j) = vRK(j-1) + (h/6)*(k1 + 2*(k2 + k3) + k4);
            tRK(j) = tRK(j-1) + h;
        end

% plot the approximation
        plot( tRK, vRK, [RKColor 'x'], 'MarkerSize',10, ...
              'LineWidth',2 );
        annotation( figh, 'textbox', [0.55 0.1 0.26 0.072],...
                    'String',{['RungeKutta approximation, n=', ...
                            num2str(numSteps)]},...
                    'FontSize',20, 'Color',RKColor, 'FitBoxToText','on');

        if pauseForGraphs == 1
            disp('Press Enter to continue');
            pause;
        end
    end

end

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