% ENM 7.11
% Denna kod skrevs på övningen 2017-10-05 i SF1514/18/19

clear

% Gör framåt Euler för att hitta en approximativ lösning till ODE:n
dt = 0.01; % steglängd
t0 = 0;
t_end = 200;
num_steps = (t_end - t0)/dt;
t = t0:dt:t_end;
y = zeros(2, length(t)); % en rad per komponent i y, alltså 2 rader nu!
y(1,1) = 2000; % starthöjd
y(2,1) = 0; % starthastighet

% Själva metoden:
for n=1:num_steps
    % Nu uppdaterar vi en kolonn i matrisen åt gången, för själva lösningen
    % är en kolonnvektor
    y(:,n+1) = y(:,n) + dt*f(t(n), y(:,n));
end
h = y(1,:); % plocka ut höjden vid varje tidpunkt
v = y(2,:); % plocka ut hastigheten vid varje tidpunkt
plot(t,h)
grid on

% När når personen marken?
i = find(h < 0, 1) % bara första indexet!
p = polyfit([t(i-1) t(i)], [h(i-1) h(i)], 1) % rät linje mellan punkterna!
t_zero = roots(p) % hitta nollstället


% Funktion för högerledet i y' = f(t,y)
function dydt = f(t,y)
    m = 75; g = 9.81; vf = 50; vs = 5;
    if t < 30 % fritt fall
        k = m*g/vf^2;
    else % fallskärm
        k = m*g/vs^2;
    end
    dydt = [y(2); (k/m)*y(2)^2 - g];
end
