Appendix 3.1.2Matlab Files Used for Trajectory Calculations

Appendix 3.1.2Matlab files used for trajectory calculations

% Habitat Trajectory determination and display

clear all; close all;

AU = 1.495978e8;

% All orbital information is taken from MIDAS files

% Venus orbit information

av = .72333160*AU;

ev = .67663667E-02;

iv = 3.3945649*pi/180;

lanv = 76.641599*pi/180;

aphv = 54.930823*pi/180;

thsv1 = 343.50998*pi/180;

thsv2 = (175.07957+1080)*pi/180;

thsv = linspace(thsv1,thsv2, 793);

% Calculate Venus orbit wrt solar ecliptic

pv = av*(1-ev^2);

thv = thsv+aphv;

rv = pv./(1+ev.*cos(thsv));

xvf = rv.*cos(thv);

yvf = rv.*sin(thv);

yvi = yvf*cos(iv);

rvf = sqrt(xvf.^2+yvi.^2);

xv = rvf.*cos(thv+lanv);

yv = rvf.*sin(thv+lanv);

zv = sqrt(rv.^2-xv.^2-yv.^2).*sign(pi-mod(thv,2*pi));

% Earth orbit

earthorb

% Mars Orbit

marsorb

% Habitat Trajectory - Earth to Mars

a1 = 1.2991554*AU;

e1 = .24858019;

i1 = .47247910*pi/180;

lan1 = 293.34263*pi/180;

aph1 = 164.43720*pi/180;

ths11 = 15.757718*pi/180;

ths12 = 143.38146*pi/180;

ths1 = linspace(ths11, ths12, 172);

p1 = a1*(1-e1^2);

th1 = ths1+aph1;

r1 = p1./(1+e1.*cos(ths1));

x1f = r1.*cos(th1);

y1f = r1.*sin(th1);

y1i = y1f*cos(i1);

r1f = sqrt(x1f.^2+y1i.^2);

x1 = r1f.*cos(th1+lan1);

y1 = r1f.*sin(th1+lan1);

z1 = sqrt(r1.^2-x1.^2-y1.^2).*sign(pi-mod(th1,2*pi));

% Habitat Trajectory - Mars to Venus

a2 = 1.1333948*AU;

e2 = .36473009;

i2 = .96953978*pi/180;

lan2 = 263.80166*pi/180;

aph2 = 170.84250*pi/180;

ths21 = 166.51516*pi/180;

ths22 = (3.5842208+360)*pi/180;

ths2 = linspace(ths21, ths22, 255);

p2 = a2*(1-e2^2);

th2 = ths2+aph2;

r2 = p2./(1+e2.*cos(ths2));

x2f = r2.*cos(th2);

y2f = r2.*sin(th2);

y2i = y2f*cos(i2);

r2f = sqrt(x2f.^2+y2i.^2);

x2 = r2f.*cos(th2+lan2);

y2 = r2f.*sin(th2+lan2);

z2 = sqrt(r2.^2-x2.^2-y2.^2).*sign(pi-mod(th2,2*pi));

% Habitat Trajectory - Venus to Earth

a3 = .84188691*AU;

e3 = .20834487;

i3 = .95981576E-01*pi/180;

lan3 = 359.42170*pi/180;

aph3 = 23.281097*pi/180;

ths31 = 55.526393*pi/180;

ths32 = (156.63321+360)*pi/180;

ths3 = linspace(ths31, ths32, 366);

p3 = a3*(1-e3^2);

th3 = ths3+aph3;

r3 = p3./(1+e3.*cos(ths3));

x3f = r3.*cos(th3);

y3f = r3.*sin(th3);

y3i = y3f*cos(i3);

r3f = sqrt(x3f.^2+y3i.^2);

x3 = r3f.*cos(th3+lan3);

y3 = r3f.*sin(th3+lan3);

z3 = sqrt(r3.^2-x3.^2-y3.^2).*sign(pi-mod(th3,2*pi));

% Display trajectory in b&w

x = [x1 x2 x3];

y = [y1 y2 y3];

z = [z1 z2 z3];

px = []; py = []; pz = [];

for i = [1:length(x)]

if (mod(i,5) == 1)

px = [px x(i)];

py = [py y(i)];

pz = [pz z(i)];

end

end

plot3(px/AU,py/AU,pz/AU,'k:');

hold

plot3(x1(1)/AU,y1(1)/AU,z1(1)/AU,'ko');

plot3(x2(1)/AU,y2(1)/AU,z2(1)/AU,'ko');

plot3(x3(1)/AU,y3(1)/AU,z3(1)/AU,'ko');

plot3(x3(end)/AU,y3(end)/AU,z3(end)/AU,'ko');

plot3(xv/AU,yv/AU,zv/AU,'k');

plot3(xm/AU,ym/AU,zm/AU,'k');

plot3(xe/AU,ye/AU,ze/AU,'k');

% ERA Trajectory Earth to Mars

clear all; close all;

AU = 1.495978e8;

% All orbital information taken from MIDAS files

% Earth orbit

ae = 1.0000002*AU;

ee = .016704123;

ie = .0015476714*pi/180;

lane = 174.84774*pi/180;

aphe = 288.13136*pi/180;

thse1 = 303.54483*pi/180;

thse2 = (245.78624 + 360)*pi/180;

thse = linspace(thse1,thse2, 307);

pe = ae*(1-ee^2);

the = thse+aphe;

re = pe./(1+ee.*cos(thse));

xef = re.*cos(the);

yef = re.*sin(the);

yei = yef*cos(ie);

ref = sqrt(xef.^2+yei.^2);

xe = ref.*cos(the+lane);

ye = ref.*sin(the+lane);

ze = sqrt(re.^2-xe.^2-ye.^2).*sign(pi-mod(the,2*pi));

% Mars Orbit

am = 1.5236915*AU;

em = .93415785E-01;

im = 1.8487363*pi/180;

lanm = 49.524141*pi/180;

aphm = 286.58789*pi/180;

thsm1 = 136.17361*pi/180;

thsm2 = 278.70649*pi/180;

thsm = linspace(thsm1,thsm2, 307);

pm = am*(1-em^2);

thm = thsm+aphm;

rm = pm./(1+em.*cos(thsm));

xmf = rm.*cos(thm);

ymf = rm.*sin(thm);

ymi = ymf*cos(im);

rmf = sqrt(xmf.^2+ymi.^2);

xm = rmf.*cos(thm+lanm);

ym = rmf.*sin(thm+lanm);

zm = sqrt(rm.^2-xm.^2-ym.^2).*sign(pi-mod(thm,2*pi));

% Trajectory information

a = 1.2524394*AU;

e = .21018949;

i = 1.6689473*pi/180;

lan = 46.568310*pi/180;

aph = 7.2399656*pi/180;

ths1 = 352.71606*pi/180;

ths2 = (201.00870 + 360)*pi/180;

ths = linspace(ths1, ths2, 307);

p = a*(1-e^2);

th = ths+aph;

r = p./(1+e.*cos(ths));

xf = r.*cos(th);

yf = r.*sin(th);

yi = yf*cos(i);

rf = sqrt(xf.^2+yi.^2);

x = rf.*cos(th+lan);

y = rf.*sin(th+lan);

z = sqrt(r.^2-x.^2-y.^2).*sign(pi-mod(th,2*pi));

% Display Trajectory in b&w

px = []; py = []; pz = [];

for i = [1:length(x)]

if (mod(i,5) == 1)

px = [px x(i)];

py = [py y(i)];

pz = [pz z(i)];

end

end

plot3(xe/AU,ye/AU,ze/AU,'k');

hold

plot3(xm/AU,ym/AU,zm/AU,'k');

plot3(px/AU,py/AU,pz/AU,'k:');

plot3(x(1)/AU,y(1)/AU,z(1)/AU,'ko');

plot3(x(end)/AU,y(end)/AU,z(end)/AU,'ko');

plot3(0,0,0,'k.',0,0,0,'ko');

% ERV Trajectory Mars to Earth

clear all; close all;

AU = 1.495978e8;

% All orbital information taken form MIDAS files

% Earth orbit

ae = 1.0000002*AU;

ee = .016702386;

ie = .0020865941*pi/180;

lane = 174.83777*pi/180;

aphe = 288.15468*pi/180;

thse1 = 351.04569*pi/180;

thse2 = (217.38320 + 360)*pi/180;

thse = linspace(thse1,thse2, 230);

pe = ae*(1-ee^2);

the = thse+aphe;

re = pe./(1+ee.*cos(thse));

xef = re.*cos(the);

yef = re.*sin(the);

yei = yef*cos(ie);

ref = sqrt(xef.^2+yei.^2);

xe = ref.*cos(the+lane);

ye = ref.*sin(the+lane);

ze = sqrt(re.^2-xe.^2-ye.^2).*sign(pi-mod(the,2*pi));

% Mars Orbit

am = 1.5236915*AU;

em = .93419578E-01;

im = 1.8483974*pi/180;

lanm = 49.511935*pi/180;

aphm = 286.61844*pi/180;

thsm1 = 195.49286*pi/180;

thsm2 = 311.74200*pi/180;

thsm = linspace(thsm1,thsm2, 230);

pm = am*(1-em^2);

thm = thsm+aphm;

rm = pm./(1+em.*cos(thsm));

xmf = rm.*cos(thm);

ymf = rm.*sin(thm);

ymi = ymf*cos(im);

rmf = sqrt(xmf.^2+ymi.^2);

xm = rmf.*cos(thm+lanm);

ym = rmf.*sin(thm+lanm);

zm = sqrt(rm.^2-xm.^2-ym.^2).*sign(pi-mod(thm,2*pi));

% Trajectory information

a = 1.3425413*AU;

e = .25115112;

i = 3.0169324*pi/180;

lan = 140.39486*pi/180;

aph = 195.90620*pi/180;

ths1 = 195.36978*pi/180;

ths2 = 344.07345*pi/180;

ths = linspace(ths1, ths2, 230);

p = a*(1-e^2);

th = ths+aph;

r = p./(1+e.*cos(ths));

xf = r.*cos(th);

yf = r.*sin(th);

yi = yf*cos(i);

rf = sqrt(xf.^2+yi.^2);

x = rf.*cos(th+lan);

y = rf.*sin(th+lan);

z = sqrt(r.^2-x.^2-y.^2).*sign(pi-mod(th,2*pi));

% Display trajectory in b&w

px = []; py = []; pz = [];

for i = [1:length(x)]

if (mod(i,5) == 1)

px = [px x(i)];

py = [py y(i)];

pz = [pz z(i)];

end

end

plot3(xe/AU,ye/AU,ze/AU,'k');

hold

plot3(xm/AU,ym/AU,zm/AU,'k');

plot3(px/AU,py/AU,pz/AU,'k:');

plot3(x(1)/AU,y(1)/AU,z(1)/AU,'ko');

plot3(x(end)/AU,y(end)/AU,z(end)/AU,'ko');

plot3(0,0,0,'k.',0,0,0,'ko');

% earthorb.m

% determine Earth orbit for entire mission span

AU = 1.495978e8;

% earth orbital constants

ae = 1.0000002*AU;

ee = .016704483;

ie = .0014359985*pi/180;

lane = 174.84981*pi/180;

aphe = 288.12652*pi/180;

thse1 = 357.07811*pi/180;

thse2 = (357.54476 + 6*360)*pi/180;

thse = linspace(thse1,thse2, 2192);

pe = ae*(1-ee^2);

the = thse+aphe;

re = pe./(1+ee.*cos(thse));

xef = re.*cos(the);

yef = re.*sin(the);

yei = yef*cos(ie);

ref = sqrt(xef.^2+yei.^2);

xe = ref.*cos(the+lane);

ye = ref.*sin(the+lane);

ze = sqrt(re.^2-xe.^2-ye.^2).*sign(pi-mod(the,2*pi));

% marsorb.m

% determine Mars orbit for entire mission span

AU = 1.495978e8;

% Mars orbital constants

am = 1.5236915*AU;

em = .093414999;

im = 1.8488065*pi/180;

lanm = 49.526670*pi/180;

aphm = 286.58157*pi/180;

thsm1 = 317.81442*pi/180;

thsm2 = (39.820555 + 4*360)*pi/180;

thsm = linspace(thsm1,thsm2, 2192);

pm = am*(1-em^2);

thm = thsm+aphm;

rm = pm./(1+em.*cos(thsm));

xmf = rm.*cos(thm);

ymf = rm.*sin(thm);

ymi = ymf*cos(im);

rmf = sqrt(xmf.^2+ymi.^2);

xm = rmf.*cos(thm+lanm);

ym = rmf.*sin(thm+lanm);

zm = sqrt(rm.^2-xm.^2-ym.^2).*sign(pi-mod(thm,2*pi));

% Calculate impact of cutting tether on approach trajectory at Mars

AU = 1.495978e8;

mum = 42809.6920936;

rm = 3393.15984;

mue = 398600;

re = 6378.14;

vinf = 6.767978;

rcut = rm + 600000;

rp = rm + 3849.1464;

vp = sqrt(vinf^2 + 2*mum/rp);

h = rp*vp;

a1 = -mum/(vp^2-2*mum/rp);

e1 = 1 - rp/a1;

p1 = h^2/mum;

vcut = sqrt(vinf^2 + 2*mum/rcut);

gcut = acos(h/(rcut*vcut));

thcut = acos((p1/rcut - 1)/e1);

vr1 = vcut*(-sin(gcut));

vv1 = vcut*cos(gcut);

dv = .04325;

vr2 = vr1 + dv*(-cos(gcut));

vv2 = vv1 + dv*(-sin(gcut));

vcut2 = sqrt(vv2^2 + vr2^2);

gcut2 = acos(vv2/vcut2);

h2 = rcut*vcut2*cos(gcut2);

p2 = h2^2/mum;

a2 = -mum/(vcut2^2-2*mum/rcut);

e2 = sqrt(1- p2/a2);

rp2 = a2*(1-e2);

alt2 = rp2 - rm;

vinf2 = sqrt(vcut2^2-2*mum/rcut);

ren = rm+100;

ven = sqrt(vinf^2+2*mum/ren);

gen = acos(h2/(ren*ven));

gend = gen*180/pi;