% Initialization commands
% -----------------------
format short e;
options = [];
turntype = 'c';

% Display header; welcome to ACTRIM!
% ----------------------------------
clc
disp('The FDC toolbox - ACTRIMFOR');
disp('========================');
disp(' ');
disp('This program searches determines a steady-state trimmed-flight condition');
disp('for a non-linear aircraft model in Simulink, for defined V and H, changing pz.');
disp(' ');
disp(' ');

% Parameters for aircraft model:
%
% AM, EM, GM1, and GM2: matrices, containing parameters for the aircraft
%                       model. See also the source-code of MODBUILD.M
%                       (contained in the directory AIRCRAFT)
%
% The other model parameters will either be explained when they are used,
% or are self-explaining (at least if you know something about aircraft
% stability and control).
%
% Check if AM, EM, GM1, and GM2 have been defined in the Matlab workspace.
% If not, run DATLOAD to load them from file.
% -------------------------------------------------------------------------
if exist('AM')==0 | exist('EM')==0 | exist('GM1')==0 | exist('GM2')==0
    h=newMsgBox(['First, the model parameters need to be retrieved from file ', ...
                 '(e.g. AIRCRAFT.DAT). Click ''OK'' or press Enter to continue.']);
    uiwait(h);
    datload('aircraft');
end

% If model parameters are still not present in the workspace,
% e.g. due to an incorrect datafile, force an abort of ACTRIM.
% ------------------------------------------------------------
if exist('AM')==0 | exist('EM')==0 | exist('GM1')==0 | exist('GM2')==0
    error(['ERROR: the model parameters are still not present in the workspace! ', ...
                 'Aborting ACTRIM.']);
end


% Set xinco (= initial value of the state vector). Note: the value of xinco
% which is defined here doesn't really matter, because this variable won't
% be used during the trim process. The variable xinco must be defined in the
% Matlab Workspace, however, because Simulink needs it in order to define
% the system internally. The variable xinco is called by the Integrator
% block in the heart of the aircraft model (level 3: 'Aircraft Equations of
% Motion').
%
% In the aircraft model, some variables are made dimensionless by dividing
% by the dynamic pressure, 1/2*rho*V^2. In order to avoid 'Division By Zero'
% errors, the first element of xinco (=V) has been equalled to 45 [m/s],
% which is the mean velocity of the region in which the aerodynamic model
% of the 'Beaver' is valid. It could have been set to any value unequal to
% zero, however.
% -------------------------------------------------------------------------
xinco = [45 0 0 0 0 0 0 0 0 0 0 0]';

% The system BEAVER uses a gain block for arbitrarily setting time-
% derivatives of the state variables to zero. This may be useful for
% some purposes, for instance to eliminate longitudinal-lateral cross-
% coupling, or to fix the airspeed to its initial value if an altitude
% controller is to be evaluated, but no autothrottle is available yet.
% The gain value is xfix, which is set to 1 if all states may vary, i.e.,
% if the complete six degree-of-freedom model is used without restrictions.
% When states have to be fixed, the routine FIXSTATE.M can be called
% (type HELP FIXSTATE for more info).
%
% Here, all state variables are allowed to vary, so xfix will be set to
% one. Replace this line by:
%
%   fixstate;
%
% if you want to set some time-derivatives to zero, even when trimming
% the aircraft model.
% -------------------------------------------------------------------------
xfix = 1;
%%%% fixstate;

% The name of the aircraft model to be evaluated will be stored in the
% stringvariable sysname. The aircraft model must use the same definitions
% of the state and input vectors as the system BEAVER in order to function
% properly. Edit this program if your model uses other definitions!
% ------------------------------------------------------------------------
disp('Give name of system with aircraft model');
disp('default = beaver');
sysname = input('> ','s');
if isempty(sysname)
   sysname = 'beaver';
end

% Display menu in which the user can choose between a number of trimmed-
% flight conditions and quitting.
% ----------------------------------------------------------------------
clc
opt = txtmenu('Select type of steady-state flight',...
              'Steady wings-level flight','Steady turning flight',...
              'Steady pull-up','Steady roll','Quit');

skip  = 0; % Do not skip iteration block, unless option 'quit' is used
           % (then skip will be set to 1).

% Define flight condition, depending upon trim option chosen
% ----------------------------------------------------------
if opt == 1                                       % STEADY WINGS-LEVEL FLIGHT
                                                  % -------------------------
   clc
   disp('Steady wings-level flight.');
   disp('==========================');

   V = input('Give desired airspeed [m/s], default = 45: ');
   if isempty(V)
      V = 45;
   end

   H = input('Give (initial) altitude [m], default = 0: ');
   if isempty(H)
      H = 0;
   end

   psi  = input('Give heading [deg], default = 0: ')*pi/180;
   if isempty(psi)
      psi = 0;
   end

   ok = 0;
   while ok~= 1
      gammatype = input('Use specified manifold pressure or flight-path angle ([m]/f)? ','s');
      if isempty(gammatype)
         gammatype = 'm';
      end
      if gammatype == 'f'
         gamma = input('Give flightpath angle [deg], default = 0: ')*pi/180;
         if isempty(gamma)
            gamma = 0;
            ok = 1;
         else
             ok=1;
         end
      elseif gammatype == 'm'
         ok = 1;
      else
         disp('   Invalid entry. Please choose either ''m'', or ''f''');
      end
   end
      

   phidot   = 0;
   psidot   = 0;
   thetadot = 0;

   rolltype = 'b';   % No rolling, so for reasons of simplicity the default setting, 
                     % i.e. a body-axes roll, will be used.

elseif opt == 2                                       % STEADY TURNING FLIGHT
                                                      % ---------------------
   clc
   disp('Steady turning flight.');
   disp('======================');

   ok = 0;
   while ok~= 1
      turntype = input('Do you want a coordinated or uncoordinated turn ([c]/u)? ','s');
      if isempty(turntype)
         turntype = 'c';
      end
      if turntype == 'c' | turntype == 'u'
         ok = 1;
      else
         disp('   Invalid entry. Please choose either ''c'', or ''u''');
      end
   end

   V = input('Give desired airspeed [m/s], default = 45: ');
   if isempty(V)
      V = 45;
   end

   H = input('Give initial altitude [m], default = 0: ');
   if isempty(H)
      H = 0;
   end

   psi = input('Give initial heading [deg], default = 0: ')*pi/180;
   if isempty(psi)
      psi = 0;
   end

   ok = 0;
   while ok~= 1
      gammatype = input('Use specified manifold pressure or flight-path angle ([m]/f)? ','s');
      if isempty(gammatype)
         gammatype = 'm';
      end
      if gammatype == 'f'
         gamma = input('Give flightpath angle [deg], default = 0: ')*pi/180;
         if isempty(gamma)
            gamma = 0;
            ok = 1;
         end
      elseif gammatype == 'm'
         ok = 1;
      else
         disp('   Invalid entry. Please choose either ''m'', or ''f''');
      end
   end

   phidot = 0;
   thetadot = 0;

   ok = 0;
   while ok~= 1
      answ = input('Use specified turnrate or radius ([t]/r)? ','s'); 
      if isempty(answ)
         answ = 't';
      end
      if answ == 't'
         psidot = input('Give desired rate of turn [deg/s], default = 0: ')*pi/180;
         if isempty(psidot)
            psidot = 0;
         end
         ok = 1;
      elseif answ == 'r'
         R = input('Give desired turn radius (>0) [m], default = straight & level: ');
         if isempty(R) ~= 1
            psidot = V/R;
         else
            psidot = 0;
         end
         ok = 1;
      else
         disp('   Invalid entry. Please choose either ''t'', or ''r''');
      end
   end

   rolltype = 'b';   % No rolling, so for reasons of simplicity the default setting, 
                     % i.e. a body-axes roll, will be used.

elseif opt == 3                                              % STEADY PULL-UP
                                                             % --------------
   clc
   disp('Steady pull-up.');
   disp('===============');

   V = input('Give desired airspeed [m/s], default = 45: ');
   if isempty(V)
      V = 45;
   end

   H = input('Give initial altitude [m], default = 0: ');
   if isempty(H)
      H = 0;
   end

   psi = input('Give initial heading [deg], default = 0: ')*pi/180;
   if isempty(psi)
      psi = 0;
   end

   gammatype = 'f';  % Use specified flightpath angle and numerically
                     % adjust manifold pressure
   gamma = input('Give initial flightpath angle [deg], default = 0: ')*pi/180;
   if isempty(gamma)
      gamma = 0;
   end

   phidot = 0;
   psidot = 0;

   thetadot = input('Give pull-up rate [deg/s], default = 0: ')*pi/180;
   if isempty(thetadot)
      thetadot = 0;
   end

   rolltype = 'b';   % No rolling, so for reasons of simplicity the default setting, 
                     % i.e. a body-axes roll, will be used.

elseif opt == 4                                                 % STEADY ROLL
                                                                % -----------
   clc
   disp('Steady roll.');
   disp('============');

   V = input('Give desired airspeed [m/s], default = 45: ');
   if isempty(V)
      V = 45;
   end

   H = input('Give initial altitude [m], default = 0: ');
   if isempty(H)
      H = 0;
   end

   psi = input('Give initial heading [deg], default = 0: ')*pi/180;
   if isempty(psi)
      psi = 0;
   end

   ok = 0;
   while ok~= 1
      gammatype = input('Use specified manifold pressure or flight-path angle ([m]/f)? ','s');
      if isempty(gammatype)
         gammatype = 'm';
      end
      if gammatype == 'f'
         gamma = input('Give flightpath angle [deg], default = 0: ')*pi/180;
         if isempty(gamma)
            gamma = 0;
            ok = 1;
         end
      elseif gammatype == 'm'
         ok = 1;
      else
         disp('   Invalid entry. Please choose either ''m'', or ''f''');
      end
   end

   thetadot = 0;
   psidot = 0;

   phidot = input('Give desired roll-rate [deg/s], default = 0: ')*pi/180;
   if isempty(phidot)
      phidot = 0;
   end

   if phidot ~= 0
      ok = 0;
      while ok~= 1
         rolltype = input('Roll in body or stability axes reference frame ([b]/s)? ','s');
         if isempty(rolltype)
            rolltype = 'b';
         end
         if rolltype == 'b' | rolltype == 's'
            ok = 1;
         else
            disp('   Invalid entry. Please choose either ''b'', or ''s''');
         end
      end
   else
      rolltype = 'b';   % No rolling, so for reasons of simplicity the default setting, 
                        % i.e. a body-axes roll, will be used.
   end

else
   % Set helpvariable skip = 1, to ensure that the aircraft configuration
   % does not have to be entered if the user chooses option 5 = QUIT.
   % --------------------------------------------------------------------
   skip = 1;
end

%%%%

if skip ~= 1         % DEFINE CONFIGURATION OF THE AIRPLANE, IF NOT QUITTING
                     % -----------------------------------------------------

   % For the 'Beaver' model, the flap angle and engine speed define the
   % configuration of the aircraft (engine speed is selected by the pilot
   % and maintained by the regulator of the propeller). Other aircraft
   % may have other definitions, like gear in/out, slats, settings of trim
   % surfaces, multiple engines, etc., so change the following lines if
   % required!
   % ---------------------------------------------------------------------
   deltaf = input('Give flap angle [deg], default = 0: ')*pi/180;
   if isempty(deltaf)
      deltaf = 0;
   end

   n = input('Give engine speed [RPM], default = 1800: ');
   if isempty(n)
      n = 1800;
   end

   % For the 'Beaver' aircraft, the engine power is determined by the engine
   % speed, which has been set above, and the manifold pressure, which will
   % be involved in the trim process if the flightpath angle gamma is user-
   % specified, or kept constant if pz is user-specified (defined in string-
   % variable gammatype). If the manifold pressure is to be adjusted by the
   % numerical trim algorithm, it is important to specify a meaningful esti-
   % mation of the manifold pressure, because otherwise the trim process may
   % not converge, or may give unrealistic results if the optimization pro-
   % cess converges to a LOCAL minimum.
   % -----------------------------------------------------------------------
   if gammatype == 'f'
      pz = 27; %input('Give estimate of manifold pressure pz ["Hg], default = 20: ');
   else  % gammatype == 'm'
      pz = input('Give manifold pressure pz ["Hg], default = 20: ');
   end
   if isempty(pz)
      pz = 20;
   end

%   % Hdot is the rate-of-climb, which is a function of the flightpath angle
%   % (change the program if you want to specify Hdot itself in stead of
%   % the flightpath angle gamma).
%   % ----------------------------------------------------------------------
%   Hdot = V*sin(gamma);

   % G is the centripetal acceleration, used for the coordinated turn
   % constraint.
   % ----------------------------------------------------------------
   G = psidot*V/9.80665;

   % If steady, uncoordinated, turning condition must be determined, the
   % roll angle can be specified freely; equilibrium will be obtained for
   % a trimmed-flight sideslip angle, which in this case usually will be
   % quite large.
   % --------------------------------------------------------------------
   phi = [];
   if opt == 2 & turntype == 'u';  % Steady turn, uncoordinated.
      phi = input('Give desired roll angle phi [deg], default = 0: ')*pi/180;
   end

   if isempty(phi)
      phi = 0;
   end

   % All variables which specify the flight condition will be put in the
   % vector ctrim (constants for trim process). The variables which are
   % varied by the minimization routine FMINS will be put into the vector
   % vtrim (variables for trim process). These vectors will have to be
   % redefined for aircraft which use other state or input vectors.
   %
   % The variable phi in ctrim is used for uncoordinated turns only; if a
   % coordinated turn is required, the coordinated turn constraint determines
   % the values of phi and beta. Type HELP ACCONSTR for more info.
   %
   % Definition:
   %
   % if gammatype =='f' (flight-path angle user-specified and manifold
   % pressure numerically adjusted):
   %    vtrim = [alpha beta deltae deltaa deltar pz]'
   % else
   %    vtrim = [alpha beta deltae deltaa deltar gamma]'
   %
   % Note: vtrim contains the first estimation of the non-constant states
   % and inputs. The final values will be determined iteratively.
   % -----------------------------------------------------------------------
   
   MAP = [];
   deltae = [];
   deltar = [];
   deltaa = [];
   xdotnostro=eye(0);

   for pz = 5:1:35;
   if gammatype == 'f' % gamma in ctrim, pz in vtrim
      ctrim = [V H psi gamma G psidot thetadot phidot deltaf n phi]';
      vtrim = [0 0 0 0 0 pz]';
   else % gamma in vtrim, pz in ctrim
      ctrim = [V H psi pz G psidot thetadot phidot deltaf n phi]';
      vtrim = [0 0 0 0 0 0]';  % <-- initial guess for gamma = 0!
   end

   % The Simulink system BEAVER or an equivalent aircraft model is called by
   % the function call xdot = feval('beaver',t,x,u,'outputs'), followed by
   % xdot = feval('beaver',t,x,u,'derivs') to obtain the time-derivatives of
   % the state vector x. Here t = 0 (system is time-independent). The
   % function call is created here, and stored in the variable modfun.
   % Note: the aircraft model itself will be compiled before applying these
   % function calls!
   % -----------------------------------------------------------------------
   modfun   = ['xdot = feval(''',sysname,''',0,x,u,''outputs''); '];
   modfun   = [modfun 'xdot = feval(''',sysname,''',0,x,u,''derivs'');'];

   % First pre-compile the aircraft model (temporarily suppress warnings
   % to prevent harmless, but inconvenient messages when the model is 
   % already compiled).
   %--------------------------------------------------------------------
   warning off
   feval(sysname,[],[],[],'compile');
   clear ans % (the above feval statement somehow generates an 'ans' variable)
   warning on

   % Call minimization routine FMINSEARCH. A scalar cost function is contained
   % in the M-file ACCOST, which also calls the constraints for coordinated
   % turns and rate-of-climb of the aircraft. FMINSEARCH returns the trim values
   % of alpha, beta, deltae, deltaa, deltar, and pz in the vector vtrimmed.
   % The convergence tolerance has been set to 1e-10 (default 1e-4), the 
   % maximum number of function evaluations has been set to 5000 (default 1200),
   % and the maximum number of iterations has been set to 3000 (default 1200)
   % ----------------------------------------------------------------------
   
   
%    clc;
%    disp('Searching for stable solution. Wait a moment...');
%    disp(' ');
   
   options = optimset('Display','off','TolX',1e-30,'MaxFunEvals',5e5,'MaxIter',5e5);
   [vtrimmed,fval,exitflag,output] = fminsearch('accost',vtrim,options,...
      ctrim,rolltype,turntype,gammatype,modfun)
   
   % Display error message when maximum number of iterations is 
   % reached before finding converged solution
   % ----------------------------------------------------------
   if exitflag == 0
      warning('Maximum number of iterations was exceeded!');
      disp(['- number of function evaluations: ' num2str(output.funcCount)]);
      disp(['- number of iterations: ' num2str(output.iterations)]);
   else
      disp('Converged to a trimmed solution...');
   end

%    disp(' ');
%    disp('<<< Press a key to get results >>>');
%    pause

   % Call subroutine ACCONSTR, which contains the flight-path constraints
   % once more to obtain the final values of the inputvector and states.
   % --------------------------------------------------------------------
   [x,u] = acconstr(vtrimmed,ctrim,rolltype,turntype,gammatype);
  
   MAP = [MAP u(6)];
   deltaa = [deltaa u(2)];
   deltar = [deltar u(3)];
   deltae = [deltae u(1)];
   
   % Call a/c model once more, to obtain the time-derivatives of the states
   % in trimmed-flight condition. By examining xdot, the user can decide if
   % the trimmed-flight condition is accurate enough. If not, either the
   % error tolerance of the minimization routine needs to be tightened, or
   % the cost function in ACCOST.M needs to be changed.
   % -----------------------------------------------------------------------
   eval(modfun);
   
   % Release compiled aircraft model now that all results are known
   % --------------------------------------------------------------
   feval(sysname,[],[],[],'term');

xdotnostro=[xdotnostro; xdot'];
   end
   % Display the final results. x, u, and xdot contain the inputvector,
   % state vector, and time-derivative of the state vector, respectively.
   % --------------------------------------------------------------------
%    clc
%    disp('State vector (trimmed): ');
%    x
%    disp(' ');
%    disp(' ');
%    disp('<<< Press a key >>>');
%    pause
% 
%    clc
%    disp('Input vector (trimmed): ');
%    u
%    disp(' ');
%    disp(' ');
%    disp('<<< Press a key >>>');
%    pause
%    
%    clc
%    disp('Time derivative of state vector (trimmed): ');
%    xdot
%    disp(' ');
%    disp(' ');
%    disp('<<< Press a key >>>');
%    pause

%    clc
%    disp('Vettore che serve a noi comuni mortali: ');
%    disp('FMprop, ypow')
%    ENGINE_Prop
%    disp(' ');
%    disp(' ');
%    disp('<<< Press a key >>>');
%    pause
   
   % The integrator block in the aircraft model uses the Matlab variable
   % xinco to store the initial value of x. The models from the library
   % EXAMPLES use the variables uaero0 and uprop0 to store the initial values
   % of the inputs to the aerodynamic and engine models, respectively.
   % Usually, trimmed-flight conditions are used as initial conditions.
   % Therefore, the trimmed-flight conditions will be stored in these'
   % variables, to ensure compatibility with the Simulink systems. The
   % trimmed-flight value of xdot = dx/dt will be stored in the variable
   % xdot0, which indicates that this value of xdot is obtained for the
   % INITIAL flight condition, defined by xinco, uaero0, and uprop0.
   %-------------------------------------------------------------------
%    xinco  = x;          % Initial value of state vector
%    xdot0  = xdot;       % Initial value of time-derivative of state vector
%    uaero0 = u(1:4);     % Initial value of input vector for aerodynamic model
%    uprop0 = u(5:6);     % Initial value of input vector for engine
                        %                            forces and moments model
   MAP
   deltae
   deltaa
   deltar
   xdotnostro
   % Now, a string-matrix will be created, which contains a description of
   % the trimmed-flight condition. Note: the function NUM2STR2 is just a
   % sligtly changed version of NUM2STR. Type HELP NUM2STR2 at the Matlab
   % command line for more info.
   % ---------------------------------------------------------------------
%    line1  = ['Trimmed-flight condition for S-function ' sysname ':'];
%    line2 = '';
% 
%    if opt == 1
%       line3 = ['Steady-state wings-level flight'];
%    elseif opt == 2
%       line3 = ['Steady-state turning flight'];
%    elseif opt == 3
%       line3 = ['Steady pull-up'];
%    elseif opt == 4
%       if rolltype == 'b'                            % Body-axis roll
%          line3 = ['Steady roll along the body-axis'];
%       else                                          % Stability-axis roll (default)
%          line3 = ['Steady roll along the stability-axis'];
%       end
%    end
% 
%    line4  = '';
%    line5  = ['User-specified definition of aircraft and flight condition:'];
%    line6  = '';
%    line7  = ['Flap angle:       deltaf    =  ' num2str2(deltaf,10) ' [deg]  '];
%    line8  = ['Engine speed:     n         =  ' num2str2(n,10)      ' [RPM]  '];
%    line9  = ['Airspeed:         V         =  ' num2str2(V,10)      ' [m/s]  '];
%    line10 = ['Altitude:         H         =  ' num2str2(H,10)      ' [m]    '];
%    line11 = ['Yaw-angle:        psi       =  ' num2str2(psi,10)    ' [deg]  '];
% 
%    if gammatype == 'f'
%       line12 = ['Flightpath angle: gamma     =  ' num2str2(gamma,10)  ' [deg]  '];
%    else
%       line12 = ['Manifold pressure: pz       =  ' num2str2(pz,10)     ' ["Hg]  '];
%    end
%    line13 = ['Yaw-rate:         psi dot   =  ' num2str2(psidot,10)   ' [deg/s]'];
%    line14 = ['Pitch-rate:       theta dot =  ' num2str2(thetadot,10) ' [deg/s]'];
%    line15 = ['Roll-rate:        phi dot   =  ' num2str2(phidot,10)   ' [deg/s]'];
% 
%    if opt == 2 & turntype == 'u' % uncoordinated turn
%       line16 = ['Uncoordinated turn,'];
%       line17 = ['Roll-angle:       phi       =  ' num2str2(phi,10) ' [deg]'];
%    elseif opt == 2 & turntype == 'c' % coordinated turn
%       line16 = ['Coordinated turn'];
%       line17 = '';
%    else % no turning flight
%       line16 = '';
%       line17 = '';
%    end
% 
%    % Some general remarks, need to be changed manually by the user!
%    % --------------------------------------------------------------
%    line18 = 'Definitions of state and input vectors as in system BEAVER:';
%    line19 = 'x = [V alpha beta p q r psi theta phi xe ye H]''';
%    line20 = 'u = [deltae deltaa deltar deltaf n pz uw vw ww uwdot vwdot wwdot]''';
%    line21 = '';
% 
%    % Add info about the variables in the Matlab Workspace.
%    % -----------------------------------------------------
%    line22 = 'Variables: xdot=initial state, uaero0=initial inputs to aeromodel';
%    line23 = 'uprop0=initial input to engine forces/moments model, xdot0=dx/dt(0)';
%    line24 = '';
% 
%    % Add time and date to text matrix.
%    % ---------------------------------
%    t  = clock;
%    t1 = num2str(t(4));
%    t2 = num2str(t(5));
%    t3 = num2str(t(6));
% 
%    line25 = ['date: ', date, '   time: ', t1, ':', t2, ':', t3];
%    
%    % Enhance explanatory lines to 80 characters and collect the results in the 
%    % string matrix 'trimdef'. Clear the individual line variables.
%    % -------------------------------------------------------------------------
%    trimdef = [];
%    for ii = 1:25
%       linestr = ['line' num2str(ii)];
%       eval([linestr '= [' linestr ' blanks(79-length(' linestr '))];']);
%       eval(['trimdef = [trimdef; ' linestr '];']);
%       eval(['clear ' linestr ]);
%    end
%    
%    % Results can now be saved to a file, which will be filled with the
%    % trim condition xinco, uaero0, uprop0, xdot0, and the explanation matrix
%    % trimdef.
%    %
%    % Note: files with steady-state trimmed flight condition have
%    % extension .TRI!
%    % -----------------------------------------------------------------
%    clc
%    answ = input('Save trimmed condition to file (y/[n])? ','s');
%    if isempty(answ)
%       answ = 'n';
%    end
%    if answ == 'y'
% 
%       % Destination directory for storing the results is the default FDC 
%       % data-directory, which is obtained from the function DATADIR.M. 
%       % If that directory does not exist, the current directory will be 
%       % used instead.
%       % ----------------------------------------------------------------
%       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 filename and path
%       % ------------------------
%       [filename,dirname] = uiputfile('*.tri','Save trimmed flight condition');
% 
%       % Save results
%       % ------------
%       save([dirname,filename],'xinco','uaero0','uprop0','xdot0','trimdef');
%      
%       % Display file location
%       % ---------------------
%       clc;
%       disp('Results saved to the file: ');
%       disp(' ');
%       disp(['     ',dirname,filename]);
%       disp(' ');
%       disp('<<< Press a key >>>');
%       pause
%    end
% 
%    % Display user-information.
%    % -------------------------
%    clc
%    disp('The results have been stored in the following variables:');
%    disp(' ');
%    disp('xinco = [V alpha beta p q r psi theta phi xe ye H]'' = state vector');
%    disp('xdot0 = dx/dt(0)');
%    disp('uaero0= [deltae deltaa deltar deltaf]'' = initial input vector for');
%    disp('                                          aerodynamic model');
%    disp('uprop0= [n pz]'' = initial input vector for engine model');
%    disp(' ');
%    disp('The text-matrix ''trimdef'' contains more info about the trimmed');
%    disp('flightcondition.');
%    disp(' ');
%    
%    % Repeat warning message if the maximum number of iterations was 
%    % reached before finding a converged solution
%    % --------------------------------------------------------------
%    if exitflag == 0
%    disp(' ');
%       disp('WARNING!');
%       disp('  The solution may not be accurate, as the maximum number of iterations'); 
%       disp('  was reached before finding a converged solution. You can increase the');
%       disp('  number of iterations by changing the MaxIter parameter of the simulation');
%       disp('  options. This involves editing actrim.m!');
%    end
% 
end                % END OF TRIM-LOOP (SKIPPED IF OPTION 5 = QUIT IS SELECTED)
% 
% 
% 
% % Clear unneccessary variables
% % ----------------------------
clear h ctrim vtrim vtrimmed modfun output options exitflag fval
clear opt V H psi gamma phidot psidot thetadot rolltype turntype gammatype
clear deltaf n pz Hdot G phi sysname x xdot u exitflag t t1 t2 t3 ii linestr
clear ok answ dirname filename savecmmnd defdir currentdir skip
disp('========================');
disp('         End, ready');
disp('========================');
disp(' ');
% disp(' ');
% disp('Ready.');
% disp(' ');
% 
% 
% %-----------------------------------------------------------------------
% % 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: September 7, 2005.  
% %-----------------------------------------------------------------------
