%% Skydive2Demo.m
%% Overview
% This script illustrates the solution of a first-order IVP
%
% $$ v'(t) = g - \frac{c}{m} v^2(t),  \quad v(0) = 0 $$
%
% for finding the terminal velocity of a skydiver. The analytical solution
% is given by
%
% $$ v(t) = \sqrt{\frac{g m}{c}} \tanh \left( t \sqrt{\frac{g c}{m}}\right). $$
%
% Here $t$ stands for time, $v$ for velocity, and
%
% * $g$: gravitational constant 
% * $m$: mass
% * $c$: seond-order drag coefficient
%
% The right-hand side function is coded in |Skydive2.m|

%% Initialization 
clear all
clc
close all
g = 9.81; m = 68.1; c = 0.225;              % specific constants
v = @(t) sqrt(g*m/c)*tanh(t*sqrt(g*c/m));   % analytical solution
f = @Skydive2;                              % function handle for right-hand side function
t0 = 0; y0 = 0;                             % initial time and initial velocity
tend = 15;                                  % end time
N = 100;                                    % number of time steps
h = (tend-t0)/N;                            % stepsize

%% Call Euler
% We will discuss this function later
[t,y] = Euler(t0,y0,f,h,N);

%% Plot
hold on 
set(gca,'Fontsize',14)
xlabel('t','FontSize',14);
ylabel('Velocity','FontSize',14);
plot(t,y,'r','LineWidth',2);                % Euler solution 

%%
pause
tt=linspace(t0,tend,201);
plot(tt,v(tt),'g','LineWidth',2);           % analytical solution
legend('Euler solution','Analytical solution','Location','NorthWest');
hold off    
disp(sprintf('Terminal velocity: %f m/s.', sqrt(g*m/c)))









