function RocketSim(NonlinearFlag)

% Hovering rocket simulation with output feedback control.
% NonlinearFlag says whether to run the nonlinear simulation or the linear simulation.

if ~exist('NonlinearFlag', 'var')
    NonlinearFlag = true;
end

G = 6.673e-11; % gravitational constant (m^3/kg/s^2)
M = 5.98e24; % Earth mass (kg)
M = 7.35e22; % Moon mass (kg)
m0 = 32000; % initial rocket mass (kg) (lunar excursion module)
R = 6.37e6; % Earth radius (meters)
R = 1.74e6; % Moon radius (meters)
g = 50; % drag constant (kg/s)
mu = 1000; % rocket thrust proportionality constant (m/s)

% Get the linearized A, B, C, D matrices for the plant.
GA = [0 1 0; 2*G*M/R^3 -g/m0 -G*M/R^2/m0; 0 0 0];
GB = [0; mu/m0; -1];
GC = [1 0 0];
GD = [0];
% Get the transfer function of the plant.
[numG, denG] = ss2tf(GA, GB, GC, GD);
% Define the controller transfer function.
% The controller was designed using the root locus and the Bode plot.
numK = 3.2 * conv([1 0.001],[1 0.01]);
denK = conv([1 0],[1 1]);
% Get the open loop controller-plus-plant transfer function.
numOL = conv(numG, numK);
denOL = conv(denG, denK);

xnom = [0; 0; m0]; % nominal initial state
ynom = xnom(1); % nominal initial output
unom = G * M * xnom(3) / mu / R^2; % nominal control
if NonlinearFlag
    x = xnom; % initial nonlinear plant state
    y = GC * x; % initial nonlinear plant output
    dy = y - ynom; % initial delta-y
else
    dx = [0; 0; 0]; % initial linearized plant state
    dy = GC * dx; % initial linearized plant output
end

% Get the A, B, C, D matrices for the controller.
Ksys = tf(numK, denK);
[KA, KB, KC, KD] = ssdata(Ksys);
nK = size(KA, 1); % number of controller states
xK = zeros(nK, 1); % initial controller state vector

dt = 0.001; % simulation step size
tSim = 100; % simulation length
dtPlot = 0.1; % how often to plot data
toldPlot = -inf;

% Initialize arrays for saving plot data.
tPlot = [0];
if NonlinearFlag
    xPlot = x;
    yPlot = y;
    uPlot = unom;
else
    xPlot = xnom + dx;
    yPlot = ynom + dy;
    uPlot = unom + 0;
end

for t = dt : dt : tSim
    % nominal trajectory
    xnom = [0; 0; m0*exp(-G*M*t/mu/R^2)];
    unom = G * M * xnom(3) / mu / R^2; 
    % Define the reference input
    r = 10; % hover at a constant altitude
    %r = t; % rise at a constant rate
%    omega = 10; r = sin(omega * t); % track a sinusoidal altitude
    du = r - dy; % controller input
    xdot = KA * xK + KB * du;
    xK = xK + xdot * dt;
    dy = KC * xK + KD * du; % controller output
    du = dy; % linearized plant input    
    if NonlinearFlag
        % Nonlinear system simulation
        % Plant input must not be negative
        u = unom + du;
        if u < 0
            u = 0;
        end 
        xdot(1) = x(2);
        xdot(2) = (mu * u - g * x(2)) / x(3) - G * M / (R + x(1))^2;
        xdot(3) = -u;
        x = x + xdot * dt; % updated plant state
        % Altitude must not be negative
        if x(1) < 0
            x(1) = 0;
        end
        y = x(1); % plant output
        dy = y - ynom; % compute delta-y for the controller
    else
        % Linearized system simulation
        % Plant input must not be negative
        if unom + du < 0
            du = -unom;
        end
        xdot = GA * dx + GB * du;
        dx = dx + xdot * dt; % updated plant state
        % Altitude must not be negative
        if xnom(1) + dx(1) < 0
            dx(1) = -xnom(1);
        end
        dy = GC * dx + GD * du; % plant output
    end
    if t > toldPlot + dtPlot
        toldPlot = t;
        tPlot = [tPlot t];
        if NonlinearFlag
            xPlot = [xPlot x];
            yPlot = [yPlot y];
            uPlot = [uPlot u];
        else
            xPlot = [xPlot xnom+dx];
            yPlot = [yPlot ynom+dy];
            uPlot = [uPlot unom+du];
        end
   end
end
% Plot results
close all;
figure;
subplot(2,2,1);
plot(tPlot, uPlot); legend('control (kg/s)');
if NonlinearFlag
    title('Nonlinear simulation');
else
    title('Linear simulation');
end
subplot(2,2,2);
plot(tPlot, xPlot(1,:)); legend('altitude (m)');
subplot(2,2,3);
plot(tPlot, xPlot(2,:)); legend('velocity (m/s)');
xlabel('time (seconds)');
subplot(2,2,4);
plot(tPlot, xPlot(3,:)); legend('mass (kg)');
xlabel('time (seconds)');