Contents

function [ simOut, passFail ] = RunF16Sim(initialState, t_vec, ...
    orient,F16_model,flightLimits,ctrlLimits,autopilot, analysisOn, ...
    printOn, plotOn)
%Simulates and analyzes autonomous F-16 maneuvers using ode45
%
%   Function Calls:
%       [ simOut ] = NonlinSim_F16_func(initialState, t_vec, orient, ...
%                       plantModel, printOn, plotOn)
%       [ simOut, passFail ] = NonlinSim_F16_func(initialState, ...
%           0:0.01:15, 4, 'morelli', true, true)
%
%   Inputs:
%       initialState - initial F-16 state vector in feet, rad, etc.
%                       [Vt alpha beta phi theta psi P Q R pn pe h pow]
%       t_vec    - desired time vector (1xN)
%       orient   - desired trim condition (values 1, 2, 3, or 4)
%                       1:    Wings Level (gamma = 0)
%                       2:    Wings Level (gamma <> 0)
%                       3:    Constant Altitude Turn
%                       4:    Steady Pull Up
%       F16_model - Optional string defining which function to use to
%           calculate state derivatives: {'morelli','stevens','linear'}.
%           Defaults to 'morelli'.
%       flightLimits - struct of pass fail limits
%       ctrlLimits - struct of control limits
%       analysisOn - if true, runs full analysis on sim results. If false,
%           only returns vector as output
%       printOn  - if true, prints results to console
%       plotOn   - if true, plots results in figures
%
%   Outputs:
%       simOut  - structure of analyzed results
%
%   Comments:
%   This is a functionalized version of "NonlinSim_F16", meant to be used
%   for monte-carlo style simulations or analysis if needed. All of the
%   calculations are the same. If this function is called with no
%   arguments, it runs my default test case.
%
%   The outer-loop control methods implemented are contained in the
%   function "getOuterLoopCtrl". If you desire to use hard-coded outer-loop
%   control or maneuvers, or want to change the implementation of GCAS,
%   that is the function which needs to be changed.
%
% <a href="https://github.com/pheidlauf/AeroBenchVV">AeroBenchVV</a>
% Copyright: GNU General Public License 2017
%
% See also: CONTROLLEDF16, GETAUTOPILOTCOMMANDS, ODE45

% If called with no arguments, run "Main.m".
if(nargin==0)
    Main;
end

% Check that enough arguments are given
if(nargin<7)
    error('Insufficient Arguments given')
end

% Default analysisOn to true
if(~exist('analysisOn','var'))
    analysisOn = true;
end
% Default printOn to true
if(~exist('printOn','var'))
    printOn = true;
end

% Default plotOn to true
if(~exist('plotOn','var'))
    plotOn = true;
end

% If true, runs the profiler for streamlining & debugging
profilerOn = false;
% Start the clock to track processing time
timerStart = tic;
% Add necessary sub folder to path
addpath(genpath('F16_Model'));
flightLimits = 

      altitudeMin: 0
      altitudeMax: 10000
     maneuverTime: 15
            NzMax: 9
            NzMin: -2
    psMaxAccelDeg: 500
             vMin: 300
             vMax: 900
      alphaMinDeg: -10
      alphaMaxDeg: 45
       betaMaxDeg: 30


ctrlLimits = 

       ThrottleMax: 1
       ThrottleMin: 0
    ElevatorMaxDeg: 25
    ElevatorMinDeg: -25
     AileronMaxDeg: 21.500000000000000
     AileronMinDeg: -21.500000000000000
      RudderMaxDeg: 30
      RudderMinDeg: -30
        MaxBankDeg: 60
             NzMax: 6
             NzMin: -1


autopilot = 

                    title: 'Default Simulation'
        basicSpeedControl: 1
    steadyLevelFlightHold: 0
         levelTurnControl: 0
               simpleGCAS: 0
            turnToHeading: 0
     timeTriggeredControl: 0

Check flightLimits against model requirements

Limit airspeed from 300 to 900 ft/s

if(flightLimits.vMin < 300 || flightLimits.vMax > 900)
    error('flightLimits: Airspeed limits outside model limits (300 to 900)');
end
% Limit alpha from -10 to 45 deg
if(flightLimits.alphaMinDeg < -10 || flightLimits.alphaMaxDeg > 45)
    error('flightLimits: Alpha limits outside model limits (-10 to 45)');
end
% Limit beta from -30 to 30 deg (positive only, symmetric)
if(abs(flightLimits.betaMaxDeg) > 30)
    error('flightLimits: Beta limit outside model limits (30 deg)');
end

Check ctrlLimits against model limits (user-proofing)

Limit throttle from 0 to 1

if(ctrlLimits.ThrottleMin < 0 || ctrlLimits.ThrottleMax > 1)
    error('ctrlLimits: Throttle Limits (0 to 1)');
end

% Limit elevator from -25 to 25 deg
if(ctrlLimits.ElevatorMaxDeg > 25 || ctrlLimits.ElevatorMinDeg < -25)
    error('ctrlLimits: Elevator Limits (-25 deg to 25 deg)');
end

% Limit aileron from -21.5 to 21.5 deg
if(ctrlLimits.AileronMaxDeg > 21.5 || ctrlLimits.AileronMinDeg < -21.5)
    error('ctrlLimits: Aileron Limits (-21.5 deg to 21.5 deg)');
end

% Limit rudder from -30 to 30 deg
if(ctrlLimits.RudderMaxDeg > 30 || ctrlLimits.RudderMinDeg < -30)
    error('ctrlLimits: Rudder Limits (-30 deg to 30 deg)');
end

Get Trim / Equilibrium Conditions

if(printOn)
    disp('------------------------------------------------------------');
    disp('F-16 Decoupled LQR Controller for Nz, P_s, and Ny+r tracking');
    disp('------------------------------------------------------------');
end

% Format initial Conditions for Simulation and append integral error states
x0 = [initialState'; 0; 0; 0];

% Define Control Guess
uguess = [.2 0 0 0];

% Format inputs for trimmerFun
inputs = [initialState(1), initialState(12), 0, 0, 0];

if(printOn)
    printmat(inputs,'Operator Inputs',[],'Vt h gamma psidot thetadot')
    fprintf('Trim Orientation Selected:   ');
    switch orient
        case 1
            disp('Wings Level (gamma = 0)');
        case 2
            disp('Wings Level (gamme <> 0)');
        case 3
            disp('Constant Altitude Turn');
        case 4
            disp('Steady Pull Up');
        otherwise
            error('Invalid Orientation for trimmerFun');
    end
    printmat(initialState,'Inititial Conditions',[],...
        'Vt alpha beta phi theta psi p q r pn pe alt pow');
    printmat(uguess,'Control Guess',[],...
        'throttle elevator aileron rudder');
end

% Get Equilibrium Values
[xequil,uequil] = trimmerFun(initialState, uguess, orient, inputs, false);

if(printOn)
    disp('------------------------------------------------------------');
    disp('Equilibrium / Trim Conditions');
    printmat(xequil','State Equilibrium',[],...
        'Vt alpha beta phi theta psi p q r pn pe alt pow');
    printmat(uequil','Control Equilibrium',[],...
        'throttle elevator aileron rudder');
end
------------------------------------------------------------
F-16 Decoupled LQR Controller for Nz, P_s, and Ny+r tracking
------------------------------------------------------------
 
Operator Inputs = 
                        Vt            h        gamma       psidot     thetadot
                 540.00000   3500.00000            0            0            0
 
Trim Orientation Selected:   Steady Pull Up
 
Inititial Conditions = 
                        Vt        alpha         beta          phi        theta
                 540.00000      0.03703            0      0.78540     -1.25664
 
                       psi            p            q            r           pn
                  -0.78540            0            0            0            0
 
                        pe          alt          pow
                         0   3500.00000      9.00000
 
 
Control Guess = 
                  throttle     elevator      aileron       rudder
                   0.20000            0            0            0
 
------------------------------------------------------------
Equilibrium / Trim Conditions
 
State Equilibrium = 
                        Vt        alpha         beta          phi        theta
                 540.00000      0.03441            0            0      0.03441
 
                       psi            p            q            r           pn
                         0            0            0            0            0
 
                        pe          alt          pow
                         0   3500.00000     10.35403
 
 
Control Equilibrium = 
                  throttle     elevator      aileron       rudder
                   0.15944     -0.77048            0            0
 

Build Linear LQR Controller from Workspace Data

% Get linearized ss-model (in case F16_method=='linear');
lin_f16 = getLinF16(xequil,uequil,printOn);

% Load Decoupled Controllers from workspace
load('lateralCtrlData.mat')
load('longitudinalCtrlData.mat')

% Append decoupled LQR gain matrices
K_lqr = blkdiag(K_long, K_lat);
if(printOn)
    printmat(K_lqr,'Decoupled LQR Controller Gains', ...
        'elevator aileron rudder',...
        'alpha q Nz_i beta p r ps_i Ny_r_i');
end
------------------------------------------------------------
Running jacobFun.m
Linearized F-16 SS Model

lin_f16 =
 
  A = 
                  Vt       alpha        beta         phi       theta
   Vt       -0.01658        10.6  -8.377e-06  -2.766e-06      -32.17
   alpha  -0.0002195     -0.9855   3.055e-06  -1.488e-06  -1.025e-05
   beta            0           0     -0.3124     0.05954           0
   phi             0           0           0           0           0
   theta           0           0           0           0           0
   psi             0           0           0           0           0
   p               0           0       -31.6           0           0
   q       6.662e-13       0.858           0           0           0
   r               0           0       8.904           0           0
   pn              1  -6.607e-06           0  -1.592e-06  -6.607e-06
   pe              0           0         540      -18.58           0
   alt    -6.579e-16        -540           0    4.64e-06         540
   pow             0           0           0           0           0
 
                 psi           p           q           r          pn
   Vt              0           0     -0.4842           0           0
   alpha           0           0      0.9146           0           0
   beta            0     0.03387           0     -0.9925           0
   phi             0           1           0     0.03443           0
   theta           0           0           1           0           0
   psi             0           0           0       1.001           0
   p               0      -3.574   0.0002627      0.6331           0
   q               0  -8.795e-06      -1.045   -0.002858           0
   r               0    -0.02249    0.002539     -0.4621           0
   pn              0           0           0           0           0
   pe            540           0           0           0           0
   alt             0           0           0           0           0
   pow             0           0           0           0           0
 
                  pe         alt         pow
   Vt              0    3.81e-05      0.3717
   alpha           0   1.774e-06   -2.37e-05
   beta            0           0           0
   phi             0           0           0
   theta           0           0           0
   psi             0           0           0
   p               0           0           0
   q               0  -5.338e-15           0
   r               0           0           0
   pn              0           0           0
   pe              0           0           0
   alt             0           0           0
   pow             0           0          -1
 
  B = 
           Throttle   Elevator    Aileron     Rudder
   Vt             0     0.1848          0          0
   alpha          0  -0.002085          0          0
   beta           0          0  0.0002863  0.0007816
   phi            0          0          0          0
   theta          0          0          0          0
   psi            0          0          0          0
   p              0          0     -0.765     0.1376
   q              0    -0.1833          0          0
   r              0          0   -0.03332   -0.06472
   pn             0          0          0          0
   pe             0          0          0          0
   alt            0          0          0          0
   pow        64.94          0          0          0
 
  C = 
                  Vt       alpha        beta         phi       theta
   Az       0.003698        16.9  -5.296e-08           0           0
   q               0           0           0           0           0
   alpha           0           1           0           0           0
   theta           0           0           0           0           1
   Vt              1           0           0           0           0
   Ay              0           0       -5.24           0           0
   p               0           0           0           0           0
   r               0           0           0           0           0
   beta            0           0           1           0           0
   phi             0           0           0           1           0
 
                 psi           p           q           r          pn
   Az              0    -4.1e-08      0.9454   -0.001336           0
   q               0           0           1           0           0
   alpha           0           0           0           0           0
   theta           0           0           0           0           0
   Vt              0           0           0           0           0
   Ay              0   -0.008953           0      0.1154           0
   p               0           1           0           0           0
   r               0           0           0           1           0
   beta            0           0           0           0           0
   phi             0           0           0           0           0
 
                  pe         alt         pow
   Az              0  -2.979e-05           0
   q               0           0           0
   alpha           0           0           0
   theta           0           0           0
   Vt              0           0           0
   Ay              0           0           0
   p               0           0           0
   r               0           0           0
   beta            0           0           0
   phi             0           0           0
 
  D = 
          Throttle  Elevator   Aileron    Rudder
   Az            0  -0.05062         0         0
   q             0         0         0         0
   alpha         0         0         0         0
   theta         0         0         0         0
   Vt            0         0         0         0
   Ay            0         0  0.004801   0.01311
   p             0         0         0         0
   r             0         0         0         0
   beta          0         0         0         0
   phi           0         0         0         0
 
Name: Linearized F-16 SS Model
Continuous-time state-space model.

 
Decoupled LQR Controller Gains = 
                     alpha            q         Nz_i         beta            p
     elevator   -156.88015    -31.03701    -38.72983            0            0
      aileron            0            0            0     38.02751     -5.65497
       rudder            0            0            0     17.56400      1.58391
 
                         r         ps_i       Ny_r_i
     elevator            0            0            0
      aileron    -14.08804    -34.06416     -9.95406
       rudder    -41.43509      6.29550    -53.86016
 

Simulate Using ODE Solver

if(printOn)
    disp('------------------------------------------------------------');
    disp('Running Nonlinear Simulation in ODE45');
end

% USER EDIT: Set ode options if desired
options = [];

% Call ode45 to simulate controlled nonlinear system
[time,x_f16_hist] = ode45(@(t,y) controlledF16(t,y,xequil, uequil,...
    K_lqr, F16_model, lin_f16, flightLimits,ctrlLimits,autopilot),...
    t_vec, x0, options);

% Record time to complete simulation
t_simmed = toc(timerStart);

if(printOn)
    fprintf('\n');
    disp('Simulation Complete');
    disp('Time to simulate (MM:SS.ms):');
    fprintf('%s\n',datestr((t_simmed)/(24*60*60),'MM:SS.FFF'))
    fprintf('\n');
    disp('Begin back-calculation of controls');
end
------------------------------------------------------------
Running Nonlinear Simulation in ODE45

Simulation Complete
Time to simulate (MM:SS.ms):
00:00.731

Begin back-calculation of controls

Back-calculate & format needed intermediate states

Format Output for easier indexing

x_f16_hist = x_f16_hist';

if(~analysisOn)
    simOut = x_f16_hist;
    passFail = 'NO ANALYSIS COMPLETED';
    warning('No analysis completed');
    return
end

% Initialize "history" vectors
u_hist = zeros(7,length(time));         % Control
Nz_hist = zeros(1,length(time));        % Nz (full nonlinear calculation)
ps_hist = zeros(1,length(time));        % ps (full nonlinear calculation)
Ny_r_hist = zeros(1,length(time));      % Ny_r
lin_Nz_hist = zeros(1,length(time));    % Nz (linear approximation)
lin_Ny_r_hist = zeros(1,length(time));    % Ny (linear approximation)

% Recall GCAS time steps from persistent memory of getAutopilotCommands
% [~,t_maneuver] = getAutopilotCommands(time(end), x0, xequil, uequil,...
%     ctrlLimits, false);
[~,t_maneuver] = getAutopilotCommands(time(end), x0, xequil, uequil, ...
    flightLimits, ctrlLimits, autopilot, false);

% Reset Persistent Variables in getOuterLoopControl
% getAutopilotCommands(0, x0, xequil, uequil, ctrlLimits, true);
getAutopilotCommands(0, x0, xequil, uequil, ...
    flightLimits, ctrlLimits, autopilot, true);

%{
Note: Several important intermediate values are not outputted from ode45
automatically. Below, I use the outputted states and time vector to
recalculate the values and store them for data analysis and visualization.
Because ode45 integrates with an non-fixed time-step, the recalled &
recalculated results may be slightly different. This primarily effects the
ps_ref and Nz_ref history.
%}

% Back calculate controlled states for each time step
for i = 1:length(time)
    [~,u,Nz,ps,Ny_r] = controlledF16(time(i),x_f16_hist(:,i),...
        xequil, uequil, K_lqr, F16_model, lin_f16,...
        flightLimits, ctrlLimits, autopilot);
    Nz_hist(i) = Nz;
    ps_hist(i) = ps;
    Ny_r_hist(i) = Ny_r;
    u_hist(:,i) = u;

    % Linear approx of Nz
    lin_Nz_hist(i) = (lin_f16.c(1,:)*(x_f16_hist(1:13,i) - xequil)+ ...
        lin_f16.d(1,:)*(u(1:4) + uequil));
    %     lin_Nz_hist(i) = (lin_f16.c(1,:)*(x_f16_hist(1:13,i)) + ...
    %         lin_f16.d(1,:)*(u(1:4)));

    % Linear approx of (Ny + r)
    lin_Ny_r_hist(i) = (lin_f16.c(6,:)*(x_f16_hist(1:13,i)-xequil)+...
        lin_f16.d(6,:)*(u(1:4) + uequil)) + x_f16_hist(9,i);
end

% Determine ps acceleration (rad/s/s)
ps_accels = zeros(1,length(ps_hist));
ps_accels(1) = 1;
for i = 2:length(ps_hist)
    ps_accels(i) = (ps_hist(i) - ps_hist(i-1))/(time(i) - time(i-1));
end

if(printOn)
    % Record processing time
    t_calced = toc(timerStart);
    disp('Back-Calculation of Controls Complete');
    disp('Time to calculate (MM:SS.ms):');
    fprintf('%s\n',datestr((t_calced-t_simmed)/(24*60*60),...
        'MM:SS.FFF'))
    fprintf('\n');
    disp('Begin struct generation & pass/fail analysis');
end
Back-Calculation of Controls Complete
Time to calculate (MM:SS.ms):
00:10.432

Begin struct generation & pass/fail analysis

Save Results to Struct

Create struct for state and error history

stateHistory = struct([]);
% F-16 Nonlinear Model States
stateHistory(1).time = t_vec;
stateHistory.VT = x_f16_hist(1,:);
stateHistory.alpha = x_f16_hist(2,:);
stateHistory.beta = x_f16_hist(3,:);
stateHistory.phi = x_f16_hist(4,:);
stateHistory.theta = x_f16_hist(5,:);
stateHistory.psi = x_f16_hist(6,:);
stateHistory.P = x_f16_hist(7,:);
stateHistory.Q = x_f16_hist(8,:);
stateHistory.R = x_f16_hist(9,:);
stateHistory.pn = x_f16_hist(10,:);
stateHistory.pe = x_f16_hist(11,:);
stateHistory.h = x_f16_hist(12,:);
stateHistory.pow = x_f16_hist(13,:);
stateHistory.x_f16_hist = x_f16_hist;
% Integral Error States
stateHistory.Nz_e_i = x_f16_hist(14,:);
stateHistory.ps_e_i = x_f16_hist(15,:);
stateHistory.Ny_r_e_i = x_f16_hist(16,:);
% Calculated outputs & controlled states
stateHistory.Nz_hist = Nz_hist;
stateHistory.ps_hist = ps_hist;
stateHistory.Ny_r_hist = Ny_r_hist;
stateHistory.lin_Nz_hist = lin_Nz_hist;
stateHistory.lin_Ny_hist = lin_Ny_r_hist;

% Create struct for back-calculated control history
controlHistory = struct([]);
controlHistory(1).throttle = u_hist(1,:);
controlHistory.elevator = u_hist(2,:);
controlHistory.aileron = u_hist(3,:);
controlHistory.rudder = u_hist(4,:);
controlHistory.Nz_ref = u_hist(5,:);
controlHistory.ps_ref = u_hist(6,:);
controlHistory.Ny_r_ref = u_hist(7,:);

% Create output struct to store results.
simOut = struct('states',stateHistory,'ctrls',controlHistory);

Numeric Data Analysis

Calculate Desired Performance Specs

minAltitude = min(x_f16_hist(12,:));            % ft
altitudeLost = x_f16_hist(12,1) - minAltitude;  % ft
maxAirspeed = max(x_f16_hist(1,:));             % ft/s
minAirspeed = min(x_f16_hist(1,:));             % ft/s
maxNz = max(Nz_hist);                           % g's (0==level flight)
minNz = min(Nz_hist);                           % g's
maxSideForce = max(abs(Ny_r_hist));             % g's
max_ps_accel = max(abs(ps_accels));             % rad/s/s

% Add analysis results to output struct.
simOut.ps_accels = ps_accels;
simOut.t_maneuver = t_maneuver;
simOut.x0 = x0;
simOut.trimState = xequil;
simOut.trimCtrl = uequil;
simOut.lin_ss_model = lin_f16;
simOut.plantModel = F16_model;
simOut.flightLimits = flightLimits;
simOut.controlLimits = ctrlLimits;

Calculate Pass/Fail Conditions

% Initialize output struct
passFail = struct([]);
passFail(1).stable = true;

% Check airspeed limits
if(maxAirspeed > flightLimits.vMax || minAirspeed < flightLimits.vMin)
    passFail.airspeed = false;
    passFail.stable = false;
else
    passFail.airspeed = true;
end

% Check alpha limits
if(max(x_f16_hist(2,:)) > deg2rad(flightLimits.alphaMaxDeg) || ...
        min(x_f16_hist(2,:)) < deg2rad(flightLimits.alphaMinDeg))
    passFail.alpha = false;
    passFail.stable = false;
else
    passFail.alpha = true;
end

% Check beta limits
if(abs(x_f16_hist(3,:)) > deg2rad(flightLimits.betaMaxDeg))
    passFail.beta = false;
    passFail.stable = false;
else
    passFail.beta = true;
end

% Check Nz limits
if(minNz < flightLimits.NzMin || maxNz > flightLimits.NzMax)
    passFail.Nz = false;
else
    passFail.Nz = true;
end

% Check Ps_rate limits
if(max_ps_accel > flightLimits.psMaxAccelDeg)
    passFail.psMaxAccelDeg = false;
else
    passFail.psMaxAccelDeg = true;
end

% Check altitude limits
if(min(x_f16_hist(12,:)) < flightLimits.altitudeMin ||...
        max(x_f16_hist(12,:)) > flightLimits.altitudeMax)
    passFail.altitude = false;
else
    passFail.altitude = true;
end

% Check maneuver time limits (GCAS?)
if(t_maneuver(2)-t_maneuver(1) > flightLimits.maneuverTime ||...
        t_maneuver(2) < 0)
    passFail.maneuverTime = false;
else
    passFail.maneuverTime = true;
end

% Add passFail to the main output too just in case
simOut.passFail = passFail;

if(printOn)
    % Record processing time
    t_analysis = toc(timerStart);
    disp('Pass/Fail Analysis Complete');
    disp('Time to calculate (MM:SS.ms):');
    fprintf('%s\n',datestr((t_analysis - t_calced)/(24*60*60),...
        'MM:SS.FFF'))
    fprintf('\n');
end
Pass/Fail Analysis Complete
Time to calculate (MM:SS.ms):
00:00.063

Print Values to Console

if(printOn)
    disp('----------------------------------------------------------');
    disp('RESULTS:');
    disp('              Min & Max Values:');
    fprintf('Min Altitude:          %9.3f   ft \n', minAltitude);
    fprintf('Altitude Lost:         %9.3f   ft \n', altitudeLost);
    fprintf('Min Airspeed:          %9.3f   ft/sec \n', minAirspeed);
    fprintf('Max Airspeed:          %9.3f   ft/sec \n', maxAirspeed);
    fprintf('Max Down Force:        %9.3f   g''s \n', maxNz);
    fprintf('Min Down Force:        %9.3f   g''s \n', minNz);
    fprintf('Max Side Force:        %9.3f   g''s \n', maxSideForce);
    fprintf('Max Roll Accel:        %9.3f   deg/sec^2 \n', ...
        rad2deg(max_ps_accel));
    fprintf('\n');
    disp('              Maneuver Event Times:');
    fprintf('Maneuver Start:        %9.3f   sec \n', t_maneuver(1));
    if(length(t_maneuver)>2)
        for(i = 3:(length(t_maneuver)))
            fprintf('Checkpoint %2d:         %9.3f   sec \n',...
                i-2,t_maneuver(i));
        end
    end
    fprintf('Maneuver Complete:     %9.3f   sec \n', t_maneuver(2));
    fprintf('Final Time:            %9.3f   sec \n', time(end));
    disp('(if Checkpoint 2 < Checkpoint 1, it is due to ODE45');
    disp('stepping backwards in time)');
    fprintf('\n');
    disp('              F-16 Attitude at Start:');
    fprintf('Roll   (phi):          %9.3f   deg \n', ...
        rad2deg(x_f16_hist(4,1)));
    fprintf('Pitch  (theta):        %9.3f   deg \n', ...
        rad2deg(x_f16_hist(5,1)));
    fprintf('Yaw    (psi):          %9.3f   deg \n', ...
        rad2deg(x_f16_hist(6,1)));
    fprintf('\n');
    disp('              F-16 Attitude at End:');
    fprintf('Roll   (phi):          %9.3f   deg \n', ...
        rad2deg(x_f16_hist(4,end)));
    fprintf('Pitch  (theta):        %9.3f   deg \n', ...
        rad2deg(x_f16_hist(5,end)));
    fprintf('Yaw    (psi):          %9.3f   deg \n', ...
        rad2deg(x_f16_hist(6,end)));
    fprintf('\n');
    disp('              Pass Fail Conditions:');
    disp(passFail);
    disp('----------------------------------------------------------');
end
----------------------------------------------------------
RESULTS:
              Min & Max Values:
Min Altitude:             96.178   ft 
Altitude Lost:          3403.822   ft 
Min Airspeed:            540.000   ft/sec 
Max Airspeed:            639.104   ft/sec 
Max Down Force:            5.222   g's 
Min Down Force:           -0.426   g's 
Max Side Force:            0.061   g's 
Max Roll Accel:          233.522   deg/sec^2 

              Maneuver Event Times:
Maneuver Start:            2.000   sec 
Checkpoint  1:             3.395   sec 
Maneuver Complete:         8.000   sec 
Final Time:               15.000   sec 
(if Checkpoint 2 < Checkpoint 1, it is due to ODE45
stepping backwards in time)

              F-16 Attitude at Start:
Roll   (phi):             45.000   deg 
Pitch  (theta):          -72.000   deg 
Yaw    (psi):            -45.000   deg 

              F-16 Attitude at End:
Roll   (phi):             -0.000   deg 
Pitch  (theta):            2.249   deg 
Yaw    (psi):            -33.514   deg 

              Pass Fail Conditions:
           stable: 1
         airspeed: 1
            alpha: 1
             beta: 1
               Nz: 1
    psMaxAccelDeg: 1
         altitude: 1
     maneuverTime: 1

----------------------------------------------------------

Create Plots

if(plotOn)
    if(printOn)
        disp('Begin plotting');
    end
    % ENU Path
    figure(10);
    hold on;
    grid on;
    fig1titlestr = sprintf('F-16 Path: t_f = %4.2f seconds',time(end));
    title(fig1titlestr);
    xlabel('East (ft)');
    ylabel('North (ft)');
    zlabel('Altitude (ft)');
    hold on;
    scatter3(x_f16_hist(11,1),x_f16_hist(10,1),x_f16_hist(12,1));
    plot3(x_f16_hist(11,:),x_f16_hist(10,:),x_f16_hist(12,:));
    legend('Start','Path')

    % Get indices of maneuver times in result time steps
    time_maneuver = zeros(size(t_maneuver));
    for (i=1:length(t_maneuver))
        if(t_maneuver(i) >= 0)
            time_maneuver(i) = find(time<=t_maneuver(i),1,'last');
        else
            time_maneuver(i) = 1; %
            warning('Maneuver Not Completed');
        end
    end

    scatter3(x_f16_hist(11,time_maneuver),...
        x_f16_hist(10,time_maneuver),...
        x_f16_hist(12,time_maneuver),'*');
    axis image
    view(3)


    % Attitude
    figure(2);
    hold on;
    grid on;
    title('Attitude History');
    xlabel('Time (sec)');
    ylabel('Attitudes & Rates (deg, deg/s)');
    hold on;
    plot(time, rad2deg(x_f16_hist(4,:)),'r');       % phi (deg)
    plot(time, rad2deg(x_f16_hist(5,:)),'g');       % theta (deg)
    plot(time, rad2deg(x_f16_hist(6,:)),'b');       % psi (deg)
    plot(time, rad2deg(x_f16_hist(7,:)),'r--');     % p (deg/s)
    plot(time, rad2deg(x_f16_hist(8,:)),'g--');     % q (deg/s)
    plot(time, rad2deg(x_f16_hist(9,:)),'b--');     % r (deg/s)
    legend('Roll \phi',' Pitch \theta','Yaw \psi',...
        'Roll Rate p','Pitch Rate q','Yaw rate r',...
        'location','SouthEast')


    % Alpha and Beta
    figure(3);
    hold on;
    grid on;
    title('Alpha & Beta');
    xlabel('Time (sec)');
    ylabel('Deflection (deg)');
    hold on;
    plot(time, rad2deg(x_f16_hist(2,:)),'r-');      % alpha (deg)
    plot(time, rad2deg(x_f16_hist(3,:)),'b-');      % beta (deg)
    plot(time, ones(size(time))*(-10),'r--');       % alpha min (deg)
    plot(time, ones(size(time))*45,'r-.');          % alpha max (deg)
    plot(time, ones(size(time))*(-30),'b--');       % beta min (deg)
    plot(time, ones(size(time))*30,'b-.');          % beta max (deg)
    legend('\alpha','\beta','\alpha_{min}','\alpha_{max}',...
        '\beta_{min}','\beta_{max}',...
        'location','SouthEast')

    % Vt
    figure(4);
    grid on;
    title('Airspeed');
    xlabel('Time (sec)');
    ylabel('KIAS (ft/sec)');
    hold on;
    plot(time, ones(size(time))*xequil(1),'k:');
    plot(time,x_f16_hist(1,:));
    plot(time, ones(size(time))*300,'r--');
    plot(time, ones(size(time))*900,'r-.');
    legend('V_{trim}','V','V_{min}','V_{max}','location','SouthEast')


    % All Control Signals
    figure(5);
    hold on;
    grid on;
    title('F-16 Control');
    xlabel('Time (sec)');
    ylabel('Control (deg & percent)');
    hold on;
    plot(time,u_hist(1,:));
    plot(time,rad2deg(u_hist(2,:)));
    plot(time,rad2deg(u_hist(3,:)));
    plot(time,rad2deg(u_hist(4,:)));
    legend('Throttle','Elevator','Aileron','Rudder','location',...
        'NorthWest')


    % Longitudinal States & Control
    figure(7)
    hold on;
    title('Longitudinal States & Control');
    xlabel('Time (sec)');
    ylabel('Nz (g''s), States (rad), States & Ctrl (rad)');
    % State Order:  alpha, q, Nz
    plot(time, (x_f16_hist(2,:)),'g');              % Alpha (rad)
    plot(time, (x_f16_hist(8,:)),'b');              % q (rad/s)
    plot(time, u_hist(5,:),'k-.');                  % Nz_r
    plot(time, Nz_hist,'r-');                       % Nz (g's)
    plot(time, lin_Nz_hist,'r--');                  % Linear Nz (g's)
    plot(time, x_f16_hist(14,:),'r:');              % Nz_i
    plot(time, (u_hist(2,:)),'c');                  % Elevator (rad)
    legend('\alpha','q','Nz_{ref}','Nz','Linear Nz','Nz_i','Elevator');


    % Lateral States & Control
    figure(8)
    hold on;
    title('Lateral States & Control');
    xlabel('Time (sec)');
    ylabel('ps (rad/sec), Ny_r (g''s), States & Ctrl (rad)');
    % State Order:  beta p r ps_i Ny_r_i
    plot(time, (x_f16_hist(3,:)),'g');              % Beta (rad)
    plot(time, u_hist(6,:),'k-.');                  % ps_ref
    plot(time, ps_hist,'r');                        % ps
    plot(time, x_f16_hist(15,:),'r:');              % ps_i
    plot(time, Ny_r_hist,'b-');                     % Ny_r (g's)
    % plot(time, lin_Ny_r_hist,'r--');              % Linear Nz_r (g's)
    plot(time, x_f16_hist(16,:),'b:');              % Ny_r_i
    plot(time, (u_hist(3,:)),'c');                  % Aileron (rad)
    plot(time, (u_hist(4,:)),'m');                  % Rudder (rad)
    legend('\beta','p_{s,ref}','p_s','p_{s,i}',...
        'Ny+r','Ny+r_i','Aileron','Rudder');

    if(printOn)
        % Record processing time
        t_plotted = toc(timerStart);

        disp('Plot Generation Complete');
        disp('Time to calculate (MM:SS.ms):');
        fprintf('%s\n',datestr((t_plotted-t_analysis)/(24*60*60),...
            'MM:SS.FFF'))
        fprintf('\n');
    end
end

if(profilerOn)
    profile viewer;
    profile off;
end
Begin plotting
Plot Generation Complete
Time to calculate (MM:SS.ms):
00:01.971

Error using RunF16Sim (line 57)
Insufficient Arguments given
end