% Fall 2006 16.06 - 16.07 MATLAB and Simulink Tutorials
% Instructor: Violeta Ivanova, violeta@mit.edu
% 16.06 TA: Thomas Gray 
% 16.07 TA: Shannon Cheng

% Example 1-B: Longitudinal Response Modes for the F-8 Aircraft

% B1. State variables
% q(t)      = pitch rate (deg/s)
% u(t)      = perturbation from horizontal velocity (ft/s)
% alpha(t)  = perturbed angle of attack from trim (deg)
% theta(t)  = perturbed pitch angle from trim (deg)

% B2. Control variable
% delta(t)  = elevator deflection from trim (deg) unit input

global deltaE 
deltaE = 5*pi/180.;

% B3. State space: linearized longitudinal equations
% x' = Ax + Bv
% y  = Cx + Dv
% x = [q; u; alpha; theta]
% v = [delta]
% y = x

global A B

% Define matrices A, B, C, D for the F-8 aircraft

% CREATE MATRICES A AND B BASED ON THE F-8 EQUATIONS FROM YOUR HANDOUT:








C = eye(4);
D = zeros(4, 1);

% B4. Modal solutions - Short Period mode and Phugoid Mode
% x' = Ax

% COMPUTE EIGENVALUES AND EIGENVECTORS OF MATRIX A (one line of code):


      
% B5. Compute transfer functions for state variables from the state space

% COMPUTE NUMERATOR AND DENUMERATOR FOR THE TRANSFER FUNCTION
% FROM THE MATRICES A, B, C, AND D (one line of code):




% B6. Time interval and initial values of state variables

t0 = 0;
tf = 1000;
X0 = D;

% B7. Compute response to unit input of elevator deflection from trim

options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6, 'InitialStep', 1., 'MaxStep', 1.);

% COMPUTE THE ELEMENTS OF X VS. TIME T WITH THE ODE45 SOLVER, USING FUNCTION 
% LONGTIMERESPONSE, AND CONDITIONS & OPTIONS DEFINED ABOVE (one line of code):




% B8. Plot results

% Pitch rate perturbation from trim
subplot(221); 
plot(t, X(:,1)); 
xlabel('Time(sec)');
ylabel('q(t) [deg/s]');
grid;

% Velocity perturbation from trim
subplot(222); 
plot(t, X(:,2)); 
xlabel('Time(sec)');
ylabel('u(t) [ft/s]');
grid;

% Angle of attack perturbation from trim
subplot(223); 
plot(t, X(:,3)); 
xlabel('Time(sec)');
ylabel('alpha [deg]');
grid;

% Pitch angle perturbation from trim

% ADD A SUBPLOT FOR PUTCH ANGLE (THETA) VS. TIME: