% O18
close all 
clear
constants

z_0 = [1 2 3].*1e5; % Initial altitude.

rho_0 = 1e-7; % [kg m^-3]
H   = 36e3; % m
m   = 1000; % [kg] 
A = 1; % [m^2]
C_D = 2.2; % Drag coeff.



td = (H.*(exp(z_0./H)-1).*m)./(rho_0.*sqrt(G .* M_E .* R_E).*C_D.*A)

td_h = td./60/60;
td_d = td_h./24

% Plot descent over time 
for caseFlag = 1:3
switch(caseFlag)
    case 1
    t = 0:60:td(1);
    case 2
    t = 0:60:td(2);
    case 3
    t = 0:60:td(3);
end

inner = -rho_0*sqrt(G*M_E)*C_D*A/m*t./(H/sqrt(R_E)) + exp(z_0(caseFlag)/H);
z = H*log(inner);

plot(t, z)
ylabel("Altitude above Earth's surface z");
xlabel("Time t");
title("Orbital decay due to atmospheric drag");
hold on
grid on
end
