function LunaOcampo_Pedro_Interplanetary(n,z)
% This function provides the following data for interplanetary missions:
% 1. Total ?v for the Hohmann transfer to and back from the planet
% 2. Total Travel Time
% 3. Synodic Time Period
% 4. Wait time
% 5. Departure characteristics
% a. Angular momentum, eccentricity of departure hyperbola from a ?z? km altitude circular
% parking orbit around earth
% b. ?v from parking orbit to departure hyperbola
% 6. Arrival characteristics
% a. Angular momentum, eccentricity of arrival hyperbola to get into an optimal circular
% orbit around the planet. If an optimal circular orbit does not exist, use a capture orbit
% with altitude ?z?.
% b. ?v to get into the circular orbit and radius, time period of the capture orbit
%
% The labeling of the planets (n) is as follow:
% 1 Mercury
% 2 Venus
% 3 Earth
% 4 Mars
% 5 Jupiter
% 6 Saturn
% 7 Uranus
% 8 Neptune
%
% z is the altitude of the parking orbit around Earth
clc
n=4;z=650;
if n == 3
display(' ')
error('Earth is not a valid destination!')
end
if n > 3
n = n-1;
end
rade = 6378; %km
re = 149.6e6 ; %km
r = [57.91e6 108.2e6 227.9e6 778.6e6 1.433e9 2.872e9 4.495e9]; %km
rap = [2440 6052 3396 71490 60270 25560 24760]; %km
SR = 23.9345; %hrs
T1 = 365.256 * SR; %hrs
T2 = [87.97*SR 224.7*SR 1.881*T1 11.86*T1 29.46*T1 84.01*T1 164.8*T1]; %hrs
mu = 132712000000; %km^3/s^2
mup = [22030 324900 42828 126686000 37931000 5794000 6835100]; %km^3/s^2
mue = 398600; %km^3/s^2
dVa = abs(sqrt(2 * (mu / re - mu / (re + r(n)))) - sqrt(mu / re)); %kms/s
dVb = abs(sqrt(2 * (mu / r(n) - mu / (re + r(n)))) - sqrt(mu / r(n))); %kms/s
tdV = 2*(dVa + dVb); %kms/s
Tsyn = T1 * T2(n) / abs(T1 - T2(n)); %hrs
Tsynd = Tsyn / 24 ; %days
n1 = 2 * pi / T1; %rev/hr
n2 = 2 * pi / T2(n); %rev/hr
a = (re + r(n)) / 2;
t12 = (pi * a^1.5) / sqrt(mu);%sec
t12 = t12/3600;
t12d = t12/24; %days
phi_f = pi - n1 * t12;
N = 1;
if n1 > n2
twait = (-2*phi_f - 2 * pi * N) / (n2 - n1);
while twait <= 0
twait = (-2*phi_f - 2 * pi * N) / (n2 - n1);
N = N + 1;
end
else
twait = (-2*phi_f + 2 * pi * N) / (n2 - n1);
while twait <= 0
twait = -2*phi_f + 2 * pi * N / (n2 - n1);
N = N + 1;
end
end
twait = twait/24;
disp(' ')
disp('Interplanetary Trajectory from Earth to Planet')
disp('----------------------------------------------')
disp(' ')
fprintf('Total Delta V for the Hohmann transfer = %4.2f km/s ',tdV)
disp(' ')
fprintf('Total travel time = %4.2f days',t12d)
disp(' ')
fprintf('Synodic Period = %4.2f days',Tsynd)
disp(' ')
fprintf('Wait time = %4.2f days',twait)
disp(' ')
%%% Departure %%%
rc = rade + z;
Vinfd = dVa;
ed = 1 + (rc * Vinfd ^ 2) / mue;
hd = mue * sqrt(ed ^ 2 - 1) / Vinfd;
Vp = hd / rc;
Vc = sqrt(mue/rc);
dVd = abs(Vc - Vp);
disp(' ')
disp('Departure Characteristics')
disp('-------------------------')
disp(' ')
fprintf('Angular momentum of departure hyperbola = %4.2f km2/s2 ',hd)
disp(' ')
fprintf('Eccentricity of departure hyperbola = %4.2f ',ed)
disp(' ')
fprintf('Delta V = %4.2f km/s',dVd)
disp(' ')
%%% Arrival %%%
Vinfa = dVb;
ra = 2 * mup(n) / Vinfa^2;
if ra < rap(n)
ra = rap(n) + z;
end
ea = 1 + ra * Vinfa^2 / mup(n);
ha = mup(n) * sqrt(ea^2 - 1) / Vinfa;
Tpc = 2 * pi * ra ^ (1.5) / sqrt(mup(n));
dVaR = abs(sqrt(Vinfa^2 + 2 * mup(n) / ra) - sqrt(mup(n) / ra));
disp(' ')
disp('Arrival Characteristics')
disp('-------------------------')
disp(' ')
fprintf('Angular momentum of arrival hyperbola = %4.2f km2/s2 ',ha)
disp(' ')
fprintf('Eccentricity of arrival hyperbola = %4.2f ',ea)
disp(' ')
fprintf('Radius of optimal capture orbit = %4.2f km ',ra)
disp(' ')
fprintf('Time period of capture orbit = %4.2f km ',Tpc)
disp(' ')
fprintf('Delta V = %4.2f km/s',dVaR)
disp(' ')
end
