%	SURVEY --  Stability and Control of a Linear (LTI) Aircraft Model
%   2021
%	July 5, 2023  
%	===============================================================
%	Copyright 2023 by ROBERT F. STENGEL.  All rights reserved.
    clear
    disp('  ')
    disp('**======================================================**')
    disp('** SURVEY of Linear, Time-Invariant (LTI) Dynamic Model **')
    disp('**======================================================**')
    disp(' ')
    disp(['Date and Time are ', num2str(datestr(now))]);
    disp('===================================')
    
%	TABLE OF CONTENTS
%   =================
%   1.  USER INPUT
%       .   Load Linear, Time-Invariant Dynamic Model (from FLIGHTv2.m)
%       A.  ANALYSIS FLAGS
%       B.  SELECTION FLAGS
%   2.  NOMINAL FLIGHT CONDITION [Input from FmodelFile.mat]
%   3.  Specification of SYSTEM MODEL DIMENSION AND AXES
%   4.  NATURAL FREQUENCIES, DAMPING RATIOS, and REAL ROOTS using NatFreq.m
%	5.  EIGENVALUES and EIGENVECTORS using StabMode.m
%   6.  INITIAL CONDITION RESPONSE
%   7.  SINGLE-INPUT/SINGLE-OUTPUT (SISO) ANALYSES
%       A.  CONTROLLABILITY and OBSERVABILITY with Single Input or Output
%       B.  Response to Control STEP INPUT
%       C.  Response to 0.5-s CONTROL IMPULSE
%       D.  Response to I.C. and CONTROL PROFILE
%       E.  EQUILIBRIUM RESPONSE to Control Perturbations
%   8.  TRANSFER FUNCTION ANALYSES
%       A.  PILOT Transfer Function
%       B.  AIRCRAFT Transfer Functions
%       C.  BODE Plots
%       D.  NYQUIST Plots
%       E.  NICHOLS Charts
%       F.  ROOT LOCUS Plots
%   9.  NOMENCLATURE
%       ==================================================================   
        D2R =   pi/180; %   Degrees to Radians
        R2D =   180/pi; %   Radians to Degrees
%       Output Matrices, Hx and Hu (called C and D elsewhere):
        Hx  =   eye(12);
        Hu  =   zeros(12,7);
        
%       =================================================================
%   1.  USER INPUT
%       ==========
        disp('1.   USER INPUT')
        disp('     ==========')
%       LOAD LTI generated and saved by FLIGHT.m
%       ========================================
%       Specify Default Full-State (12 x 12) and and Control (12 x 7) Model 
        load FmodelFile
        load GmodelFile
%       =================================================================
%       FmodelFile and GmodelFile contain the (12 x 12), (12 x 7), and (12 x 3) 
%       Jacobian Matrices Fmodel, Gmodel, (called A and B elsewhere) and
%       Lmodel (disturbance effects).  
%       The Aircraft's Stability, Control, and Disturbance Response 
%       Derivatives are expressed in Metric Units with Radian Measure.
%       FmodelFile also contains the Trim Conditions: 
%           TASms, hm, TrimAlphaDeg, and AeroModel
%       Output Matrices, Hx and Hu (called C and D elsewhere):
        Hx  =   eye(12);
        Hu  =   zeros(12,7);
        
%       1A. ANALYSIS FLAGS 
%           ==============
            NF      =	'Y';	%   Natural Frequencies, Damping Datios
                                %   'Y':  ON
                                %   'N':  OFF 
            EVEC    =	'N';	%   Eigenvalues and Eigenvectors
                                %   'Y':  ON
                                %   'N':  OFF 
            IC      =	'Y';	%   Time Response to Initial Condition
                                %   'Y':  ON
                                %   'N':  OFF 
            CO      =	'N';	%   Controllability and Observability
                                %   'Y':  ON
                                %   'N':  OFF 
            STEP    =	'N';	%   Time Response to Control Step
                                %   'Y':  ON
                                %   'N':  OFF 
            CONPR   =   'N';    %   Time Response to Control Profile 
                                %   'N':  OFF 
                                %   'I':  ON, 0.5-sec Impulse
                                %   'U':  ON, User-specified Profile
            EQUIL   =	'N';	%   Equilibrium Response to Fixed Control 
                                %   'Y':  ON
                                %   'N':  OFF 
            TRF     =	'Y';	%   Aircraft Transfer Function
                                %   'Y':  ON (Required for FREQUENCY-DOMAIN
                                %         analysis)
                                %   'N':  OFF 
            PILOT   =   'N';    %   Human-Pilot Transfer Function
                                %   'Y':  ON (Required for PILOT-AIRCRAFT
                                %         analysis)
                                %   'N':  OFF 
            STAB    =   'N';    %   Stability Margins for SISO Unit Feedback, 
                                %   Positive and Negative Gains.  
                                %   Gm(linear), Pm(deg), Wcg(rad/s),
                                %   Wcp(rad/s)
                                %   'Y':  ON
                                %   'N':  OFF 
            BODE    =	'N';	%   BODE Plot
                                %   'Y':  ON
                                %   'A':  ON, with amplitude ratio and phase
                                %         angle aysmptotes (requires asymp.m)
                                %   'N':  OFF 
            NYQUIST =	'N';	%   NYQUIST Chart
                                %   'Y':  ON
                                %   'N':  OFF 
            NICHOLS =	'N';	%   NICHOLS Chart
                                %   'Y':  ON
                                %   'N':  OFF 
            RLOCUS  =   'N';    %   ROOT LOCUS Plot
                                %   'Y':  ON
                                %   'N':  OFF 
                            
            tSim    =   20;     %   Time Duration for LTI Simulation, sec
            tInt    =   0.1;    %   Sampling Interval for LTI Simulation, sec
                            
            
%       1B. SELECTION FLAGS
%           ===============
            ORIG    =   'Y';    %   Show Original Model & E-values
                                %   'Y':  Yes
                                %   'N':  No 
            ROM     =	'Y';    %   State Order
                                %   'N':  Original State Order 
                                %   'Y':  Re-order State
            LonLD   =	'Lon';	%   Decoupled Reference Frame
                                %   'Lon':  Longitudinal
                                %   'LaD':  Lateral-Directional 
            SBA     =	'H';	% 	Original or Transformed State Elements
                                %   'B':  Body Axes 
                                %   'H':  Hybrid/Stability Axes                        
            DIM     =	4;	    % 	Model Dimension
                                %   12: 12th-order Coupled Model 
                                %   6:  6th-order Decoupled Model
                                %   4:  4th-order Decoupled Model
 %      =================================================================
 
        disp(['     Natural Frequencies: ',num2str(NF)])
        disp(['     EigenValues/Vectors: ',num2str(EVEC)])
        disp('  ')
        disp(['     Controllability/Observability: ',num2str(CO)])
        disp(['     Initial-Condition Response: ',num2str(IC)])
        disp(['     Step Response = ',num2str(STEP)])
        disp(['     Control-Profile Response: ',num2str(CONPR)])
        disp(['     Equilibrium Response: ',num2str(EQUIL)])
        disp('  ')
        disp(['     Transfer Functions: ',num2str(TRF)])
        disp(['     Human Pilot Dynamics: ',num2str(PILOT)])
        disp(['     Bode Plot: ',num2str(BODE)])
        disp(['     Nyquist Plot: ',num2str(NYQUIST)])
        disp(['     Nichols Chart: ',num2str(NICHOLS)])
        disp(['     Root Locus: ',num2str(RLOCUS)])
        disp('     --------------------------------------------------------')
        disp(['     AeroModel is "',MODEL,'"'])
        disp(['     Show Original Model: ',num2str(ORIG)])
        disp(['     Re-Order Model: ',num2str(ROM)])
        disp(['     Lon/LatDir: ',num2str(LonLD)])
        disp(['     Model Order: ',num2str(DIM)])
        disp(['     Body or Hybrid/Stability Axes: ',num2str(SBA)])
        disp('     ========================================================')
    
%   2.  Nominal Flight Condition INPUT [from FmodelFile.mat]
%       ====================================================
        TAS         =   TASms;          % True Air Speed, TAS, m/s	(relative to air mass)
        hTrim       =   hm;             % Altitude above Sea Level, m
        AlphaTrim   =   TrimAlphaDeg;   % Trim Angle of Attack, deg
        BetaTrim    =   0;              % Trim Sideslip Angle, deg
            
        disp('  ')
        disp('2.   NOMINAL FLIGHT CONDITION [Input from FmodelFile, GmodelFile]')
        disp('     ========================')
        disp(['     TAS   = ',num2str(TAS),' m/s      Altitude = ',num2str(hTrim),' m'])
        disp(['     Alpha = ',num2str(AlphaTrim),' deg       Beta = ',num2str(BetaTrim),' deg'])
        Alpha0  =	AlphaTrim * D2R;
        Beta0   =	BetaTrim * D2R;

        if ORIG == 'Y'
            disp('  ')
            disp('  ====================================================')
            disp('  ORIGINAL Stability, Control, & Disturbance Variables')
            disp('     Body-Axis Velocity        del(u,v,w)') 
            disp('     Inertial Position         del(x,y,z)')
            disp('     Body-Axis Angular Rate    del(p,q,r)') 
            disp('     Euler Angles              del(Phi,Theta, Psi)')
            disp('     Controls                  del(dE,dA,dR,dT,dAS,dF,dS)')
            disp('     Disturbances              del(u_w,v_w,w_w)')
            F   =	Fmodel
            G   =	Gmodel
            L   =   Lmodel
            Sys =   ss(F,[G L],Hx,0);
            disp('=============================')
            disp('EIGENVALUES of ORIGINAL Model')
            E   =	eig(F)
            disp('=============================')
        end
        
%   3.  Specification of SYSTEM MODEL DIMENSION AND AXES
%       ================================================
        disp('  ')
        disp('3.    Specification of SYSTEM MODEL DIMENSION AND AXES')
        disp('      ================================================')
    			
%       ------------------
%       MODEL for ANALYSIS
%       ------------------
        [F,G,L]	=	LonLatDir(Fmodel,Gmodel,Lmodel,TAS,Beta0,Alpha0,ROM,LonLD,DIM,SBA);
    
        if DIM == 12
            Hx  =   eye(12);
        end
        if DIM == 6
            Hx  =   eye(6);
        end
        if DIM == 4
            Hx  =   eye(4);
        end
        Hu  =   0;

        Sys  =	ss(F,G,Hx,Hu);      %   Without Disturbance Effects
        SysL =  ss(F,[G L],Hx,Hu);  %   Disturbance Effects Included
        disp('  ')
        disp('==========================')
        disp('MODEL CHOSEN for ANALYSIS')
        disp('==========================')

        if LonLD == 'Lon'
            disp('  . LONGITUDINAL Axes')
        else
            disp('  . LATERAL-DIRECTIONAL AXES')
        end
        if SBA == 'B'
            disp('  . BODY reference frame')
        else
            disp('  . HYBRID reference frame')
        end
        if DIM == 6
            disp('  . 6th Order')
        end
        if DIM == 4
            disp('  . 4th Order')
        end
        disp('--------------------------')
        disp('  ')
        disp('Stability Matrix, F')
        F
        disp('  ')
        disp('Control Matrix, G')
        G
        disp('  ')
        disp('Disturbance Matrix, L')
        L
        disp('  ')
        disp('State Output Matrix, Hx')
        Hx
        disp('  ')
        disp('Control Output Matrix, Hu')
        Hu
    
%   4.  NATURAL FREQUENCIES, DAMPING RATIOS, and REAL ROOTS using NatFreq.m
%       ===================================================================
        if NF == 'Y'
            disp('4.    NATURAL FREQUENCIES, DAMPING RATIOS, and REAL ROOTS')
            disp('      ===================================================')
            nf	=	NatFreq(F,Sys);
        end
    
%	5.  EIGENVALUES and EIGENVECTORS using StabMode.m
%       =============================================
        if EVEC == 'Y'
            disp('  ')
            disp('5.     EIGENVALUES and EIGENVECTORS')
            disp('       ============================')
            sm	=	StabMode(F,TAS,ROM,LonLD,SBA,DIM);
        end
        
%   6.  INITIAL CONDITION RESPONSE
%       ==========================
        disp('  ')
        disp('6.    INITIAL CONDITION RESPONSE')
        disp('      ==========================')
        if IC == 'Y'
            disp('      See plots')
            delX0   =   zeros(DIM,1);
            t       = 	0:tInt:tSim;	% Time Vector for Initial-Condition Response, sec
            for i = 1:DIM
                delX0(i,1)    =   1;    % Scalar unit initial condition
                figure
                initial(Sys,delX0,t),grid,ylabel('State'),
                    title(['Initial-Condition Response, LonLD = ',...
                    num2str(LonLD),', DIM = ',num2str(DIM),', SBA = ',...
                    num2str(SBA),' x0(i) = ',num2str(i)])            
                delX0   =   zeros(DIM,1);
            end
        end

%   7.  SINGLE-INPUT/SINGLE-OUTPUT ANALYSES
%       ===================================
        disp('  ')
        disp('7.    SINGLE-INPUT/SINGLE-OUTPUT ANALYSES')
        disp('      ===================================')

%       7A. Controllability and Observability with Single Input or Output
%           =============================================================
            if CO == 'Y'
                disp('  ')
                disp('  7A. CONTROLLABILITY with Single Control Input')
                disp('      -----------------------------------------')
                if LonLD == 'LaD'
                    index   =   3;
                end
                if LonLD == 'Lon'
                    index   =   4;
                end
                for i = 1:index
                    Gcon    =   G(:,i)
                    Con		=	ctrb(F,Gcon);
                    unCon	=	length(F) - rank(Con);
                    disp(['Uncontrollable Elements: ',num2str(unCon)])
                end
                disp('  ')
                disp('  7A.   OBSERVABILITY with Single Output')
                disp('        --------------------------------')
                for i = 1:DIM
                    Hobs    =   Hx(i,:)
                    Obs		=	obsv(F,Hobs);
                    unObs	=	length(F) - rank(Obs);			
                    disp(['Unobservable Elements: ',num2str(unObs)])
                end
            end
        

%       7B. Response to Control STEP INPUTS
%           ===============================
            if STEP == 'Y'
                disp('  ')
                disp('  7B. Control STEP-INPUT RESPONSE')
                disp('      ---------------------------')
                disp('      See plots')
                t       = 	0:tInt:tSim;	% Time Vector for Step-Input Response, sec      
                tLength =   length(t);
                if LonLD == 'LaD'
                    delU0   =   zeros(tLength,3);
                end
                if LonLD == 'Lon'
                    delU0   =   zeros(tLength,4);
                end
                [M,N]   =   size(delU0);
                for i = 1:N
                    for j = 1:tLength
                        delU0(j,i)  =   1;  % Scalar unit control step
                    end
                    figure
                    lsim(Sys,delU0,t),grid,ylabel('State'),
                        title(['STEP-INPUT RESPONSE, LonLD = ',...
                        num2str(LonLD),', DIM = ',num2str(DIM),...
                        ', SBA = ',num2str(SBA),', Control = ',num2str(i)])            
                    delU0   =   zeros(tLength,N);
                end
            end
    
%       7C. Response to SCALAR CONTROL IMPULSES
%           ===================================
            if CONPR == 'I'
                disp('  ')
                disp('  7C. System Response to 0.5-s CONTROL IMPULSE')
                disp('      ----------------------------------------')
                disp('      See plots')
                t       = 	0:tInt:tSim;	% Time Vector, 0.5-sec Impulse, sec      
                tLength =   length(t);
                uPro    =   [1,1,1,1,1,zeros(1,(tLength - 5))]'; % Column
                if LonLD == 'LaD'
                    N   =   3;
                end
                if LonLD == 'Lon'
                    N   =   4;
                end
                delU0   =   zeros(tLength,N);
                for i = 1:N
                    for j = 1:tLength
                        delU0(j,i)  =   uPro(j);
                    end
                    figure
                    lsim(Sys,delU0,t),grid,ylabel('State'),
                        title(['0.5-sec IMPULSE RESPONSE, LonLD = ',...
                        'LonLD',' DIM = ',num2str(DIM),', SBA = ','SBA',...
                        ', Control = ',num2str(i)])            
                    delU0   =   zeros(tLength,N);
                end
                figure
                plot(t,uPro),grid, ylabel('Control Profile'), ...
                    xlabel('Time, s'), title('Control Profile')
            end

  %     7D. Multivariate Response to I.C. and CONTROL PROFILE
  %         =================================================    
            if CONPR == 'U'
                disp('  ')
                disp('7D.   System Response to I.C. and CONTROL PROFILE')
                disp('      -------------------------------------------')
                disp('      See plot')
                t       = 	0:tInt:tSim;	% Time for I.C. and Control Profile, sec      
                tLength =   length(t);
%               ======================================
                if LonLD == 'LaD'
                    N   =   3;
                end
                if LonLD == 'Lon'
                    N   =   4;
                end
%               USER INPUT of I.C. and Control Profile
%               --------------------------------------
%               REPLACE Default State and Control with DESIRED VALUES
%               -----------------------------------------------------
                delX0 =     zeros(DIM,1);
                delU0 =     zeros(tLength,N);
%               -----------------------------------------------------
                figure
                lsim(Sys,delU0,t,delX0),grid,ylabel('State'),
                     title(['I.C. and CONTROL PROFILE Response, LonLD = ',...
                     'LonLD',' DIM = ',num2str(DIM),', SBA = ',num2str(SBA)])            
            end

%       7E. EQUILIBRIUM RESPONSE to Control Perturbations
%           [4th-Order Models only]
%           =============================================
            if EQUIL == 'Y'
                disp('  7E. EQUILIBRIUM RESPONSE to Control Perturbation')
                if DIM ~= 4
                    disp('      No CONTROL Equilibrium: DIM ~= 4')
                end
                if DIM == 4
                    disp('      --------------------------------------------')

                    if det(F) ~= 0
                        disp('      4th-Order EQUILIBRIUM RESPONSE to Constant CONTROL Input')
                        uEq =   [];
                        if LonLD == 'LaD'
                            index   =   3;
                        end
                        if LonLD == 'Lon'
                            index   =   4;
                        end
                        uEq =   zeros(index,1);
                        for i = 1:index
                            uEq(i,1)    =   1;
                            xEq         =   -inv(F)*G*uEq
                            uEq
                            uEq         =   zeros(index,1);
                        end
                    else
                        disp('  ')
                        disp('  No CONTROL Equilibrium: F is Singular')
                    end
                end
            end
        
        disp('  ')
        disp('8.    TRANSFER FUNCTION ANALYSES')
        disp('      ==========================')

%   8.  TRANSFER FUNCTION ANALYSES
%       ==========================
%       8A. HUMAN PILOT Transfer Function
%           ================================
            if PILOT == 'Y'
                        disp('  ')
                        disp('8A.    HUMAN PILOT Transfer Function')
                        disp('       -----------------------------')
%           1st-Order Pilot Model, TP = Time Constant, sec
                TP          =   0.25;   % Pilot time constant, sec
                SysPilot    =   tf((1/TP),[1 1/TP])
            end
%       8B. AIRCRAFT Transfer Functions       
%           =========================== 
            if TRF == 'Y'
                disp('  ')
                disp(' 8B.  Single-Input/Single-Output AIRCRAFT Transfer Functions')
                disp('      ------------------------------------------------------')
                if LonLD == 'LaD'
                    N   =   3;
                end
                if LonLD == 'Lon'
                    N   =   4;
                end
                for i = 1:N
                    disp('  ')
                    disp(['TRANSFER FUNCTION, Input # = ',num2str(i)])
                    disp('-----------------  ============')
                    Gcon    =   G(:,i)
                    for j = 1:DIM
                        disp('  ')
                        disp(['OUTPUT # = ',num2str(j)])
                        disp('------------')
                        Hobs    =   Hx(j,:);
                        Sys     =	ss(F,Gcon,Hobs,0);
                        if PILOT == 'Y'
                            Sys     =   SysPilot*Sys;
                            ZPKsys  =	zpk(Sys)
                        end
                            DC_Gain	=	dcgain(Sys);
                            disp(['DC GAIN = ',num2str(DC_Gain)])
%                         if STAB == 'Y'                   
% %                           Stability Margins for SISO Unit Feedback, Positive 
% %                           and Negative Gains.  Gm(linear), Pm(deg), Wcg(rad/s),
% %                           Wcp(rad/s)
%                             [Gm,Pm,Wcg,Wcp] =   margin(ZPKsys)
%                             [Gm,Pm,Wcg,Wcp] =   margin(-ZPKsys)
%                         end
                    
%       8C. BODE Plot
%           =========
                        if  BODE == 'Y'
                            disp('i.   BODE Plot')
                            figure
                            bode(Sys),grid,
                                title(['LonLD = ',num2str(LonLD),', y(j) = ',num2str(j),...
                                ' u(i) = ',num2str(i)])
                        end
                        if  BODE == 'A'
                            if  DC_Gain ~= 0
                                disp('i.  BODE Plot with Asymptotes. T. Andresen, MATLAB Central, 2023')
                	            figure
                                asymp(Sys),grid on
                                    title(['BODE Plot with Asymptotes. LonLD = ',num2str(LonLD),', y(j) = ',...
                                    num2str(j),' u(i) = ',num2str(i)])
                            else
                                disp('i.  BODE Plot')
                                figure
                                bode(Sys),grid,
                                    title(['BODE Plot. LonLD = ',num2str(LonLD),', y(j) = ',...
                                    num2str(j),' u(i) = ',num2str(i)])
                            end
                        end
                    
%       8D. NYQUIST Plot
%           ============
                        if  NYQUIST == 'Y'
            	            disp('ii.  NYQUIST Plot')
                            ZPKsys  =	zpk(Sys)
                            if  DC_Gain ~= 0
                                ZPKabs	=	sign(DC_Gain) * ZPKsys;
                            else
                            ZPKabs	=	ZPKsys;
                            end
                            figure
                            nyquist(ZPKabs), grid
                                title(['NYQUIST Plot. LonLD = ',num2str(LonLD),', y(j) = ',num2str(j),...
                    	        ' u(i) = ',num2str(i)])
                        end

%           iii.    NICHOLS Chart
                        if NICHOLS == 'Y'
                            disp('iii.  NICHOLS Chart')
                            ZPKsys  =	zpk(Sys)
                            if DC_Gain ~= 0
                                ZPKabs	=	sign(DC_Gain) * ZPKsys; % Impose positivity
                            else
                            ZPKabs	=	ZPKsys;
                            end
                            figure
                               nichols(ZPKabs),grid
                                    title(['Nichols Chart. LonLD = ',num2str(LonLD),', y(j) = ',num2str(j),...
                                    ' u(i) = ',num2str(i)])
                        end
%           iv. ROOT LOCUS
                        if RLOCUS == 'Y'
                            disp('iv.  ROOT LOCUS')
                            figure
                            rlocus(Sys,-Sys),grid
                                title(['ROOT LOCUS. LonLD = ',num2str(LonLD),', y(j) = ',num2str(j),...
                                ' u(i) = ',num2str(i)])
                        end
                    end
                end
            end            

    AeroModel   =   MODEL;        
    save('ControlData.mat', 'F','G','L','Hx','Hu','LonLD','DIM','SBA', ...
        'TASms','hm','TrimAlphaDeg','AeroModel')
                
%   =====================================================================		
%   9.  NOMENCLATURE
%       Functions used by SURVEY
%           LonLatDir.m     Reduced-Order Models derived from Fmodel and Gmodel
%           NatFreq.m       Natural Frequencies, Damping Ratios, and Time Constants
%           StabMode.m      Modes of Motion [Eigenvalues and Eigenvectores]
%           Control.m       LTI Flight Control System Design

%       Loaded Data Files: 
%           FmodelFile.mat
%           GmodelFile.mat

%       Saved Data Files:
%           ControlData.mat

%	DEFINITION OF THE STATE VECTOR FOR Fmodel Rows and Columns 
%   and Gmodel Rows
%		x(1) = 		ub, m/s:    Body-axis x Inertial Velocity
%		x(2) =		vb, m/s:    Body-axis y Inertial Velocity 
%		x(3) =		wb, m/s:    Body-axis z Inertial Velocity 
%		x(4) =		xe, m:      North Position of Center of Mass (c.m.) WRT Earth 
%		x(5) =		ye, m:      East Position of Center of Mass WRT Earth 
%		x(6) =		ze, m:      Negative of c.m. Altitude (-h) WRT Earth 
%		x(7) =		p, rad/s:  Body-axis Roll Rate
%		x(8) =		q, rad/s:  Body-axis Pitch Rate
%		x(9) =		r,rad/s:   Body-axis Yaw Rate
%		x(10)=		Phi, rad:  Roll Angle of Body WRT Earth 
%		x(11)=		Pheta, rad:    Pitch Angle of Body WRT Earth 
%		x(12)=		Psi, rad:  Yaw Angle of Body WRT Earth 
	
%	DEFINITION OF THE CONTROL VECTOR FOR Gmodel Columns
%		u(1) = 		dEr, rad:   Elevator 
%		u(2) = 		dAr, rad:   Aileron
%		u(3) = 		dRr, rad:   Rudder 
%		u(4) = 		dT, %/100:  Throttle
%		u(5) =		dASr, rad:  Asymmetric Spoiler
%		u(6) =		dFr, rad:   Flap
%		u(7) =		dSr, rad:   Stabilator

%	DEFINITION OF THE DISTURBANCE VECTOR FOR Lmodel Columns
%		w(1) = 		u_w, m/s:   Body-axis Axial Wind Component 
%		w(2) = 		v_w, m/s:   Body-axis Side Wind Component
%		w(3) = 		w_w, m/s:   Body-axis Normal Wind Component

%	DEFINITION OF THE STATE, CONTROL, and DISTURBANCE VECTORs for 
%   REDUCED-ORDER MODELs

%   6th-Order Lateral-Directional Body-Axis Model
%       xLD = [v p r Phi Psi y]'
%       uLD = [dA dR dAS]'
%       wLD = [v_w]'

%   4th-Order Lateral-Directional Body-Axis Model
%       xLD = [v p r Phi]'
%       uLD = [dA dR dAS]'
%       wLD = [v_w]'

%   6th-Order Lateral-Directional Stability-Axis Model
%       xLD = [rS Beta pS Phi Psi y]'
%       uLD = [dA dR dAS]'
%       wLD = [Beta_w]'

%   4th-Order Lateral-Directional Stability-Axis Model
%       xLD = [rS Beta pS Phi]'
%       uLD = [dA dR dAS]'
%       wLD = [Beta_w]'
%   6th-Order Longitudinal Body-Axis Model
%       xLon = [u w q Theta x z]'
%       uLon = [dE dT dF dS]'
%       wLon = [u_w w_w]'

%   4th-Order Longitudinal Body-Axis Model
%       xLon = [u w q Theta]'
%       uLon = [dE dT dF dS]'
%       wLon = [u_w w_w]'

%   6th-Order Longitudinal Hybrid-Axis Model
%       xLon = [V Gamma q Alpha x z]'
%       uLon = [dE dT dF dS]'
%       wLon = [V_w Alpha_w]'

%   4th-Order Longitudinal Hybrid-Axis Model
%       xLon = [V Gamma q Alpha]'
%       uLon = [dE dT dF dS]'
%       wLon = [V_w Alpha_w]'

    disp('  ')
    disp('End of SURVEY Analysis of Open-Loop System')
    disp([' AeroModel is "',MODEL,'"'])
    disp('==========================================')

%   END of Program
