% EX 7.11

% Fallskärmshopp mh'' = k(h')^2 - mg, h(0) = H, h'(0) = 0 
% Lös med ode45, standardfelmarginaler 

% Substitution
% u1 = h, u2 = h'
% u1' = u2             , u1(0) = H
% u2' = k*(u2)^2/m - g , u2(0) = 0

% Konstanter
g = 9.81;
H = 2000;
m = 75;

% Luftmotstånd, beräknas från sluthastigheten (acc = 0)
% k = mg/v_slut^2
k_fritt = m*g/50^2;
k_skarm = m*g/5^2;

% Definiera funktion
f_fritt = @(t, u) [u(2); k_fritt*u(2).^2/m - g];
f_skarm = @(t, u) [u(2); k_skarm*u(2).^2/m - g];

figure(1)
clf
subplot(2,1,1)
hold on
axis([0 200 0 H])
xlabel('Tid från hopp', 'FontSize', 14)
ylabel('Höjd över marken', 'FontSize', 14)
subplot(2,1,2)
hold on
axis([0 200 -50 0])
xlabel('Tid från hopp', 'FontSize', 14)
ylabel('Vertikal hastighet', 'FontSize', 14)

% Fas 1: Fritt fall från höjden H
[t_fritt, u_fritt] = ode45(f_fritt, [0 30], [H 0]);

% Höjd och fart vid skärmens utvecklande
H_skarm = u_fritt(end, 1); 
V_skarm = u_fritt(end, 2); 

subplot(2,1,1)
plot(t_fritt, u_fritt(:, 1), 'LineWidth', 2)
subplot(2,1,2)
plot(t_fritt, u_fritt(:, 2), 'LineWidth', 2)


% Fas 1: Fritt fall från höjden H
[t_skarm, u_skarm] = ode45(f_skarm, [30 200], [H_skarm V_skarm]);

% Tidpunkt och fart vid landning, linjärinterpol
landat = find(u_skarm(:, 1) <= 0); landat = landat(1);
bakatvikt = (0 - u_skarm(landat, 1))/...
	(u_skarm(landat - 1, 1) - u_skarm(landat, 1));
framatvikt = (u_skarm(landat - 1, 1) - 0)/...
	(u_skarm(landat - 1, 1) - u_skarm(landat, 1));
t_landning = bakatvikt*t_skarm(landat - 1) + framatvikt*t_skarm(landat);
v_landning = bakatvikt*u_skarm(landat - 1, 2) + framatvikt*u_skarm(landat, 2);
% Hastigheten blir 0 vid landning
u_skarm(landat:end, 2) = 0;

subplot(2,1,1)
plot(t_skarm, u_skarm(:, 1), 'g', 'LineWidth', 2)
subplot(2,1,2)
plot(t_skarm, u_skarm(:, 2), 'g', 'LineWidth', 2)



% Kolla resultatets noggrannhet genom att göra om beräkningen med annan
% relativfeltolerans (1e-3 standard) 
tolerans = odeset('RelTol', 1e-6);
[t_fritt2, u_fritt2] = ode45(f_fritt, [0 30], [H 0], tolerans);
H_skarm2 = u_fritt2(end, 1); 
V_skarm2 = u_fritt2(end, 2); 
[t_skarm2, u_skarm2] = ode45(f_skarm, [30 200], [H_skarm2 V_skarm2]);
landat = find(u_skarm2(:, 1) <= 0); landat = landat(1);
bakatvikt = (0 - u_skarm2(landat, 1))/...
	(u_skarm2(landat - 1, 1) - u_skarm2(landat, 1));
framatvikt = (u_skarm2(landat - 1, 1) - 0)/...
	(u_skarm2(landat - 1, 1) - u_skarm2(landat, 1));
t_landning2 = bakatvikt*t_skarm2(landat - 1) + framatvikt*t_skarm2(landat);
v_landning2 = bakatvikt*u_skarm2(landat - 1, 2) + framatvikt*u_skarm2(landat, 2);

H_skarmfel = abs(H_skarm - H_skarm2)
t_landningfel = abs(t_landning - t_landning2)
v_landningfel = abs(v_landning - v_landning2)

fprintf(['Fallskärmen vecklas ut vid höjden %0.' ...
	num2str(round(-1 - log10(H_skarmfel))) 'f m.\n'], H_skarm)
fprintf(['Landning sker efter %0.' num2str(round(-1 - log10(t_landningfel))) ...
	'f s från hopp, med den vertikala hastigheten %0.' ...
	num2str(round(-1 - log10(v_landningfel))) 'f m/s.\n'], t_landning, v_landning)


% Körning:
%H_skarmfel =
%    0.0533
%t_landningfel =
%    0.0107
%v_landningfel =
%   1.0836e-04
%Fallskärmen vecklas ut vid höjden 677 m.
%Landning sker efter 164.4 s från hopp, med den vertikala hastigheten -5.002 m/s.
