% This code integrates the equation in problem 3 of HW 2
%   v' = -a*v + g,
%  where the friction constant a abruptly changes at t=topen
%  due to the opening of the parachute.
%  vmax=80 mph  for t<topen  and  vmax=4 mph for t>topen.
% The integration method is MATLAB's  ode45.

clear all

global a1 a2 g topen

% Convert the variables from one usints to the others:
vmax1_mph = 80;                  % max speed in miles per hour before parachute opens
vmax2_mph = 4;                   % max speed in miles per hour after parachute opens
vmax1 = vmax1_mph*(1609/3600); 
vmax2 = vmax2_mph*(1609/3600);   % same speeds in meters per second
g = 9.8;                         % acceleration of gravity in m/s^2
a1 = g/vmax1;                    % friction constant before parachute opens
a2 = g/vmax2;                    % friction constant after parachute opens

% Other parameters of the problem:
tmin = 0;
tmax = 4;                        % beginning and ending times [s]
topen = 2;                       % time parachute opens [s]

tspan = [tmin tmax];             % this is used by the solvers
v0 = 0;                          % this is the initial speed
epstolA = 10^(-3)*(1609/3600);   % allowed local error converted from mph to  m/s
epstolR = 10^(-4);               % allowed relative error 


options = odeset('RelTol', epstolR, 'AbsTol', epstolA);
% options = odeset('AbsTol', epstolA);
[t_ode45,v_ode45] = ode45(@fun4_hw2_p3, tspan, v0, options);

% Convert velocity in meters per second in that in miles per hour
% to compare it with the exact solution, which below is computed, at the end, in mph:
v_ode45_mph = v_ode45*(3600/1609);

% Exact Solution
for i=1:length(t_ode45)
    if t_ode45(i) < topen
        y(i) = vmax1*( 1 - exp(-a1*t_ode45(i)) );
    else
        y(i)=(vmax1*(1 - exp(-a1*topen)) - vmax2) * exp( -a2*(t_ode45(i)-topen)) + vmax2;
    end
end
y_mph = y*(3600/1609);

figure(20401);
subplot(211);
plot(t_ode45, y_mph, t_ode45, v_ode45_mph, 'r--')
legend('Exact Solution','ode45')
title('HW 2, Problem 4 - Comparing Exact Solution to ODE45');
xlabel('time, [s]', 'fontsize', 16)
ylabel('diver`s speed, [mph]', 'fontsize', 16)
subplot(212);
plot(t_ode45, y_mph' - v_ode45_mph)   
            % y_mph needs to be transposed to have the same dimensions as v_ode45_mph
title('HW 2, Problem 4 - error of the ODE45 solution');
xlabel('time, [s]', 'fontsize', 16)
ylabel('exact - numerical, [mph]', 'fontsize', 16)