%.............................................................
%
%   Lunar Terminal Landing Simulation
%
%   Assume:
%       -Flat planet for Gravity Turn portion for sim
%       -Retrograde initial orbit (if prograde, change w2 def)
%
%.............................................................

clear all
clc

%..Initialize gobal variables

global vs
global ms
global hs
global gs
global mf

global game
global pe
global pa
global p0
global Arat
global Tav
global nang
global isp
global Cstar
global dist
global ts


%..Input initial state.........................................

md  = 67.692;            %..Dry mass, kg
mf  = 47.5;            %..Fuel mass, kg
ro  = 100;               %..Orbital radius, km

dv = -685;           %..Deorbit delta-V, m/s

tdc = 70000;             %..Terminal Descent Ceiling, m
                         %  (Must begin GT by this altitude)

%..Engine Inputs...............................................

game  = 1.2;              %..Approx. ratio of specific heats
pe    = 0.02*101325;      %..Nozzle exit pressure, Pa
pa    = 0;                %..Ambient pressure, Pa
p0    = 20*101325;        %..Tank pressure, Pa
Arat  = 50;               %..Nozzle area ratio
Tav   = 400;              %..Required max thrust, N
nang  = 15*pi/180;        %..Nozzle angle
isp   = 255;              %..Specific impulse, s
Cstar = 1350;             %..Charistic Velocity, m/s


%..............................................................
%..DO NOT EDIT BELOW THIS LINE.................................
%..............................................................
%..Constants...................................................

mu = 4.9028E12;         %..Moon Grav Parameter, m^3/s^2
rm = 1737400;           %..Radius of Moon, m
g0 = 1.6;               %..Moon Grav acceleration, m/s^2
wm = 2.661707223E-6;    %..Rotational Angular Velocity of Moon, rad/s

%..Calculate initial values for trajectory simulation...........

h0   = ro*1000;         %..Altitude above lunar surface, m
gam0 = 0;               %..Flight path angle, rad

v0   = sqrt(mu/rm);     %..Orbital Velocity (INERTIAL), m/s
w0   = v0/(rm+h0);      %..Angular Velocity (INERTIAL), rad/s
w2   = w0+wm;           %..Angular Velocity (RELATIVE), rad/s
v2   = w2*(rm+h0);      %..Orbital Velocity (RELATIVE), m/s
vdv  = v2 + dv;         %..Initial Ballistic Velocity (RELATIVE), m/s

%..Output Initial State Data...................................

disp(sprintf('-------------------------------------------------'))
disp(sprintf('LUNAR TERMINAL LANDING SIMULATION'))
disp(sprintf('-------------------------------------------------'))
disp(sprintf('                                                 '))
disp(sprintf('Initial State Data'))
disp(sprintf('-------------------------------------------------'))
disp(sprintf('                                                 '))
disp(sprintf('Dry Mass              %g     kg',md))
disp(sprintf('Fuel Mass             %g       kg',mf))
disp(sprintf('Total Entry Mass      %g    kg',md+mf))
disp(sprintf('Orbital Radius        %g        km',ro))
disp(sprintf('Rel Orbital Velocity  %g    km/s',v2/1000))
disp(sprintf('Deorbit Delta-V       %g     km/s',dv/1000))
disp(sprintf('Ballistic Velocity    %g   km/s',vdv/1000))
disp(sprintf('Terminal Thrust Avail %g        N',Tav))
disp(sprintf('                                                 '))
disp(sprintf('-------------------------------------------------'))
disp(sprintf('                                                 '))


%..Call ODE45 and integrate through ballistic trajectory

run = 1;
tstop = 1;

options = odeset('RelTol', 1E-10);

while run > 0

   [t,y] = ode45(@eom,[0:1:tstop],[vdv gam0 h0],options);

   v   = y(:,1);
   gam = y(:,2);
   h   = y(:,3);

   i = length(h);

   if h(i) < 0
       break
   end

   tstop = tstop + 1;

end

%..Sift through trajectory data to determine gravity turn initialization
%..point i.e. where Treq = Tav

run2 = 1;

for k = 2:length(h)
   R(k)       = h(k)/sin(gam(k));
   Treq(k)    = (md+mf)*(g0+v(k)^2/(2*R(k)));
   if Treq(k) > Tav
       break
   elseif h(k) < tdc
       break
   end
end

%..Output state data to gravity turn algorithm

vs = v(k);
ms = md+mf;
hs = h(k);
gs = gam(k);
ts = t(k);

%..Calculate X distance transversed

vx     = v.*cos(gam);
dist   = trapz((vx(1:k)))*1;

%..Call gravity turn algorithm

gravturn


%..Plot Ballistic Trajectory..................

% figure
% plot(gam*180/pi,h)
% title('Altitude vs. Flight Path Angle')
% xlabel('Flight Path Angle, deg')
% ylabel('Altitude, m')
%
% figure
% plot(v,h)
% title('Altitude vs. Velocity')
% xlabel('Velocity, m/s')
% ylabel('Altitude, m/s')
%
% figure
% plot(t,h)
% title('Altitude vs. Time')
% xlabel('Time, s')
% ylabel('Altitude, m')

Published with MATLAB® 7.4