2.2.8Appendix
aae450atrajs.m – code responsible for propagating the Crew Transfer Vehicle’s arrival into the Earth’s atmosphere
%Nicholas Saadah
%AAE 450
%Presentation #1 (1/16/01)
%CTV arrival to Earth
clear all
global fcflag fhflag faflag feflag fvflag
%Constants
umars = -398600; %km^3/s^2
g = 9.8; %m/s^2
rmars = 6378; %km
msc = 5000; %kg
areasc = 25.66; %m^2
fcflag = 0;
fhflag = 0;
faflag = 0;
feflag = 0;
fvflag = 0;
%Initial Conditions
xi = rmars+200; %km
yi = 0; %km
vi = 11.65; %km/s
gammaideg= input('gamma: '); %deg
maxgload = input('maxgload: ');
betadeg = input('beta: '); %deg
gammai = gammaideg*pi/180; %radians
beta = betadeg*pi/180; %radians
vxi = vi*sin(gammai); %km/s
vyi = vi*cos(gammai); %km/s
p0 = [xi vxi yi vyi];
ti = [0 1e5];
options = odeset('Reltol',1e-8);
[t,xx] = ode45('aae450eoms',ti,p0,options,umars,rmars,msc,areasc,g,beta,maxgload);
x = xx(:,1);
vx = xx(:,2);
y = xx(:,3);
vy = xx(:,4);
r = sqrt(x.^2 + y.^2);
alt = r-rmars;
clapp = min(alt);
if min(alt) < 0 | fvflag == 1
[a b] = min(abs(alt));
x = x(1:b);
vx = vx(1:b);
y = y(1:b);
vy = vy(1:b);
t = t(1:b);
r = r(1:b);
alt = alt(1:b);
contact = 1;
else
contact = 0;
end
if fvflag == 1
contact = 0;
end
if fhflag == 1
[c d] = size(xx);
for i = 1:1:c-1
if vx(i+1) - vx(i) == 0
htrj = i;
break
end
end
x = x(1:htrj);
vx = vx(1:htrj);
y = y(1:htrj);
vy = vy(1:htrj);
t = t(1:htrj);
r = r(1:htrj);
alt = alt(1:htrj);
end
v = sqrt(vx.^2 + vy.^2);
ro = .047*exp(-.1*((r-rmars)-49));
q = .5.*ro.*(v.*1000).^2;
q = q + 1e-12;
k = .125*((msc*maxgload*g)./(q*areasc)).^2;
%options = optimset('Display','off');
for i = 1:1:length(q)
if k(i) > 2
alpha(i) = pi/2;
elseif k(i) > 0 & k(i) < 2
alpha(i) = fzero('aae450fzeros',[0 pi/2],[],[],k(i));
end
end
alpha = alpha';
cl = 2*sin(2*alpha);
lift = cl*areasc.*sin(alpha).*q;
cd = 2*(1-cos(2*alpha));
drag = cd*areasc.*sin(alpha).*q;
taf = sqrt(lift.^2 + drag.^2);
adt = taf/msc;
alt = r-rmars;
gs = adt/9.8;
tmin = t/60;
thrs = t/3600;
alpha = alpha*180/pi;
rvec = [x y zeros(length(x),1)];
vvec = [vx vy zeros(length(vx),1)];
h = cross(rvec,vvec);
hmag = sqrt(h(:,1).^2 + h(:,2).^2 + h(:,3).^2);
scripte = v.^2/2 + umars./r;
ecc = sqrt(1 + 2*scripte.*hmag.^2/umars^2);
[e f] = max(alt);
alts = alt(1:f);
[e f] = min(alts);
alts = alts(f:length(alts));
[a100err c100] = min(abs(alts-100));
c100 = c100+f-1;
rvecexit = [x(c100) y(c100) 0];
vvecexit = [vx(c100) vy(c100) 0];
hvecexit = cross(rvecexit,vvecexit);
hexitmag = sqrt(hvecexit(1)^2 + hvecexit(2)^2 + hvecexit(3)^2);
gammaexit = acos(hexitmag/(r(c100)*v(c100)));
gammaexit = gammaexit*180/pi;
vexit = v(c100);
theta = linspace(0,2*pi,100);
xms = rmars.*cos(theta);
yms = rmars.*sin(theta);
figure
plot(x,y)
hold on
axis equal
fill(xms,yms,'g')
if contact == 1
plot(x(b),y(b),'r*')
end
aae450atrajr.m – code responsible for propagating the Earth Return Vehicle’s arrival into the Martian atmosphere
%Nicholas Saadah
%AAE 450
%Presentation #1 (1/16/01)
%ERV arrival at Mars
clear all
global fcflag fhflag faflag feflag
%Constants
umars = -42828; %km^3/s^2
g = 9.8; %m/s^2
rmars = 3397; %km
msc = 55200; %kg
areasc = 81.71; %m^2
fcflag = 0;
fhflag = 0;
faflag = 0;
feflag = 0;
%Initial Conditions
xi = rmars+100; %km
yi = 0; %km
vi = 5.64; %km/s
gammaideg= input('gamma: '); %deg
maxgload = input('maxgload: ');
betadeg = input('beta: '); %deg
gammai = gammaideg*pi/180; %radians
beta = betadeg*pi/180; %radians
vxi = vi*sin(gammai); %km/s
vyi = vi*cos(gammai); %km/s
p0 = [xi vxi yi vyi];
ti = [0 3e5];
options = odeset('Reltol',1e-8);
[t,xx] = ode45('aae450eomt',ti,p0,options,umars,rmars,msc,areasc,g,beta,maxgload);
x = xx(:,1);
vx = xx(:,2);
y = xx(:,3);
vy = xx(:,4);
r = sqrt(x.^2 + y.^2);
alt = r-rmars;
clapp = min(alt);
if min(alt) < 0
contact = 1;
[a b] = min(abs(alt));
x = x(1:b);
vx = vx(1:b);
y = y(1:b);
vy = vy(1:b);
t = t(1:b);
r = r(1:b);
alt = alt(1:b);
else
contact = 0;
end
if fhflag == 1
[c d] = size(xx);
for i = 1:1:c-1
if vx(i+1) - vx(i) == 0
htrj = i;
break
end
end
x = x(1:htrj);
vx = vx(1:htrj);
y = y(1:htrj);
vy = vy(1:htrj);
t = t(1:htrj);
r = r(1:htrj);
alt = alt(1:htrj);
end
v = sqrt(vx.^2 + vy.^2);
ro = .00047*exp(-.1*((r-rmars)-49));
q = .5.*ro.*(v.*1000).^2;
q = q + 1e-12;
k = .125*((msc*maxgload*g)./(q*areasc)).^2;
%options = optimset('Display','off');
for i = 1:1:length(q)
if k(i) > 2
alpha(i) = pi/2;
elseif k(i) > 0 & k(i) < 2
alpha(i) = fzero('aae450fzeros',[0 pi/2],[],[],k(i));
end
end
alpha = alpha';
cl = 2*sin(2*alpha);
lift = cl*areasc.*sin(alpha).*q;
cd = 2*(1-cos(2*alpha));
drag = cd*areasc.*sin(alpha).*q;
taf = sqrt(lift.^2 + drag.^2);
adt = taf/msc;
alt = r-rmars;
gs = adt/9.8;
tmin = t/60;
thrs = t/3600;
alpha = alpha*180/pi;
rvec = [x y zeros(length(x),1)];
vvec = [vx vy zeros(length(vx),1)];
h = cross(rvec,vvec);
hmag = sqrt(h(:,1).^2 + h(:,2).^2 + h(:,3).^2);
scripte = v.^2/2 + umars./r;
ecc = sqrt(1 + 2*scripte.*hmag.^2/umars^2);
[e f] = max(alt);
alts = alt(1:f);
[e f] = min(alts);
alts = alts(f:length(alts));
[a100err c100] = min(abs(alts-100));
c100 = c100+f-1;
rvecexit = [x(c100) y(c100) 0];
vvecexit = [vx(c100) vy(c100) 0];
hvecexit = cross(rvecexit,vvecexit);
hexitmag = sqrt(hvecexit(1)^2 + hvecexit(2)^2 + hvecexit(3)^2);
gammaexit = acos(hexitmag/(r(c100)*v(c100)));
gammaexit = gammaexit*180/pi;
vexit = v(c100);
theta = linspace(0,2*pi,100);
xms = rmars.*cos(theta);
yms = rmars.*sin(theta);
figure
plot(x,y)
hold on
axis equal
fill(xms,yms,'r')
if contact == 1
plot(x(b),y(b),'g*')
end
aae450atrajt.m – code responsible for propagating the Habitation Module’s arrival into the Martian atmosphere
%Nicholas Saadah
%AAE 450
%Presentation #1 (1/16/01)
%Hab arrival at Mars
clear all
global fcflag fhflag faflag feflag
%Constants
umars = -42828; %km^3/s^2
g = 9.8; %m/s^2
rmars = 3397; %km
msc = 75000; %kg
areasc = 65; %m^2
fcflag = 0;
fhflag = 0;
faflag = 0;
feflag = 0;
%Initial Conditions
xi = rmars+100; %km
yi = 0; %km
vi = 8.39; %km/s
gammaideg= input('gamma: '); %deg
maxgload = input('maxgload: ');
betadeg = input('beta: '); %deg
gammai = gammaideg*pi/180; %radians
beta = betadeg*pi/180; %radians
vxi = vi*sin(gammai); %km/s
vyi = vi*cos(gammai); %km/s
p0 = [xi vxi yi vyi];
ti = [0 200];
options = odeset('Reltol',1e-8);
[t,xx] = ode45('aae450eomt',ti,p0,options,umars,rmars,msc,areasc,g,beta,maxgload);
x = xx(:,1);
vx = xx(:,2);
y = xx(:,3);
vy = xx(:,4);
r = sqrt(x.^2 + y.^2);
alt = r-rmars;
clapp = min(alt);
if min(alt) < 0
contact = 1;
[a b] = min(abs(alt));
x = x(1:b);
vx = vx(1:b);
y = y(1:b);
vy = vy(1:b);
t = t(1:b);
r = r(1:b);
alt = alt(1:b);
else
contact = 0;
end
if fhflag == 1
[c d] = size(xx);
for i = 1:1:c-1
if vx(i+1) - vx(i) == 0
htrj = i;
break
end
end
x = x(1:htrj);
vx = vx(1:htrj);
y = y(1:htrj);
vy = vy(1:htrj);
t = t(1:htrj);
r = r(1:htrj);
alt = alt(1:htrj);
end
v = sqrt(vx.^2 + vy.^2);
ro = .00047*exp(-.1*((r-rmars)-49));
q = .5.*ro.*(v.*1000).^2;
q = q + 1e-12;
k = .125*((msc*maxgload*g)./(q*areasc)).^2;
%options = optimset('Display','off');
for i = 1:1:length(q)
if k(i) > 2
alpha(i) = pi/2;
elseif k(i) > 0 & k(i) < 2
alpha(i) = fzero('aae450fzeros',[0 pi/2],[],[],k(i));
end
end
alpha = alpha';
cl = 2*sin(2*alpha);
lift = cl*areasc.*sin(alpha).*q;
cd = 2*(1-cos(2*alpha));
drag = cd*areasc.*sin(alpha).*q;
taf = sqrt(lift.^2 + drag.^2);
adt = taf/msc;
alt = r-rmars;
gs = adt/9.8;
tmin = t/60;
thrs = t/3600;
alpha = alpha*180/pi;
rvec = [x y zeros(length(x),1)];
vvec = [vx vy zeros(length(vx),1)];
h = cross(rvec,vvec);
hmag = sqrt(h(:,1).^2 + h(:,2).^2 + h(:,3).^2);
scripte = v.^2/2 + umars./r;
ecc = sqrt(1 + 2*scripte.*hmag.^2/umars^2);
[e f] = max(alt);
alts = alt(1:f);
[e f] = min(alts);
alts = alts(f:length(alts));
[a100err c100] = min(abs(alts-100));
c100 = c100+f-1;
rvecexit = [x(c100) y(c100) 0];
vvecexit = [vx(c100) vy(c100) 0];
hvecexit = cross(rvecexit,vvecexit);
hexitmag = sqrt(hvecexit(1)^2 + hvecexit(2)^2 + hvecexit(3)^2);
gammaexit = acos(hexitmag/(r(c100)*v(c100)));
gammaexit = gammaexit*180/pi;
vexit = v(c100);
theta = linspace(0,2*pi,100);
xms = rmars.*cos(theta);
yms = rmars.*sin(theta);
figure
plot(x,y)
hold on
axis equal
fill(xms,yms,'r')
if contact == 1
plot(x(b),y(b),'g*')
end
aae450eoms.m – equations of motion for Earth arrival
%Nicholas Saadah
%AAE 450
%Planetary EOMs
%CTV arrival to Earth
function xdot = aae450eoms(t,x,flag,umars,rmars,msc,areasc,g,beta,maxgload)
global fcflag fhflag faflag feflag fvflag
r = sqrt(x(1)^2 + x(3)^2);
rvec = [x(1) x(3) 0];
v = sqrt(x(2)^2 + x(4)^2);
vvec = [x(2) x(4) 0];
ro = .047*exp(-.1*((r-rmars)-49));
q = .5*ro*(v*1000)^2;
if q < 1e-12
k = 3;
else
k = .125*((msc*maxgload*g)/(q*areasc))^2;
end
%options = optimset('Display','off');
if k > 2
alpha = pi/2;
elseif k > 0 & k <= 2
alpha = fzero('aae450fzeros',[0 pi/2],[],[],k);
else
fprintf('error\n')
end
cd = 2*(1-cos(2*alpha));
cl = 2*sin(2*alpha);
drag = cd*areasc*sin(alpha)*q;
lift = cl*areasc*sin(alpha)*q*cos(beta);
aero = sqrt(lift^2 + drag^2);
gload = aero/(msc*g);
h = cross(rvec,vvec);
hmag = sqrt(h(1)^2 + h(2)^2 + h(3)^2);
scripte = v^2/2 + umars/r;
ecc = sqrt(1 + 2*scripte*hmag^2/umars^2);
xdot(1) = x(2);
xdot(2) = umars*x(1)/r^3 - (.001*drag*x(2))/(v*msc) + ...
(.001*lift*x(4))/(v*msc);
xdot(3) = x(4);
xdot(4) = umars*x(3)/r^3 - (.001*drag*x(4))/(v*msc) - ...
(.001*lift*x(2))/(v*msc);
if r-rmars <= 100 & faflag == 0
fprintf('atmospheric entry\n')
faflag = 1;
end
if r-rmars > 100 & faflag == 1
fprintf('atmospheric exit\n')
faflag = 0;
end
if r < rmars & fcflag == 0
fprintf('contact with planet surface\n')
fcflag = 1;
end
if r-rmars > 200 & x(3) > 0 & ecc > 1 & fhflag == 0
fprintf('hyperbolic trajectory - e = %.4f\n',ecc)
fhflag = 1;
end
if r-rmars > 200 & x(3) > 0 & ecc < 1 & feflag == 0
fprintf('elliptic trajectory - e = %.4f\n',ecc)
feflag = 1;
end
if v <= 1
fvflag = 1;
end
if r < rmars | fhflag == 1 | v <= 1
xdot(1) = 0;
xdot(2) = 0;
xdot(3) = 0;
xdot(4) = 0;
end
xdot = xdot';
return
aae450eomt.m – equations of motion for Mars arrival
%Nicholas Saadah
%AAE 450
%Planetary EOMs
%CTV arrival to Earth
function xdot = aae450eomt(t,x,flag,umars,rmars,msc,areasc,g,beta,maxgload)
global fcflag fhflag faflag feflag
r = sqrt(x(1)^2 + x(3)^2);
rvec = [x(1) x(3) 0];
v = sqrt(x(2)^2 + x(4)^2);
vvec = [x(2) x(4) 0];
ro = .00047*exp(-.1*((r-rmars)-49));
q = .5*ro*(v*1000)^2;
if q < 1e-12
k = 3;
else
k = .125*((msc*maxgload*g)/(q*areasc))^2;
end
%options = optimset('Display','off');
if k > 2
alpha = pi/2;
elseif k > 0 & k <= 2
alpha = fzero('aae450fzeros',[0 pi/2],[],[],k);
else
fprintf('error\n')
end
cd = 2*(1-cos(2*alpha));
cl = 2*sin(2*alpha);
drag = cd*areasc*sin(alpha)*q;
lift = cl*areasc*sin(alpha)*q*cos(beta);
aero = sqrt(lift^2 + drag^2);
gload = aero/(msc*g);
h = cross(rvec,vvec);
hmag = sqrt(h(1)^2 + h(2)^2 + h(3)^2);
scripte = v^2/2 + umars/r;
ecc = sqrt(1 + 2*scripte*hmag^2/umars^2);
xdot(1) = x(2);
xdot(2) = umars*x(1)/r^3 - (.001*drag*x(2))/(v*msc) + ...
(.001*lift*x(4))/(v*msc);
xdot(3) = x(4);
xdot(4) = umars*x(3)/r^3 - (.001*drag*x(4))/(v*msc) - ...
(.001*lift*x(2))/(v*msc);
if r-rmars <= 100 & faflag == 0
fprintf('atmospheric entry\n')
faflag = 1;
end
if r-rmars > 100 & faflag == 1
fprintf('atmospheric exit\n')
faflag = 0;
end
if r < rmars & fcflag == 0
fprintf('contact with planet surface\n')
fcflag = 1;
end
if r-rmars > 200 & x(3) > 0 & ecc > 1 & fhflag == 0
fprintf('hyperbolic trajectory - e = %.4f\n',ecc)
fhflag = 1;
end
if r-rmars > 200 & x(3) > 0 & ecc < 1 & feflag == 0
fprintf('elliptic trajectory - e = %.4f\n',ecc)
feflag = 1;
end
if r < rmars | fhflag == 1
xdot(1) = 0;
xdot(2) = 0;
xdot(3) = 0;
xdot(4) = 0;
end
xdot = xdot';
return
aae450fzeros.m – flat plate equation
%Nicholas Saadah
%AAE 450 fzero solver
%2/17/01
function alpha = aae450fzero(x,k)
alpha = [(sin(x))^2*(1-cos(2*x)) - k];
return