Simulate an F16 performing a coordinated turn.

The F16 will go unstable due to an unstable longitudinal mode. ------------------------------------------------------------------------ See also ECIToNED, VTToVB, AC, ACBuild, ACInit, ACPlot, DrawAC, HUD, @acstate/acstate.m, Eul2Q, QMult, TimeGUI ------------------------------------------------------------------------

Global for the time GUI
Global for the HUD
F16 database
Load the standard atmosphere
Control
Initial state vector
Initial time and state
Initialize the model
Set up the HUD
Set up the aircraft display
Initialize the plots
Initialize the time display
Create the plots

%--------------------------------------------------------------------------
%   Copyright (c) 1997 Princeton Satellite Systems, Inc.
%   All rights reserved.
%--------------------------------------------------------------------------
%   Since version 2.0 (ACT)
%--------------------------------------------------------------------------

Global for the time GUI

%------------------------
global simulationAction
simulationAction = ' ';

Global for the HUD

%-------------------
global hUDOutput
hUDOutput = struct('pushbutton1',0,'pushbutton2',0,'checkbox1',0,...
                   'checkbox2',0,'checkbox3',0);

F16 database

%-------------
d               = ACBuild('F16');
d.theta0        = 0;
d.wPlanet       = [0;0;0];
d.actuator.name = [];
d.aero.name     = 'ACAero';
d.engine.name   = 'ACEngine';
d.rotor.name    = [];
d.sensor.name   = 'ACSensor';
d.disturb.name  = [];

Load the standard atmosphere

%-----------------------------
d.atmData       = load('AtmData.txt');
d.atmUnits      = 'eng';

Control

%--------
d.control.throttle  =   0.8349601   + 0.00177955182723;
d.control.elevator  =  -1.481766    + 0.00136225302829;
d.control.aileron   =   0.09553108  - 6.122723107626488e-04;
d.control.rudder    =  -0.4118124   + 0.00314275758041;

Initial state vector

%---------------------
alpha     = 0.23929886303038;
beta      = 5.066120644138852e-04;
vT        = 502;
v         = VTToVB( vT, alpha, beta );

cG        = [0.35;0;0];

r         = [0;0;0] + [d.equatorialRadiusPlanet;0;0];

eulInit   = [1.36428891652801;0.05049248496658;0.23546606115943];

qNEDToB   = Eul2Q( eulInit );
qECIToNED = ECIToNED( r, 'quaternion' );
q         = QMult( qECIToNED, qNEDToB );
w         = [-0.01567467524226;0.29324283670362;0.06118889106375];

wR        = 160;
engine    = 64.12363;
mass      = 1/1.57e-3;
inertia   = [9497;55814;63100;0;-982;0];
actuator  = [];
sensor    = [];
flex      = [];
disturb   = [];

Initial time and state

%-----------------------
t = 0;
x = acstate( r, q, w, v, wR, mass, inertia, cG, engine, actuator, sensor, flex, disturb );

Initialize the model

%---------------------
dT   = 0.1;
nSim = 10/dT;

d    = ACInit( x, d );

Set up the HUD

%---------------
dHUD.atmData  = d.atmData ;
dHUD.atmUnits = 'eng';

cHUD.control     = d.control;
cHUD.elevatorMax = 90;
cHUD.aileronMax  = 90;
cHUD.rudderMax   = 90;
cHUD.dT          = dT;
hHUD = HUD( 'init', dHUD, x, [], cHUD );

Set up the aircraft display

%----------------------------
gF16 = load('gF16');
hF16 = DrawAC( 'init', gF16, x, [], d.atmUnits  );

Initialize the plots

%---------------------
plots = [ 'Euler angles';...
          'Quaternion  ';...
          'Angular rate';...
          'Position ECI';...
          'Velocity    ';...
          'Alpha       '];

dPlot = ACPlot( x, 'init', plots, d, nSim, dT, nSim );

Initialize the time display

%----------------------------
tToGoMem.lastJD        = 0;
tToGoMem.lastStepsDone = 0;
tToGoMem.kAve          = 0;
[ ratioRealTime, tToGoMem ] =  TimeGUI( nSim, 0, tToGoMem, 0, dT, 'F16 Simulation' );

for k = 1:nSim

	% Display the status message
	%---------------------------
	[ ratioRealTime, tToGoMem ] = TimeGUI( nSim, k, tToGoMem, ratioRealTime, dT );

  % HUD information
  %----------------
  hHUD = HUD( 'run', dHUD, x, hHUD, cHUD );

  % Plotting
  %---------
  dPlot = ACPlot( x, 'store', dPlot, k );

  % 3D Display
  %-----------
  hF16 = DrawAC( 'run', gF16, x, hF16, d.atmUnits );

  % The simulation
  %---------------
	x     = AC( x, t, dT, d );
	t     = t + dT;

  % Time control
  %-------------
  switch simulationAction
    case 'pause'
      pause
      simulationAction = ' ';
    case 'stop'
      return;
    case 'plot'
      break;
  end

end

TimeGUI('close');

Create the plots

%-----------------
ACPlot( x, 'plot', dPlot );



%--------------------------------------