%% Problem 3.8

% Earth's gravitational parameter (km^3/s^2)
mu_E = 398600;

% Earth's radius (km)
Re = 6378;

% Orbit altitudes (km)
h_perigee = 200;    % Perigee altitude
h_apogee = 600;     % Apogee altitude
h_limit = 400;      % Altitude limit

% Calculate orbital radii from Earth's center (km)
rp = Re + h_perigee;    % Perigee radius
ra = Re + h_apogee;     % Apogee radius
r_limit = Re + h_limit; % Radius at 400 km altitude

% Calculate the semi-major axis of the elliptical orbit
a = (rp + ra) / 2;

% Calculate orbital eccentricity
e = (ra - rp) / (ra + rp);

% Orbit equation:
% r = a(1 - e*cos(E))
% Solve for eccentric anomaly E at the 400 km altitude
cosE = (a - r_limit) / (a * e);
E = acos(cosE);

% Calculate the mean anomaly using Kepler's equation
M = E - e * sin(E);

% Calculate the mean motion (rad/s)
n = sqrt(mu_E / a^3);

% Calculate the time from 400 km altitude to apogee
t_half = (pi - M) / n;

% The spacecraft crosses 400 km twice,
% so multiply by 2 to get total time above 400 km
t_total = 2 * t_half;

% Display the total time in seconds
disp(['Time above 400 km = ', num2str(t_total), ' seconds'])

% Display the total time in minutes
disp(['Which is = ', num2str(t_total / 60), ' minutes'])