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

%-------------------------------------------------------------------------------
% 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);

% we insert a pause after the first graph
    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' );
    
    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);

    disp('Press Enter to continue');
    pause;
    
%-------------------------------------------------------------------------------
% the Euler approximation

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

    deployIndex = ceil(timeOpenParachute*numEulerSteps/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):(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');

    disp('Press Enter to continue');
    pause;

% 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);
    localErrorArrowEulerv = [ lerrUseTime, vEuler(ceil(localErrorTime/h)+1) ];

% plot the corresponding point on the exact solution
    if lerrUseTime < timeOpenParachute
        lcol = closedColor;
	vexa = ClosedExactSol(lerrUseTime);
    else
        lcol = openColor;
	vexa = OpenExactSol(lerrUseTime);
    end
    localErrorArrowExactv = [ lerrUseTime, vexa ];

    plot( lerrUseTime, vexa, ['.' lcol], 'MarkerSize', 30, 'LineWidth', 2 );

% do we put the arrow between the points, or around them?
    locErr = abs(localErrorArrowExactv(2) - localErrorArrowEulerv(2));
    if locErr >= 15
        localErrorArrowMid = (localErrorArrowExactv + localErrorArrowEulerv)./2;

        a1 = arrow( localErrorArrowMid, localErrorArrowExactv,...
	            'FaceColor',annotationColor, 'EdgeColor',annotationColor );
        a2 = arrow( localErrorArrowMid, localErrorArrowEulerv,...
	            'FaceColor',annotationColor, 'EdgeColor',annotationColor );
    else
        if localErrorArrowExactv(2) - localErrorArrowEulerv(2) > 0
	    a1 = arrow( localErrorArrowExactv+[0,10], localErrorArrowExactv,...
                        'FaceColor', annotationColor, ...
	 	 	'EdgeColor', annotationColor );
	    a2 = arrow( localErrorArrowEulerv-[0,10], localErrorArrowEulerv,...
                        'FaceColor', annotationColor, ...
			'EdgeColor', annotationColor );
        else
	    a1 = arrow( localErrorArrowExactv-[0,10], localErrorArrowExactv,...
                        'FaceColor', annotationColor, ...
			'EdgeColor', annotationColor );
	    a2 = arrow( localErrorArrowEulerv+[0,10], localErrorArrowEulerv,...
    			'FaceColor', annotationColor, ...
			'EdgeColor', annotationColor );
        end
    end

% and add an annotation box explaining this
    locErr = round(100*locErr)/100;
    errAnn = annotation( figh, 'textbox', [0.6, 0.6, 0.15, 0.05],...
		         'String', {['Local Error at t=',...
                                     num2str(localErrorTime),...
                                     ': ', num2str(locErr)]},...
	                 'FontSize', 20, 'Color', annotationColor,...
             	  	 'FitBoxToText', 'on' );

    disp('Press Enter to continue');
    pause;

% clear the annotation box and error arrows
    delete(errAnn);
    delete(a1);
    delete(a2);

% finally, add the cumulative error
    cumErrorArrowExactv = [ totalTime, OpenExactSol(totalTime) ];
    cumErrorArrowEulerv = [ totalTime, vEuler(numEulerSteps+1) ];
    cumErr = cumErrorArrowExactv(2) - cumErrorArrowEulerv(2);
    if cumErr > 0
	arrow( cumErrorArrowExactv + [0,10], cumErrorArrowExactv,...
               'FaceColor',annotationColor, 'EdgeColor',annotationColor );
	arrow( cumErrorArrowEulerv - [0,10], cumErrorArrowEulerv,...
               'FaceColor',annotationColor, 'EdgeColor',annotationColor );
    else
	arrow( cumErrorArrowExactv - [0,10], cumErrorArrowExactv,...
               'FaceColor',annotationColor, 'EdgeColor',annotationColor );
	arrow( cumErrorArrowEulerv + [0,10], cumErrorArrowEulerv,...
               'FaceColor',annotationColor, 'EdgeColor',annotationColor );
    end

    cumErr = round(100*abs(cumErr))/100;
    annotation( figh, 'textbox', [0.6, 0.6, 0.15, 0.05],...
                'String', {['Cumulative Error: ', ...
                           num2str(cumErr)]},...
	        'FontSize', 20, 'Color', annotationColor,...
	  	'FitBoxToText', 'on' );
end

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