Added Part 3

TermProject_Part_III/TermProject_Part_III/4_6DoFControl/AircraftParameters.m

@@ -0,0 +1,203 @@
+%% Aircraft Model Parameters
+
+% defines the interfaces
+load BusObjects.mat
+
+% Constants for earth geoid model (used for position on earth)
+CNT_a = 6378137;
+CNT_f = 0.003352810664747;
+
+%% Trim Point
+
+%TrimPoint.States.debar = deg2rad(-6.8640);
+TrimPoint.States.debar = -0.14252;
+%TrimPoint.States.dabar = 0;;
+TrimPoint.States.dabar = -0.022319;
+%TrimPoint.States.drbar = 0;
+TrimPoint.States.drbar = 0.12421;
+%TrimPoint.States.dtbar = 0.185;
+TrimPoint.States.dtbar = 0.33985;
+
+%TrimPoint.States.X_P = 2000;
+TrimPoint.States.X_P = 5429.5159;
+
+%TrimPoint.States.Phibar   = 0;
+TrimPoint.States.Phibar   = 0.015575;
+%TrimPoint.States.Thetabar = deg2rad(4);
+TrimPoint.States.Thetabar = 0.10798;
+%TrimPoint.States.Psibar   = deg2rad(60);
+TrimPoint.States.Psibar   = 1.0472;
+
+%TrimPoint.States.Latbar   = deg2rad(53.58715);
+TrimPoint.States.Latbar   = 0.93527;
+%TrimPoint.States.Lonbar   = deg2rad(9.89862);
+TrimPoint.States.Lonbar   = 0.17276;
+%TrimPoint.States.Altbar   = 3000;
+TrimPoint.States.Altbar   = 4500;
+
+%TrimPoint.States.ubar = 120*cos(TrimPoint.States.Thetabar);
+TrimPoint.States.ubar = 115;
+%TrimPoint.States.vbar =   0;
+TrimPoint.States.vbar =   3.2615;
+%TrimPoint.States.wbar = 120*sin(TrimPoint.States.Thetabar);  
+TrimPoint.States.wbar = 12.4174;
+
+TrimPoint.States.pbar = 0;
+TrimPoint.States.qbar = 0;
+TrimPoint.States.rbar = 0;
+
+%% Inertia:
+
+% Tensor:
+
+Ixx = 105390; % inertia around body-fixed x axis in "kg*m^2"
+Iyy = 280790; % inertia around body-fixed y axis in "kg*m^2"
+Izz = 350200; % inertia around body-fixed z axis in "kg*m^2"
+Ixz = 15539;  % product of inertia in "kg*m^2"
+
+% Mass:
+
+CNT_g = [0;0;9.81]; % Gravitational acceleration vector 
+CNT_m = 15917;      % aircraft mass in "kg"
+
+%% Propulsion:
+
+Propulsion.Act_t = ss(-0.5450, 0.5450, 1, 0);  % engine dynamics
+Propulsion.Act_t.InputName = 'cmd_dt';
+Propulsion.Act_t.OutputName = 'X_Pi';
+Propulsion.Act_t.InputUnit = '-';
+Propulsion.Act_t.OutputUnit = 'N';
+
+Propulsion.Delays.tau_dt = 0.45;            % engine time delay in "s"
+Propulsion.Act_t.InputDelay = Propulsion.Delays.tau_dt;  % add time delay as input delay to ss model of engine
+
+% ---- Look-up table to compute thrust ----%
+Propulsion.LUT.LUT_dt  = [-0.8:0.1:1];  % in "-" (normalized)
+Propulsion.LUT.LUT_Alt = [2000;8000];   % in "m"
+Propulsion.LUT.LUT_V_A = [90;200];      % in "m/s"
+
+k = 1.04;   % temporary skaling factor
+Propulsion.LUT.LUT_X_P(:,:,1)=k.*[1640,-343;2816,1244];
+Propulsion.LUT.LUT_X_P(:,:,2)=k.*[1640,-343;2816,1244];
+Propulsion.LUT.LUT_X_P(:,:,3)=k.*[1640,-343;2816,1244];
+Propulsion.LUT.LUT_X_P(:,:,4)=k.*[2077,-85;2816,1244];
+Propulsion.LUT.LUT_X_P(:,:,5)=k.*[3386,1032;3475,1833];
+Propulsion.LUT.LUT_X_P(:,:,6)=k.*[4607,2042;4335,2610];
+Propulsion.LUT.LUT_X_P(:,:,7)=k.*[6044,3268;5061,3263];
+Propulsion.LUT.LUT_X_P(:,:,8)=k.*[7511,4452;5977,4119];
+Propulsion.LUT.LUT_X_P(:,:,9)=k.*[9174,5983;6870,4932];
+Propulsion.LUT.LUT_X_P(:,:,10)=k.*[10744,7271;7687,5669];
+Propulsion.LUT.LUT_X_P(:,:,11)=k.*[12458,8892;8518,6478];
+Propulsion.LUT.LUT_X_P(:,:,12)=k.*[14574,10806;9442,7423];
+Propulsion.LUT.LUT_X_P(:,:,13)=k.*[16780,12799;10179,8397];
+Propulsion.LUT.LUT_X_P(:,:,14)=k.*[18736,14701;10875,9215];
+Propulsion.LUT.LUT_X_P(:,:,15)=k.*[20883,16795;11509,9982];
+Propulsion.LUT.LUT_X_P(:,:,16)=k.*[22384,18750;12105,10744];
+Propulsion.LUT.LUT_X_P(:,:,17)=k.*[23705,20108;12712,11434];
+Propulsion.LUT.LUT_X_P(:,:,18)=k.*[25050,21432;13323,12079];
+Propulsion.LUT.LUT_X_P(:,:,19)=k.*[26243,22779;13933,12728];
+clear k
+% -----------------------------------------%
+
+Propulsion.rr_P1_f = [-1.448;-3;-0.533];   % location of left engine wrt. cg: [x,y,z] in "m"
+Propulsion.rr_P2_f = [-1.448; 3;-0.533];   % location of right engine wrt. cg: [x,y,z] in "m"
+Propulsion.CNT_kappa = 3*pi/180;        % engine inclination angle in "rad"
+
+%% Actuators:
+
+% initialize structs
+Actuators.Limits = [];
+Actuators.RateLimits = [];
+Actuators.Delays = [];
+
+%% Actuators:
+
+% elevator
+Actuators.Limits.u_max_He = [-0.3115,0.2332]; %maximum and minimum elevator deflection in "rad"
+Actuators.RateLimits.udot_max_He = [-0.85,0.85]; %maximum and minimum elevator deflection rate in "rad/s"
+
+Actuators.Act_e = ss(-21.2760, 21.2760, 1, 0);  % elevator dynamics
+Actuators.Act_e.InputName = 'cmd_de';
+Actuators.Act_e.OutputName = 'de';
+Actuators.Act_e.InputUnit = 'rad';
+Actuators.Act_e.OutputUnit = 'rad';
+Actuators.Act_e.StateName = 'de';
+Actuators.Act_e.StateUnit = 'rad';
+
+Actuators.Delays.tau_de = 0.033;    % time delay of elevator in "s"
+Actuators.Act_e.InputDelay = Actuators.Delays.tau_de; % add time delay as input delay to ss model of elevator
+
+% ---------------------------------------------------------------%
+
+% aileron
+Actuators.Limits.u_max_Ha = [-0.18,0.18]; %maximum and minimum aileron deflection in "rad"
+Actuators.RateLimits.udot_max_Ha = [-0.4,0.4]; %maximum and minimum aileron deflection rate in "rad/s"
+
+Actuators.Act_a = ss(-36.0248, 36.0248, 1, 0);  % aileron dynamics
+Actuators.Act_a.InputName = 'cmd_da';
+Actuators.Act_a.OutputName = 'da';
+Actuators.Act_a.InputUnit = 'rad';
+Actuators.Act_a.OutputUnit = 'rad';
+Actuators.Act_a.StateName = 'da';
+Actuators.Act_a.StateUnit = 'rad';
+
+Actuators.Delays.tau_da = 0.041;    % time delay of aileron in "s"
+Actuators.Act_a.InputDelay = Actuators.Delays.tau_da; % add time delay as input delay to ss model of aileron
+
+% ---------------------------------------------------------------%
+
+% rudder
+Actuators.Limits.u_max_Hr = [-0.407,0.407]; %maximum and minimum aileron deflection in "rad"
+Actuators.RateLimits.udot_max_Hr = [-0.8,0.8]; %maximum and minimum aileron deflection rate in "rad/s"
+
+Actuators.Act_r = ss(-15.6734,  15.6734,  1, 0);  % rudder dynamics
+Actuators.Act_r.InputName = 'cmd_dr';
+Actuators.Act_r.OutputName = 'dr';
+Actuators.Act_r.InputUnit = 'rad';
+Actuators.Act_r.OutputUnit = 'rad';
+Actuators.Act_r.StateName = 'dr';
+Actuators.Act_r.StateUnit = 'rad';
+
+Actuators.Delays.tau_dr = 0.044;    % time delay of aileron in "s"
+
+Actuators.Act_r.InputDelay = Actuators.Delays.tau_dr; % add time delay as input delay to ss model of aileron
+
+%% Aerodynamics:
+
+Aerodynamics.S = 62.4;      % reference area in "m^2"
+Aerodynamics.cbar = 3.36;   % mean aerodynamic chord in "m"
+
+Aerodynamics.LUT_rho = [1.16438645958276, 1.05543269917814, 0.954457178414045, 0.861045612606227, 0.774795981693757	0.695318454434115, 0.622235311201119, 0.555180865323805, 0.493801382899880, 0.437755001012510, 0.386711644273802, 0.340352939612400];
+Aerodynamics.LUT_Alt = [0:1000:11000];
+
+Aerodynamics.Ca_L0 = 0.071;  % lift 
+Aerodynamics.Ca_La = 4.55;
+Aerodynamics.Ca_Lq = 8.27; 
+Aerodynamics.Ca_Lde = 0.406;
+
+Aerodynamics.Ca_D0 = 0.01;  % drag
+Aerodynamics.Ca_Da = 0.22;
+Aerodynamics.Ca_Da2 = 0.67;
+Aerodynamics.Ca_Db2 = 0.041; % (LDM only)
+
+Aerodynamics.Ca_Yb = -0.82;  % side force (LDM only)
+Aerodynamics.Ca_Ydr = 0.126;
+Aerodynamics.Ca_Yp = 0;
+Aerodynamics.Ca_Yr = 0.0128;
+
+Aerodynamics.Ca_lb = -0.43; % roll moment (LDM only)
+Aerodynamics.Ca_lp = -11.53; 
+Aerodynamics.Ca_lr = 4.29;
+Aerodynamics.Ca_lda = -0.44;
+Aerodynamics.Ca_ldr = 0.051;
+
+Aerodynamics.Ca_m0 = -0.14; % pitching moment
+Aerodynamics.Ca_ma = -1.13;
+Aerodynamics.Ca_mq = -9.94;
+Aerodynamics.Ca_mde = -1.88;
+
+Aerodynamics.Ca_nb = 0.75;  % yaw moment (LDM only)
+Aerodynamics.Ca_np = 0.30; 
+Aerodynamics.Ca_nr = -6.68;
+Aerodynamics.Ca_nda = 0.023;
+Aerodynamics.Ca_ndr = -0.167;

TermProject_Part_III/TermProject_Part_III/4_6DoFControl/MyTrimStates.m

@@ -0,0 +1,15 @@
+%% Trim specification
+
+% all trim states refer to straight and level flight, i.e., the aircraft maintains altitude and attitude as well as its velocity  
+
+MyTrimState_1.u     = 115;      % m/s
+MyTrimState_1.q     = 0;        % rad/s
+MyTrimState_1.Alt   = 4500;     % m
+
+MyTrimState_2.q     = 0;        % rad/s
+MyTrimState_2.Alt   = 8000;     % m
+MyTrimState_2.Theta = 0.092;    % rad
+
+MyTrimState_3.q     = 0;        % rad/s
+MyTrimState_3.Alt   = 3000;     % m
+MyTrimState_3.dt    = 0.45;     % -

TermProject_Part_III/TermProject_Part_III/4_6DoFControl/checklinearmodel.m

@@ -0,0 +1,37 @@
+%% function to check step responses of a linear model
+% checklinearmodel(model,t)
+% inputs:   model -> linear system under investigation
+%           t     -> time scale
+% outputs:  figure only
+
+% This function is provided as part of the term project for the class 
+% "Flight Control Law Design and Application" at TU Hamburg.
+% by Julian Theis and Nicolas Sedlmair, 2021
+
+function [] = checklinearmodel(model,t)
+
+% get step responses
+[sig,time] = step(model*blkdiag(deg2rad(-3),deg2rad(3),deg2rad(5),0.2),t); % the inputs deviations are scaled to -3° elevator deflection, 3° aileron deflection, 5° rudder deflection and 0.2 thrust lever command
+
+% longitudinal motion
+figure(400); clf
+    subplot(3,2,1), plot(time,sig(:,4,1),'LineWidth',2), grid on, title('Input: $\Delta \delta_e$',"Interpreter","latex","FontSize",14), xlabel(''), ylabel('Output: $\tilde{u}$', "Interpreter","latex","FontSize",13)  
+    subplot(3,2,2), plot(time,sig(:,4,4),'LineWidth',2), grid on, title('Input: $\Delta \delta_t$',"Interpreter","latex","FontSize",14), ylabel('')   
+    subplot(3,2,3), plot(time,sig(:,2,1),'LineWidth',2), grid on, title(''), ylabel('Output: $\tilde{q}$', "Interpreter","latex","FontSize",13)
+    subplot(3,2,4), plot(time,sig(:,2,4),'LineWidth',2), grid on, title(''), ylabel('')
+    subplot(3,2,5), plot(time,sig(:,15,1),'LineWidth',2), grid on, title(''), xlabel('Time (s)'), ylabel('Output: $\tilde{\alpha}$', "Interpreter","latex","FontSize",13)
+    subplot(3,2,6), plot(time,sig(:,15,4),'LineWidth',2), grid on, title(''), xlabel('Time (s)'), ylabel('')
+    annotation('textbox', [0.35 0.99 0.1 0.02],'String','(Scaled) Step Inputs','LineStyle','none','Interpreter','latex','FontSize',13,'FitBoxToText','on');
+
+% lateral directional motion
+figure(401); clf
+    subplot(3,2,1), plot(time,sig(:,1,2),'LineWidth',2), grid on, title('Input: $\Delta \delta_a$',"Interpreter","latex","FontSize",14), xlabel(''), ylabel('Output: $\tilde{p}$', "Interpreter","latex","FontSize",13)  
+    subplot(3,2,2), plot(time,sig(:,1,3),'LineWidth',2), grid on, title('Input: $\Delta \delta_r$',"Interpreter","latex","FontSize",14), xlabel(''),
+    subplot(3,2,3), plot(time,sig(:,7,2),'LineWidth',2), grid on, title(''), ylabel('Output: $\tilde{\Phi}$', "Interpreter","latex","FontSize",13)
+    subplot(3,2,4), plot(time,sig(:,7,3),'LineWidth',2), grid on, title(''), ylabel('')
+    subplot(3,2,5), plot(time,sig(:,3,2),'LineWidth',2), grid on, title(''), xlabel('Time (s)'), ylabel('Output: $\tilde{r}$', "Interpreter","latex","FontSize",13)
+    subplot(3,2,6), plot(time,sig(:,3,3),'LineWidth',2), grid on, title(''), xlabel('Time (s)'), ylabel('')
+    annotation('textbox', [0.35 0.99 0.1 0.02],'String','(Scaled) Step Inputs','LineStyle','none','Interpreter','latex','FontSize',13,'FitBoxToText','on');
+
+%------------%
+clear sig time

TermProject_Part_III/TermProject_Part_III/4_6DoFControl/checkloopshape.m

@@ -0,0 +1,55 @@
+%% function to create magnitude plots of the loop transfer function
+% checkloopshape(L,Llbd,wlbd,Lubd,wubd)
+% inputs:   L -> loop transfer function (must be siso!)
+%           Llbd -> in db, lower bound for the loop gain at low frequencies
+%           wlbd -> in rad/s, frequency below which |L| shall be larger than |Llbd|
+%           Lubd -> in db, upper bound for the loop gain at high frequencies
+%           wubd -> in rad/s, frequency above which |L| shall be less than |Lubd|
+% outputs:  figure only
+
+% This function is provided as part of the term project for the class
+% "Flight Control Law Design and Application" at TU Hamburg.
+% by Julian Theis and Nicolas Sedlmair, 2021
+
+function []=checkloopshape(L,Llbd,wlbd,Lubd,wubd)
+
+% ------ preparations ----------
+Nw = 1000;
+w = logspace(-3,2,Nw)';
+
+Lmag = bode(L,w);
+Lmag = Lmag(:);
+LmagdB = 20*log10(Lmag);
+[~,Lidx] = min(abs(LmagdB));
+wL = w(Lidx);
+
+idx = (w >= wubd);
+idxwubd = find(idx,1,'first');
+
+xmin = 1e-3;
+xmax = 100;
+ymin = -50;
+ymax = 50;
+% ------------------------------
+
+lhL = semilogx([wL wL xmin],[ymin 0 0],'k');
+hold on;
+phLclb  = patch([wL/sqrt(10) wL  wL*sqrt(10) wL*sqrt(10) wL wL/sqrt(10)],...
+    [10 0 -15 ymin ymin ymin],'g');
+phLcub  = patch([wL/sqrt(10) wL  wL*sqrt(10) wL*sqrt(10) wL wL/sqrt(10)],...
+    [15 0 -10 ymax ymax ymax],'g');
+semilogx([wL/sqrt(10) wL  wL*sqrt(10)],[10 0 -15],'g');
+semilogx([wL/sqrt(10) wL  wL*sqrt(10)],[15 0 -10],'g');
+semilogx(w,LmagdB,'k','LineWidth',2);
+axis([xmin xmax ymin ymax]);
+phLlb = patch([w(1) wlbd wlbd w(1)],[ymin ymin Llbd Llbd],'b');
+phLub = patch([w(end) wubd wubd w(end)],[ymax ymax Lubd LmagdB(idxwubd)-20],'r');
+semilogx([w(1) wlbd],[Llbd Llbd],'b');
+semilogx([wubd w(end)],[Lubd Lubd-20],'r');
+set(phLlb,'FaceAlpha',0.15,'LineStyle','none');
+set(phLub,'FaceAlpha',0.15,'LineStyle','none');
+set(phLclb,'FaceAlpha',0.15,'LineStyle','none');
+set(phLcub,'FaceAlpha',0.15,'LineStyle','none');
+xlabel('Frequency, rad/sec');
+ylabel('|L|=|GK|, dB');
+hold off

TermProject_Part_III/TermProject_Part_III/4_6DoFControl/lat_compare.m

@@ -0,0 +1,41 @@
+%% function to compare the aircraft's response for different stages of the LDM SAS
+% [] = lat_compare(sys1,sys2,sys3)
+% inputs:   sysx -> systems to be compared
+% outputs:  figure only
+
+% This function is provided as part of the term project for the class
+% "Flight Control Law Design and Application" at TU Hamburg.
+% by Julian Theis and Nicolas Sedlmair, 2021
+
+function [] = lat_compare(sys1,sys2,sys3)
+
+ns = nargin();
+
+switch ns
+    case 2        
+        [sig,t] = step(sys1({'p','r','beta'},{'cmd_da','cmd_dr'})*blkdiag(deg2rad(5),deg2rad(5)),10);
+        [sig2,t2] = step(sys2({'p','r','beta'},{'cmd_da','cmd_dr'})*blkdiag(deg2rad(5),deg2rad(5)),10);
+
+        subplot(3,2,1), plot(t,180/pi.*sig(:,1,1),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,1,1),'LineWidth',2), grid, title('Step Response from input: $\Delta \delta_a$',"Interpreter","latex","FontSize",14), ylabel('Output: $\tilde{p}$ in deg/s', "Interpreter","latex","FontSize",10), legend('w/o yaw damper','with yaw damper'), ...
+        subplot(3,2,2), plot(t,180/pi.*sig(:,1,2),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,1,2),'LineWidth',2), grid, title('Step Response from input: $\Delta \delta_r$',"Interpreter","latex","FontSize",14), ...        
+        subplot(3,2,3), plot(t,180/pi.*sig(:,2,1),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,2,1),'LineWidth',2), grid, ylabel('Output: $\tilde{r}$ in deg/s', "Interpreter","latex","FontSize",10), ...
+        subplot(3,2,4), plot(t,180/pi.*sig(:,2,2),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,2,2),'LineWidth',2), grid, ...
+        subplot(3,2,5), plot(t,180/pi.*sig(:,3,1),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,3,1),'LineWidth',2), grid, ylabel('Output: $\tilde{\beta}$ in deg', "Interpreter","latex","FontSize",10), xlabel('Time in s'), ...
+        subplot(3,2,6), plot(t,180/pi.*sig(:,3,2),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,3,2),'LineWidth',2), grid, xlabel('Time in s'),
+
+    case 3
+        [sig,t] = step(sys1({'p','r','beta'},{'cmd_da','cmd_dr'})*blkdiag(deg2rad(5),deg2rad(5)),10);
+        [sig2,t2] = step(sys2({'p','r','beta'},{'cmd_da','cmd_dr'})*blkdiag(deg2rad(5),deg2rad(5)),10);
+        [sig3,t3] = step(sys3({'p','r','beta'},{'cmd_da','cmd_dr'})*blkdiag(deg2rad(5),deg2rad(5)),10);
+
+        subplot(3,2,1), plot(t,180/pi.*sig(:,1,1),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,1,1),'LineWidth',2), plot(t3,180/pi.*sig3(:,1,1),'LineWidth',2), grid, title('Step Response from input: $\Delta \delta_a$',"Interpreter","latex","FontSize",14), ylabel('Output: $\tilde{p}$ in deg/s', "Interpreter","latex","FontSize",10), legend('w/o yaw damper','with yaw damper','with ARI'), ...
+        subplot(3,2,2), plot(t,180/pi.*sig(:,1,2),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,1,2),'LineWidth',2), plot(t3,180/pi.*sig3(:,1,2),'LineWidth',2), grid, title('Step Response from input: $\Delta \delta_r$',"Interpreter","latex","FontSize",14), ...        
+        subplot(3,2,3), plot(t,180/pi.*sig(:,2,1),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,2,1),'LineWidth',2), plot(t3,180/pi.*sig3(:,2,1),'LineWidth',2), grid, ylabel('Output: $\tilde{r}$ in deg/s', "Interpreter","latex","FontSize",10), ...
+        subplot(3,2,4), plot(t,180/pi.*sig(:,2,2),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,2,2),'LineWidth',2), plot(t3,180/pi.*sig3(:,2,2),'LineWidth',2), grid, ...
+        subplot(3,2,5), plot(t,180/pi.*sig(:,3,1),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,3,1),'LineWidth',2), plot(t3,180/pi.*sig3(:,3,1),'LineWidth',2), grid, ylabel('Output: $\tilde{\beta}$ in deg', "Interpreter","latex","FontSize",10), xlabel('Time in s'), ...
+        subplot(3,2,6), plot(t,180/pi.*sig(:,3,2),'LineWidth',2), hold on, plot(t2,180/pi.*sig2(:,3,2),'LineWidth',2), plot(t3,180/pi.*sig3(:,3,2),'LineWidth',2), grid, xlabel('Time in s'),
+
+end
+
+clear t t2 t3 sig si2 sig3
+

TermProject_Part_III/TermProject_Part_III/4_6DoFControl/mode_resp_lat.m

@@ -0,0 +1,25 @@
+%% function to create mode response
+% mode_resp_lat(model,idx,ylimits)
+% inputs:   model   -> linear system under investigation
+%           idx     -> indizes corresponding to the the dynamic mode of interest, e.g. 'idx_dr' for the dutch roll mode
+%           ylimits -> limit the y-axis of the plot, e.g. [-0.1 0.9]
+% outputs:  figure only
+
+% This function is provided as part of the term project for the class 
+% "Flight Control Law Design and Application" at TU Hamburg.
+% by Julian Theis and Nicolas Sedlmair, 2021
+
+function [] = mode_resp_lat(model,idx,ylimits)
+
+[M,~] = eig(model.a); % compute the modal matrix of the linear system under investigation (columns are the right eigenvectors of the system)
+
+[X,T] = initial(model,M(:,idx(1))); % compute free response of the model: note that the first entry of 'idx' determines the dynamic mode as it defines the (right) eigenvector of the system to investigate
+X = real(X(:,[7 2 3 4])); % take only real part of the outputs of interest, i.e. the vector [7 2 3 4] picks the desired outputs of the model (see 'OutputNames' of the model) and changes the order
+X = X.*[180/pi 180/pi 180/pi 180/pi]; % transform 'rad' --> '°' and 'rad/s' --> '°/s'
+plot(T,X,'-','LineWidth',3);
+legend('beta (°)', 'p (°/s)', 'r (°/s)', '\Phi (°)','Location','best')
+xlabel('Time (s)')
+ylim(ylimits)
+grid
+clear X T M
+end