% 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, and then the local error at a specified point and cumulative
% error illustrated.
%
% 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 SkydiverNumericalMethods()

%-------------------------------------------------------------------------------
% 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;
    
% 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 insert a pause after each method
    pause on;
    
% 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' );
    
    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' );
    
    plot(tClosed, vClosed, [closedColor '-'], 'LineWidth',2);
    plot(tOpen, vOpen, [openColor '-'], 'LineWidth',2);

    disp('Press Enter to continue');
    pause;
    
    delete(eqntxt1);
    delete(eqntxt2);
    
%-------------------------------------------------------------------------------
% 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');

    disp('Press Enter to continue');
    pause;
    
%-------------------------------------------------------------------------------
% 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');

    disp('Press Enter to continue');
    pause;
    
%-------------------------------------------------------------------------------
% 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');

        disp('Press Enter to continue');
        pause;
    end
    
%-------------------------------------------------------------------------------
% add local and cumulative error values

% put in arrows showing the local error; we use the time where we have a
%    data point for the Euler solution that is approximately the correct
%    value for the time we want
    lerrUseTime = tEuler(ceil(localErrorTime/h)+1);

    if lerrUseTime < timeOpenParachute
	vexa = ClosedExactSol(lerrUseTime);
    else
	vexa = OpenExactSol(lerrUseTime);
    end

    localErrorExactvPt = [ lerrUseTime, vexa ];
    localErrorEulervPt = [ lerrUseTime, vEuler(ceil(localErrorTime/h)+1) ];
    localErrorImpEulervPt = [ lerrUseTime, vImpEuler(ceil(localErrorTime/h)+1)];
    localErrorRKvPt = [ lerrUseTime, vRK(ceil(localErrorTime/h)+1) ];

    topvVal = max( [ localErrorExactvPt(2), localErrorEulervPt(2), ...
                     localErrorImpEulervPt(2), localErrorRKvPt(2) ] );
    
    % loctarr = arrow( [ lerrUseTime, topvVal+20 ], [ lerrUseTime, topvVal ], ...
    loctarr = arrow( [ 10, -65 ], [ lerrUseTime, topvVal ], ...
                     'FaceColor', annotationColor, ...
                     'EdgeColor', annotationColor );
    
    eulLocalErr = round(abs(100*(localErrorEulervPt(2) - vexa)))/100;
    impeulLocalErr = round(abs(100*(localErrorImpEulervPt(2) - vexa)))/100;
    rkLocalErr = round(abs(100*(localErrorRKvPt(2) - vexa)))/100;
    
    % loctann = annotation( figh, 'textbox', [.25, .65, .15, .05], ...
    %                       String', { ['Local Error at t=', ...
    %                       num2str(localErrorTime), ': ', err1, '; ', ...
    %                       err2, '; ', err3]}, 'FontSize', 20 );
    text( 5, -20, {['Local Error at t=', num2str(localErrorTime)]}, ...
          'FontSize', 20, 'Color', annotationColor );
    text( 6, -32, {['Euler: ', num2str(eulLocalErr)]}, ...
          'FontSize', 20, 'Color', EulerColor );
    text( 6, -44, {['ImpEuler: ', num2str(impeulLocalErr)]}, ...
          'FontSize', 20, 'Color', ImpEulerColor );
    text( 6, -56, {['RungeKutta: ', num2str(rkLocalErr)]}, ...
          'FontSize', 20, 'Color', RKColor );
    % loctann = annotation( figh, 'textbox', 
    
    disp('Press Enter to continue');
    pause;

    exactVf = OpenExactSol(totalTime);
    cumtarr = arrow( [ 25, -50 ], [ totalTime, exactVf ], ...
                     'FaceColor', annotationColor, ...
                     'EdgeColor', annotationColor );

    format long;
    eulCumErr = round(abs(100000*(vEuler(numSteps+1) - exactVf)))/100000;
    impeulCumErr = round(abs(100000*(vImpEuler(numSteps+1) - exactVf)))/100000;
    rkCumErr = round(abs(100000*(vRK(numSteps+1) - exactVf)))/100000;
    eulCumErr = abs((vEuler(numSteps+1) - exactVf));
    impeulCumErr = abs((vImpEuler(numSteps+1) - exactVf));
    rkCumErr = abs((vRK(numSteps+1) - exactVf));
    
    text( 20, -60, {['Cumulative Error (at t=', num2str(totalTime), ')']}, ...
          'FontSize', 20, 'Color', annotationColor );
    text( 21, -72, {['Euler: ', num2str(eulCumErr)]}, ...
          'FontSize', 20, 'Color', EulerColor );
    text( 21, -84, {['ImpEuler: ', num2str(impeulCumErr)]}, ...
          'FontSize', 20, 'Color', ImpEulerColor );
    text( 21, -96, {['RungeKutta: ', num2str(rkCumErr)]}, ...
          'FontSize', 20, 'Color', RKColor );

end

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