gravturn

%.............................................................
%
%   Gravity Turn Terminal Descent Algorithm
%
%............................................................

%..Initialize global 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 m_trans
global t_trans
global i_trans

global dist
global ts

%..Inputs passed from Trajectory Simulation..................

v(1)   = vs;                %..Initial Velocity, m/s
m(1)   = ms;                %..Initial Mass, kg
h(1)   = hs;                %..Initial Altitude, m
gam(1) = -gs;               %..Initial Flight Path Angle, deg
g      = 1.62;              %..Moon gravity constant, m/s^2
ge     = 9.81;              %..Earth gravity constant, m/s^2

%..Set a cutoff altitude if other than zero
        hov    = 0;         %..Cutoff altitude, m

%..Iteration values..........................................

i      = 1;
dt     = 0.05;

%..Gravity Turn..............................................

R(i)   = -(h(i)-hov)/sin(gam(i));
x(i)   = 0;
t(i)   = 0;

T(i)   = m(i)*(g + v(i)^2/(2*R(i)));

%..Output initial gravity turn parameters....................


disp(sprintf('Conditions at Start of Gravity Turn'))
disp(sprintf('-------------------------------------------------'))
disp(sprintf('                                                 '))
disp(sprintf('Velocity              %g    km/s',v(1)/1000))
disp(sprintf('Altitude              %g    km',h(1)/1000))
disp(sprintf('FPA                   %g   deg',gam(1)*180/pi))
disp(sprintf('Slant Range           %g     km',R(1)/1000))
disp(sprintf('Thrust                %g    N',T(1)))
disp(sprintf('                                                 '))
disp(sprintf('-------------------------------------------------'))
disp(sprintf('                                                 '))

while h(i) > hov

    dv   = -g*(T(i)/(m(i)*g)+sin(gam(i)));
    dgam = (-g/v(i))*cos(gam(i));
    dh   = v(i)*sin(gam(i));
    dx   = v(i)*cos(gam(i));
    dm   = -T(i)/(ge*isp);

    v(i+1)   = v(i) + dv*dt;
    gam(i+1) = gam(i) + dgam*dt;
    h(i+1)   = h(i) + dh*dt;
    x(i+1)   = x(i) + dx*dt;
    m(i+1)   = m(i) + dm*dt;
    R(i+1)   = -(h(i+1)-hov)/sin(gam(i+1));
    T(i+1)   = m(i+1)*(g + v(i+1)^2/(2*R(i+1)));
    n(i+1)   = T(i+1)/(m(i+1)*g);

    t(i+1)   = t(i) + dt;

    i = i + 1;
end

%..Store data to transfer to engine analysis

m_trans = m;
t_trans = t;
i_trans = i;

%..Output Gravity Turn Summary Data....................................

disp(sprintf('Gravity Turn Summary Data'))
disp(sprintf('-------------------------------------------------'))
disp(sprintf('                                                 '))
disp(sprintf('Grav turn duration    %g      sec',t(i)))
disp(sprintf('Total entry time      %g      sec',t(i)+ts))
disp(sprintf('Grav turn dist        %g    km',x(i)/1000))
disp(sprintf('Total entry dist      %g    km',(x(i)+dist)/1000))
disp(sprintf('Final FPA             %g   deg', gam(i)*180/pi))
disp(sprintf('Propellant mass used  %g    kg',m(1)-m(i)))
disp(sprintf('-------------------------------------------------'))

if m(1)-m(i) > mf
    disp('***WARNING: You need more terminal descent fuel!***')
else
    disp('***NOTE: You have excess terminal descent fuel!***')
end

%..Plot Data...........................................................

figure

subplot(2,3,1)
plot(t,h);
xlabel('Time, s');
ylabel('Altitude, m')
title('Altitude vs. Time')

subplot(2,3,2)
plot(t,gam*180/pi);
xlabel('Time, s');
ylabel('Flight Path Angle, deg')
title('Flight Path Angle vs. Time')

subplot(2,3,3)
plot(t,v);
xlabel('Time, s');
ylabel('Velocity, m/s');
title('Velocity vs Time')

subplot(2,3,4)
plot(t,T);
xlabel('Time, s');
ylabel('Thrust, N');
title('Thurst vs. Time');

subplot(2,3,[5 6])
plot(t,n);
xlabel('Time, s');
ylabel('Load Factor');
title('Load Factor vs. Time');

%..Call engine analysis

engine;