% RESULTS - processes aircraft simulation results 
%
% The Matlab macro RESULTS creates time-trajectories of individual inputs 
% to and outputs from the Beaver model. It can be used after running a 
% simulation involving the aircraft model.
%
% During simulations, the time-trajectories of all inputvariables are sent 
% to the matrix In and all time-trajectories of the outputvariables are 
% sent to to the matrix Out in the Matlab workspace; see the input and 
% output definitions for the Beaver model (inputs.htm and outputs.htm in 
% the HELP subdirectory) for the definitions of these matrices. Since 
% analyzing the results directly from In and Out is possible, but not very 
% convenient, RESULTS was created to transform these matrices into time-
% trajectories of individual input and output signals with self-explaining 
% variable names.
%
% So in order to produce a plot of the elevator deflection (for instance), 
% instead of typing plot(time,In(:,1)), you can run RESULTS and then type 
% plot(time,deltae). Even easier: you can also run RESULTS first and then 
% use the Matlab macro RESPLOT to automatically plot the most important 
% time-trajectories.
%
% Note: RESULTS only works correctly if the definitions of the matrices In 
% and Out are not changed by editing the aircraft model. If the definitions 
% are changed, then the file RESULTS.M in in the subdirectory TOOLS needs 
% to be edited accordingly. Apart from some options regarding the units of 
% measurements of the variables that represent angles and/or angular 
% velocities, the macro RESULTS is very straightforward, so adapting this 
% file is not very difficult.
%
% RESULTS does not contain a 'save' option, so remember to save the 
% resulting time-trajectories yourself if you want to use them later. 
% Don't forget to store the time-axis time too; otherwise you won't be able
% to plot the time-trajectories of the variables afterwards.
%
% Type 'browse inputs' or 'browse outputs' (without quotes) for more info.


% Check if the matrices In and Out, and the time axis are present in the
% Matlab workspace. If not, display error message. 
% -----------------------------------------------------------------------
if exist('time') == 0 | exist('Out') == 0 | exist('In') == 0
   disp('Aircraft simulation results (matrices ''In'' and ''Out'' and');
   disp('time vector ''time'' are not present in the workspace!');
   error;
end

% Check if the number of inputs and outputs stored in these matrices cor-
% respond with the default definitions of In and Out used in the system
% BEAVER. Of course, the latter check does not guarantee that the defini-
% tions of the matrices In and Out as found in the Matlab workspace, are in
% accordance with the default definitions used in the system BEAVER, but if
% the number of variables does not match we know for sure that the results
% aren't correct!
%
% The user him/herself must take care of changing RESULTS.M if the defini-
% tions of the matrices In and Out need to be changed!
% ----------------------------------------------------------------------------

[m1,n1] = size(In);
[m2,n2] = size(Out);

if n1 ~= 12
   error('Number of inputs is not equal to 12. RESULTS won''t work!');
elseif n2 ~= 89
   error('Number of outputs is not equal to 89. RESULTS won''t work!');
end

clear m1 n1 m2 n2

% Header to make sure user knows what goes on...
% -----------------------------------------------
clc
disp('FDC toolbox - RESULTS');
disp('=====================');
disp(' ');
disp('Transfer simulation results from matrices In and Out into');
disp('time-trajectories of individual inputs and outputs, using');
disp('self-explaining variable names.');
disp(' ');

% Specify dimension of angles and angular velocities. While-loops ensure
% that that user gives valid answers to the questions. Default is [rad]
% and [rad/s].
% ----------------------------------------------------------------------
ok = 0;
while ok ~= 1
   disp(' ');
   angledims = input('Angles in [rad] or [deg] ([r]/d)? ','s');

   if isempty(angledims)
      angledims = 'r';
   end
   
   if angledims == 'r'
      disp(' ');
      disp('    Angles will be measured in [rad]');
      c1 = 1;                      % multiplication factor for angles in [rad]
      ok = 1;
      unit_angle = '[rad]';        % angle unit for display in RESPLOT
   else
      disp(' ');
      disp('    Angles will be measured in [deg]');
      c1 = 180/pi;                 % multiplication factor for angles in [deg]
      ok = 1;
      unit_angle = '[deg]';        % angle unit for display in RESPLOT
   end
end

ok = 0;
while ok ~= 1
   disp(' ');
   ratedims = input('Angular velocities in [rad/s] or [deg/s] ([r]/d)? ','s');

   if isempty(ratedims)
      ratedims = 'r';
   end

   if ratedims == 'r'
      disp(' ');
      disp('    Angular velocities will be measured in [rad/s]');
      c2 = 1;                 % multiplication factor for angular velocities in [rad/s]
      ok = 1;
      unit_rate = '[rad/s]';  % angular rate unit for display in RESPLOT
   else
      disp(' ');
      disp('    Angular velocities will be measured in [deg/s]');
      c2 = 180/pi;   % multiplication factor for angular velocities in [deg/s]
      ok = 1;
      unit_rate = '[deg/s]';  % angular rate unit for display in RESPLOT
   end
end

disp(' ');
disp(' ');
disp('Wait a moment please...');

% Now create individual trajectories of control inputs and atmospheric
% disturbances from the results, stored in the matrix In. See INPUTS.HLP
% for the exact definitions.
% ----------------------------------------------------------------------
deltae = In(:,1) * c1;
deltaa = In(:,2) * c1;
deltar = In(:,3) * c1;
deltaf = In(:,4) * c1;
n      = In(:,5);
pz     = In(:,6);
uw     = In(:,7);
vw     = In(:,8);
ww     = In(:,9);
uwdot  = In(:,10);
vwdot  = In(:,11);
wwdot  = In(:,12);

% Now create individual trajectories of outputs from the system BEAVER,
% from the results stored in the matrix Out. See OUTPUTS.HLP for the
% exact definitions.
% ---------------------------------------------------------------------

% State variables:

V        = Out(:,1);
alpha    = Out(:,2)* c1;
beta     = Out(:,3)* c1;
p        = Out(:,4)* c2;
q        = Out(:,5)* c2;
r        = Out(:,6)* c2;
psi      = Out(:,7)* c2;
theta    = Out(:,8)* c2;
phi      = Out(:,9)* c2;
xe       = Out(:,10);
ye       = Out(:,11);
H        = Out(:,12);

% Time-derivatives of states (see VABDOT.HLP, PQRDOT.HLP, EULERDOT.HLP, and
% XYHDOT.HLP):

Vdot     = Out(:,13);
alphadot = Out(:,14)* c1;
betadot  = Out(:,15)* c1;
pdot     = Out(:,16)* c2;
qdot     = Out(:,17)* c2;
rdot     = Out(:,18)* c2;
psidot   = Out(:,19)* c2;
thetadot = Out(:,20)* c2;
phidot   = Out(:,21)* c2;
xedot    = Out(:,22);
yedot    = Out(:,23);
Hdot     = Out(:,24);

% Body-axes velocities (see UVW.HLP):

u        = Out(:,25);
v        = Out(:,26);
w        = Out(:,27);

% Time-derivatives of body-axes velocities (see UVWDOT.HLP):

udot     = Out(:,28);
vdot     = Out(:,29);
wdot     = Out(:,30);

% Dimensionless rotational velocities in aircraft's body-axes (see
% DIMLESS.HLP):

pb_2V    = Out(:,31);
qc_V     = Out(:,32);
rb_2V    = Out(:,33);

% Some additional flightpath variables (see FLPATH.HLP):

gamma    = Out(:,34)* c1;
fpa      = Out(:,35);
chi      = Out(:,36)* c1;
Phi      = Out(:,37)* c1;

% Engine variables (see POWER.HLP):

dpt      = Out(:,38);
P        = Out(:,39);

% Specific forces and kinematic accelerations (see ACCEL.HLP):

Ax       = Out(:,40);
Ay       = Out(:,41);
Az       = Out(:,42);
axk      = Out(:,43);
ayk      = Out(:,44);
azk      = Out(:,45);

% Dimensionless aerodynamic force and moment coefficients (see AEROMOD.HLP):

CXa      = Out(:,46);
CYa      = Out(:,47);
CZa      = Out(:,48);
Cla      = Out(:,49);
Cma      = Out(:,50);
Cna      = Out(:,51);

% Dimensionless propulsive force and moment coefficients (see ENGMOD.HLP):

CXp      = Out(:,52);
CYp      = Out(:,53);
CZp      = Out(:,54);
Clp      = Out(:,55);
Cmp      = Out(:,56);
Cnp      = Out(:,57);

% Aerodynamic forces and moments (see AEROMOD.HLP and FMDIMS.HLP):

Xa       = Out(:,58);
Ya       = Out(:,59);
Za       = Out(:,60);
La       = Out(:,61);
Ma       = Out(:,62);
Na       = Out(:,63);

% Propulsive forces and moments (see ENGMOD.HLP and FMDIMS.HLP):

Xp       = Out(:,64);
Yp       = Out(:,65);
Zp       = Out(:,66);
Lp       = Out(:,67);
Mp       = Out(:,68);
Np       = Out(:,69);

% Forces due to gravity (see GRAVITY.HLP):

Xgr      = Out(:,70);
Ygr      = Out(:,71);
Zgr      = Out(:,72);

% Corrections to body-axes forces due to nonsteady atmosphere (see
% FWIND.HLP):

Xw       = Out(:,73);
Yw       = Out(:,74);
Zw       = Out(:,75);

% Atmosphere variables (see ATMOSPH.HLP):

rho      = Out(:,76);
ps       = Out(:,77);
T        = Out(:,78);
mu       = Out(:,79);
g        = Out(:,80);

% Airdata variables (see AIRDATA1.HLP, AIRDATA2.HLP, and AIRDATA3.HLP):

a        = Out(:,81);
M        = Out(:,82);
qdyn     = Out(:,83);

qc       = Out(:,84);
Ve       = Out(:,85);
Vc       = Out(:,86);

Tt       = Out(:,87);
Re       = Out(:,88);
Rc       = Out(:,89);

% Ready...
% --------
disp(' ');
disp('Ready. Press key to continue...');
pause
clc

% Give user the option to delete matrices In and Out from the workspace
% ---------------------------------------------------------------------
answ = input('Delete matrices In and Out from workspace (y/[n])? ','s');

if isempty(answ)
   answ = 'n';
end

if answ == 'y'
   clear In Out
end

% Display help message for user. Goodbye, see U next time...
% ----------------------------------------------------------
disp(' ');
disp('The workspace now contains time trajectories of the twelve inputs');
disp('(six vectors with control inputs and six with atmospheric distur-');
disp('bances) and the 89 outputs from the system BEAVER. The variable');
disp('names of these vectors are self-explaining, but if you have doubts');
disp('about the definitions, examine the help texts INPUTS.HLP and/or');
disp('OUTPUTS.HLP from the subdirectory HELP.');
disp(' ');
disp('The time-axis against which the results can be plotted is con-');
disp('tained in the variable time.');
disp(' ');
disp('Remember that currently:');
disp(' ');
if angledims == 'd'
   disp('  angles are measured in [deg]');
else
   disp('  angles are measured in [rad]');
end
disp(' ');
if ratedims == 'd'
   disp('  angular velocities are measured in [deg/s]');
else
   disp('  angular velocities are measured in [rad/s]');
end
disp(' ');
disp('Ready.');
disp(' ');


% Clear interim variables. Note: unit_angle and unit_rate are 
% not cleared, because they may be used by RESPLOT!
% -----------------------------------------------------------
clear ok answ angledims ratedims c1 c2


%-----------------------------------------------------------------------
% The Flight Dynamics and Control Toolbox version 1.4.0. 
% (c) Copyright Marc Rauw, 1994-2005. Licensed under the Open Software 
% License version 2.1; see COPYING.TXT and LICENSE.TXT for more details.
% Last revision of this file: July 17, 2005. 
%-----------------------------------------------------------------------