% MODBUILD - Builds a datafile with Aircraft model parameters
%

% The Matlab script MODBUILD is used to build a datafile with
% parameters for the aircraft model. Although the current version 
% of MODBUILD is only valid for the DHC-2 'Beaver' model (the 
% Simulink system BEAVER or one of its subsystem equivalents),
% the part in which the inertia coefficienst are computed is more 
% generic in nature (see the comment lines in the MODBUILD.M 
% source code for more details).
%
% MODBUILD needs to be used only if the aircraft model parameters 
% are to be changed. This includes changes in mass and inertia 
% properties, since these are considered to be constant during the 
% short time-intervals over which the aircraft responses are 
% considered. MODBUILD lets you save the parameters to a datafile
% (by default: AIRCRAFT.DAT), which can be loaded for future 
% sessions, using the DATLOAD utility. It is not necessary to run 
% MODBUILD more than once for the same set of model parameters if 
% the results are saved.
%
% Changing the model parameters with the current version of 
% MODBUILD can be done only by editing the source-file of the 
% program! In practice, MODBUILD won't have to be applied very 
% often, so this is not really as inconvenient as it may seem. 
% Moreover, it is straightforward to add some kind of user-
% interface if required. Currently, MODBUILD functions as a 
% simple batch-file. The user is informed about the construction 
% of the datafile by means of echo on and echo off commands.
%
% Change the aircraft-dependent parts of the program if you want 
% to implement model-parameters for different types of aircraft. 
% It is also possible to change the mass, geometry, and moments 
% and products of inertia, but be careful not to change the 
% generic inertia coefficient equations.
%
% See also DATLOAD.


% First, we'll save the data which is currently available in the workspace
% to a temporaty file. This is to make sure that no important data will be
% deleted during execution of MODBUILD.
% ------------------------------------------------------------------------
save modbuild.tmp


% MODBUILD HEADER.
%=================
clc;
disp('The FDC toolbox - MODBUILD');
disp('==========================');
disp(' ');
disp('Build datafiles with parameters for non-linear aircraft model.');
disp(' ');
disp('Results are valid for the DHC-2 ''Beaver'' aircraft. Change the');
disp('program MODBUILD.M if models of other aircraft are to be used.');
disp(' ');
disp(' ');
disp('<< Press key to continue >>');
pause
clc

% DEFINE PARAMETERS FOR THE AERODYNAMIC MODEL
% ===========================================

% CURRENT CONFIGURATION: Standard nonlinear model of the DHC-2 'Beaver',
% according to [Tjee and Mulder, 1988]. The aerodynamic force and moment
% coefficients for this aircraft can be written as nonlinear polynomial
% functions of the state and input variables. This model description is
% very compact, because only one set of stability and control derivatives
% (= polynomial coefficients) is needed.
%
% For other aircraft, it is more customary to determine the aerodynamic
% forces and moments by interpolating in large multi-dimensional tables.
% Although such models use significantly more parameter matrices, it is
% really quite straightforward to replace the following matrix definitions
% by similar definitions for other aircraft.
%
% Here, the aerodynamic stability and control derivatives are stored in
% the matrix AM. For aerodynamic models, consisting of multiple tables,
% it is recommended to use variable names like AM1, AM2, etc. to maintain
% commonality with the current 'Beaver' model BEAVER.

echo on

  % Define stability and control derivatives of the DHC-2 'Beaver'
  % for the nonlinear aerodynamic model, which is valid within the
  % 35-55 m/s TAS range (see [Tjee & Mulder, 1988]).
  % --------------------------------------------------------------
 
  CX0   = -0.01804;    CZ0    = -0.1198;     Cm0  = -0.2036;      %1int  2int  3int
  CXa   =  0.1118;    CZa    = -5.731;      Cma  = -5.831;       %1int  2int  3int
  CXa2  =  5.03296;     CZa3   =  9.719;      Cma2 =  1.288;       %1int  2int  3int
  CXa3  =  0;          CZq    =  16.94;      Cmq  = -31.35;       %1int  2dtcm 3dtcm
  CXq   =  0;          CZde   = -0.2634;     Cmde = -0.9516;      %1el   2int  3dtcm
  CXdr  =  0;          CZdeb2 =  0.0;        Cmb2 =  0.0;         %1el   2el   3el
  CXdf  = -0.01334;    CZdf   = -0.7792;     Cmr  = -0.0;         %1int  2int  3el
  CXadf =  0.4095;     CZadf  = -0.4301;     Cmdf = -0.8073;      %1calcolostrano 2calcolstrano 3dtcm

  CY0   = -0.0;        Cl0    =  0.0;        Cn0  = -0.0;         %1el   2el   3el  
  CYb   = -0.4271;     Clb    = -0.03905;    Cnb  =  0.1316;      %1dtcm 2dtcm 3dtcm
  CYp   = -0.0739;     Clp    = -0.41547;    Cnp  = -0.055921;    %1dtcm 2dtcm 3dtcm
  CYr   =  0.3076;     Clr    =  0.09279;    Cnr  = -0.37982;     %1dn   2dtcm 3dtcm
  CYda  = -0.0;        Clda   = -0.1467;     Cnda =  0.00161;     %1dn   2int  3int
  CYdr  =  0.0534;     Cldr   =  0.0033;     Cndr = -0.0349;      %1dn   2dn   3dn    %perplessitÃ  sul CLaphav nel calcolo di CYdr
  CYdra =  0.0;        Cldaa  = -0.0;        Cnq  =  0.0;         %1el   2el   3el
  CYbdot=  0.0;                              Cnb3 =  0.0;         %1dn         3el

  
echo off

AM = [ CX0     CY0     CZ0     Cl0     Cm0     Cn0  ;
       CXa     0       CZa     0       Cma     0    ;
       CXa2    0       0       0       Cma2    0    ;
       CXa3    0       CZa3    0       0       0    ;
       0       CYb     0       Clb     0       Cnb  ;
       0       0       0       0       Cmb2    0    ;
       0       0       0       0       0       Cnb3 ;
       0       CYp     0       Clp     0       Cnp  ;
       CXq     0       CZq     0       Cmq     Cnq  ;
       0       CYr     0       Clr     Cmr     Cnr  ;
       0       0       CZde    0       Cmde    0    ;
       CXdf    0       CZdf    0       Cmdf    0    ;
       0       CYda    0       Clda    0       Cnda ;
       CXdr    CYdr    0       Cldr    0       Cndr ;
       CXadf   0       CZadf   0       0       0    ;
       0       CYdra   0       0       0       0    ;
       0       0       0       Cldaa   0       0    ;
       0       0       CZdeb2  0       0       0    ;
       0       CYbdot  0       0       0       0    ];

AM = AM';


disp(' ');
disp('<< Press key to continue >>');
pause
clc


% DEFINE PARAMETERS FOR ENGINE FORCES & MOMENTS MODEL.
% ====================================================

% CURRENT CONFIGURATION: Standard nonlinear model of the DHC-2 'Beaver',
% according to [Tjee and Mulder, 1988].
%
% The engine forces and moments model of the DHC-2 'Beaver' also expresses
% the force and moment coefficients as nonlinear polynomial functions of
% external inputs and state variables.
%
% Here, the engine stability and control derivatives are stored in the
% matrix EM. For engine models, consisting of multiple tables, it is re-
% commended to use variable names like EM1, EM2, etc. to maintain com-
% monality with the current 'Beaver' model BEAVER.

echo on

  % The nonlinear engine model of the DHC-2 "BEAVER" aircraft
  % valid within the 35-55 m/sec TAS-range (see [Tjee & Mulder, 1988]).
  % -------------------------------------------------------------------
  CXdpt   =  0.0;
  CXadpt2 =  0.0;
  CZdpt   = -0.0;
  Cla2dpt = -0.0;
  Cmdpt   = -0.0;
  Cndpt3  = -0.0;

echo off

EM = [ CXdpt    0        CZdpt    0        Cmdpt    0      ;
       0        0        0        0        0        Cndpt3 ;
       CXadpt2  0        0        0        0        0      ;
       0        0        0        Cla2dpt  0        0      ];

EM = EM';

disp(' ');
disp('<< Press key to continue >>');
pause
clc


% DEFINITION OF RELEVANT AIRCRAFT GEOMETRY AND MASS-DISTRIBUTION PARAMETERS.
%
% (Note: geometry and mass properties are assumed to be constant during the
% motions of interest! If it is required to compute the influence of changes
% in these properties on-line during the simulation, many blocks in the sys-
% tems BEAVER or BEAVER1 have to be changed. First of all, a block which
% computes the geometrical and/or mass properties must be added to the sys-
% tem. Then use a text editor to find the blocks in which the parameter ma-
% trices GM1 and GM2 are used and replace these parameters by variables.
% Add corresponding Inport blocks to these blocks, and make the proper
% connections to the block in which the geometry and mass properties are
% calculated.)
% ==========================================================================

echo on

  % Aircraft data on which the aerodynamic model is based. CHANGE THE
  % VARIABLES ACCORDING TO YOUR OWN WISHES BY EDITING MODBUILD.
  % -----------------------------------------------------------------

  % Mass and mass-distribution.
  % ---------------------------
  Ix      =  Ixx_tot;               % kgm^2 in Fr
  Iy      =  Iyy_tot;               % kgm^2 in Fr
  Iz      =  Izz_tot;               % kgm^2 in Fr
  Jxy     =     0.0;                % kgm^2 in Fr
  Jxz     =     Jxz;                % kgm^2 in Fr
  Jyz     =     0.0;                % kgm^2 in Fr
  m       =    MTOW;                % kg

  % geometric data.
  % ---------------
  cbar    =     1.252;              % m
  b       =     7.76;               % m quella 
  S       =     9.51;               % m^2

echo off

% THE FOLLOWING EQUATIONS ARE VALID FOR ANY RIGID AIRCRAFT
% (ALSO VALID FOR NON-SYMMETRIC AIRCRAFT).
%
% Calculate inertia parameters (see NASA TP 2768). The formula's
% are valid for symmetric and asymmetric aircraft.
% --------------------------------------------------------------
detI = Ix*Iy*Iz - 2*Jxy*Jxz*Jyz - Ix*Jyz^2 - Iy*Jxz^2 - Iz*Jxy^2;
I1   = Iy*Iz - Jyz^2;
I2   = Jxy*Iz + Jyz*Jxz;
I3   = Jxy*Jyz + Iy*Jxz;
I4   = Ix*Iz - Jxz^2;
I5   = Ix*Jyz + Jxy*Jxz;
I6   = Ix*Iy - Jxy^2;

Pl  = I1/detI; Pm = I2/detI; Pn = I3/detI;
Ppp = -(Jxz*I2 - Jxy*I3)/detI;
Ppq = (Jxz*I1 - Jyz*I2 - (Iy-Ix)*I3)/detI;
Ppr = -(Jxy*I1 + (Ix-Iz)*I2 - Jyz*I3)/detI;
Pqq = (Jyz*I1 - Jxy*I3)/detI;
Pqr = -((Iz-Iy)*I1 - Jxy*I2 + Jxz*I3)/detI;
Prr = -(Jyz*I1 - Jxz*I2)/detI;

Ql  = I2/detI; Qm = I4/detI; Qn = I5/detI;
Qpp = -(Jxz*I4 - Jxy*I5)/detI;
Qpq = (Jxz*I2 - Jyz*I4 - (Iy-Ix)*I5)/detI;
Qpr = -(Jxy*I2 + (Ix-Iz)*I4 - Jyz*I5)/detI;
Qqq = (Jyz*I2 - Jxy*I5)/detI;
Qqr = -((Iz-Iy)*I2 - Jxy*I4 + Jxz*I5)/detI;
Qrr = -(Jyz*I2 - Jxz*I4)/detI;

Rl  = I3/detI; Rm = I5/detI; Rn = I6/detI;
Rpp = -(Jxz*I5 - Jxy*I6)/detI;
Rpq = (Jxz*I3 - Jyz*I5 - (Iy-Ix)*I6)/detI;
Rpr = -(Jxy*I3 + (Ix-Iz)*I5 - Jyz*I6)/detI;
Rqq = (Jyz*I3 - Jxy*I6)/detI;
Rqr = -((Iz-Iy)*I3 - Jxy*I5 + Jxz*I6)/detI;
Rrr = -(Jyz*I3 - Jxz*I5)/detI;

% Summarizing results in aircraft parameter matrices GM1 and GM2.
% ---------------------------------------------------------------
GM1 = [cbar b    S    Ix   Iy   Iz   Jxy  Jxz  Jyz  m];

GM2 = [ Pl   Pm   Pn   Ppp  Ppq  Ppr  Pqq  Pqr  Prr ;
        Ql   Qm   Qn   Qpp  Qpq  Qpr  Qqq  Qqr  Qrr ;
        Rl   Rm   Rn   Rpp  Rpq  Rpr  Rqq  Rqr  Rrr ];

disp(' ');
disp('<< Press key to continue >>');
pause
clc


% Save the aircraft parameters to file. The default folder for this
% action is the subfolder DATA of the FDC toolbox, which is put in
% the variable defdir. The current directory is stored in the
% variable currentdir, which is used if the default directory can't
% be found.
% =================================================================
defdir = datadir;
currentdir = chdir;

% Go to default directory if that directory exists (if not, start
% save-routine from the current directory).
% ---------------------------------------------------------------
try
  chdir(defdir);
catch
  chdir(currentdir);
end

% Obtain path (use default filename AIRCRAFT.DAT).
% ------------------------------------------------
[filename,dirname] = uiputfile('aircraft.dat','Save aircraft model parameters');

% Save results to the specified file
% ----------------------------------
save([dirname,filename],'AM','EM','GM1','GM2');

% CONCLUDING REMARKS.
%====================
clc
disp('Ready. The data can now be loaded into the Matlab workspace');
disp('by typing: LOAD AIRCRAFT.DAT -MAT');
disp(' ');
disp('Aerodynamic data is contained in the matrix AM, engine data is');
disp('contained in EM, aircraft geometry, mass, and mass distribution');
disp('data is contained in GM1 and GM2.');
disp(' ');
disp('You may also use the macro LOADER.M to load the data matrices');
disp('in the Matlab workspace. MODBUILD doesn''t have to be used again,');
disp('unless you want to make changes to the model parameters. In that');
disp('case, MODBUILD.M needs to be edited.');
disp(' ');

% Return to previous directory, clear workspace, retrieve variables from 
% temporary MAT-file and delete temporary MAT-file. Workspace will now 
% contain the same variables as before running MODBUILD.
% ----------------------------------------------------------------------
chdir(currentdir);

clear
load modbuild.tmp -mat
delete('modbuild.tmp');


% ------------------------------------------------------------
% References:
%
% R.T.H. Tjee and J.A. Mulder. Stability and Control Deriva-
% tives of the De Havilland DHC-2 "Beaver" aircraft. Report
% LR-556, Delft University of Technology, 1988.
%
% Duke, E.L., Patterson, B.P. and Antoniewicz, R.F. User's
% manual for LINEAR, a Fortran Program to derive Linear Air-
% craft Models. NASA TP 2768, 1987.
%
% Ruijgrok: Elements of Airplane Performance, Delft University
% Press, DUT (L&R), 1990.
%
% Rauw, M.O.: A Simulink environment for Aircraft Dynamics and
% Control analysis - Application to the DHC-2 'Beaver',
% MSc thesis, Delft University of Technology (L&R), 1993.
%
% Rauw, M.O.: FDC 1.2 - A Simulink Toolbox for Flight Dynamics
% and Control Analysis, 1997.
% ------------------------------------------------------------


%-----------------------------------------------------------------------
% 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: May 10, 2005.
%-----------------------------------------------------------------------