how to include variables for plotting

Hi, I have this function and script files. I would like to know how I could plot the control moments?? de, da, dr not sure how to do it. Thanks, any help will be gretly appreciated.

function xdot = STOL_EOM(t,x)
% STOL_EOM contains the nonlinear equations of motion for a rigid airplane.
% (NOTE: The aerodynamic model is linear.)
global e1 e2 e3 rho m g S b c AR WE Power ...
    GeneralizedInertia GeneralizedInertia_Inv ...
    CD0 e CL_Trim CLalpha CLq CLde Cmde ...
    CYbeta CYp CYr Cmalpha Cmq Clbeta Clp Clr Cnbeta Cnp Cnr ...
    CYda CYdr Clda Cldr Cnda Cndr ...
    de0 da0 dr0 
X = x(1:3);
Theta = x(4:6);
    phi = Theta(1);
    theta = Theta(2);
    psi = Theta(3);
V = x(7:9);
    u = V(1);
    v = V(2);
    w = V(3);    
omega = x(10:12);
    p = omega(1);
    q = omega(2);
    r = omega(3);    
alpha = atan(w/u);
beta = asin(v/norm(V));
P_dynamic = (1/2)*rho*norm(V)^2;
RIB = expm(psi*hat(e3))*expm(theta*hat(e2))*expm(phi*hat(e1));
LIB = [1, sin(phi)*tan(theta), cos(phi)*tan(theta);
    0, cos(phi), -sin(phi);
    0, sin(phi)/cos(theta), cos(phi)/cos(theta)];
VE = V + RIB'*WE;
% Kinematic equations
XDot = RIB*VE;
ThetaDot = LIB*omega;
% Weight
W = RIB'*(m*g*e3);
% Control moments
de = de0;
da = da0;
dr = dr0;
% Components of aerodynamic force (modulo "unsteady" terms)
CL = CL_Trim + CLalpha*alpha + CLq*((q*c)/(2*norm(V))) + CLde*de;
CD = CD0 + (CL^2)/(e*pi*AR);
temp = expm(-alpha*hat(e2))*expm(beta*hat(e3))*[-CD; 0; -CL];
CX = temp(1);
CZ = temp(3);
CY = CYbeta*beta + CYp*((b*p)/(2*norm(V))) + CYr*((b*r)/(2*norm(V))) + ...
    CYda*da + CYdr*dr;
%X = P_dynamic*S*CX + T0;  %Constant Thrust?
Thrust = Power/norm(V);
X = P_dynamic*S*CX + Thrust;  %Constant Thrust?
Y = P_dynamic*S*CY;
Z = P_dynamic*S*CZ;
Force_Aero = [X; Y; Z];
% Components of aerodynamic moment (modulo "unsteady" terms)
Cl = Clbeta*beta + Clp*((b*p)/(2*norm(V))) + Clr*((b*r)/(2*norm(V))) + ...
    Clda*da + Cldr*dr;
Cm = Cmalpha*alpha + Cmq*((c*q)/(2*norm(V))) + Cmde*de;
Cn = Cnbeta*beta + Cnp*((b*p)/(2*norm(V))) + Cnr*((b*r)/(2*norm(V))) + ...
    Cnda*da + Cndr*dr;
L = P_dynamic*S*b*Cl;
M = P_dynamic*S*c*Cm;
N = P_dynamic*S*b*Cn;
Moment_Aero = [L; M; N];
% Sum of forces and moments
Force = W + Force_Aero;
Moment = Moment_Aero;
% NOTE: GeneralizedInertia includes "unsteady" aerodynamic terms
% (i.e., added mass/inertia).
temp = GeneralizedInertia*[VE; omega];
LinearMomentum = temp(1:3);
AngularMomentum = temp(4:6);
clear temp
RHS = [cross(LinearMomentum,omega) + Force; ...
    cross(AngularMomentum,omega) + Moment];
temp = GeneralizedInertia_Inv*RHS;
VEDot =  temp(1:3);
VDot = VEDot + cross(omega,RIB'*WE);
omegaDot = temp(4:6);
clear temp
xdot = [XDot; ThetaDot; VDot; omegaDot];
====================================script=================
clear
close all
% GenericFixedWingScript.m solves the nonlinear equations of motion for a 
% fixed-wing aircraft flying in ambient wind with turbulence. 
global e1 e2 e3 rho m g S b c AR Inertia 
% Basis vectors
e1 = [1;0;0];
e2 = [0;1;0];
e3 = [0;0;1];
% Atmospheric and gravity parameters (Constant altitude: Sea level)
%a = 340.3;          % Speed of sound (m/s)
rho = 1.225;        % Density (kg/m^3)
g = 9.80665;        % Gravitational acceleration (m/s^2)
%u0 = Ma*a;
u0 = 13.7;
P_dynamic = (1/2)*rho*u0^2;
% Aircraft parameters (Bix3)
m = 1.202;          % Mass (slugs)
W = m*g;            % Weight (Newtons)
Ix = 0.2163;         % Roll inertia (kg-m^2)
Iy = 0.1823;         % Pitch inertia (kg-m^2)
Iz = 0.3396;         % Yaw inertia (kg-m^2)
Ixz = 0.0364;            % Roll/Yaw product of inertia (kg-m^2)
%Inertia = [Ix, 0, -Ixz; 0, Iy, 0; -Ixz, 0, Iz];
S = 0.0285;         % Wing area (m^2)
b = 1.54;           % Wing span (m)
c = 0.188;          % Wing chord (m)
AR = (b^2)/S;       % Aspect ratio
CL_Trim = W/(P_dynamic*S);
CD_Trim = 0.036;
e = 0.6;    %Contrived
CD0 = CD_Trim - (CL_Trim^2)/(e*pi*AR);
% Equilibrium power (constant)
Thrust_Trim = P_dynamic*S*CD_Trim;
Power = Thrust_Trim*u0;
% Longitudinal nondimensional stability and control derivatives
%Cx0 = 0.197;
%Cxu = -0.156;
%Cxw = 0.297;
%Cxw2 = 0.960;
%Cz0 = -0.179;
%Czw = -5.32;
%Czw2 = 7.02;
%Czq = -8.20;
%Czde = -0.308;
%Cm0 = 0.0134;
%Cmw = -0.240;
%Cmq = -4.49;
%Cmde = -0.364;
% Lateral-directional nondimensional stability and control derivatives
X0 = ones(3,1);
Theta0 = ones(3,1);
V0 = u0*e1;
omega0 = ones(3,1);
y0 = [X0; Theta0; V0; omega0];
t_final = 10;
[t,y] = ode45('STOL_EOM',[0:0.1:t_final]',y0);
figure(1)
subplot(2,1,1)
plot(t,y(:,1:2))
ylabel('Position')
subplot(2,1,2)
plot(t,y(:,3:4))
ylabel('Attitude')
figure(2)
subplot(2,1,1)
plot(t,y(:,6:8))
ylabel('Velocity')
subplot(2,1,2)
plot(t,y(:,9:11))
ylabel('Angular Velocity')
Cesar Cardenas

Answer 1

[Harsh]

Hi @Cesar Cardenas,

To visualize how the elevator, aileron, and rudder deflections evolve over time during your simulation, you’ll need to extract or reconstruct these values at each time step and then plot them using MATLAB. Here's a high-level approach you can follow:

1. Integrate your dynamics using ode45

• Use the "ode45" solver to simulate your aircraft model over a time span.
• If control inputs are constant, define "de0", "da0", "dr0" before simulation and treat them as fixed throughout.
• If you plan to vary control inputs, you’ll need to modify your "STOL_EOM" to accept them as inputs (e.g., via interpolation tables or a controller function).

2. Record control inputs at each time step

• Since "de", "da", and "dr" are set as constants inside "STOL_EOM", you can reconstruct them after the simulation using their global values.

global de0 da0 dr0
de_array = de0 * ones(size(t));
da_array = da0 * ones(size(t));
dr_array = dr0 * ones(size(t));

3. Plot the control deflections

• Use "figure", "subplot", and "plot" to display each control signal over time:

figure
subplot(3,1,1)
  plot(t, de_array)
  ylabel('Elevator deflection \delta_e')
  grid on
subplot(3,1,2)
  plot(t, da_array)
  ylabel('Aileron deflection \delta_a')
  grid on
subplot(3,1,3)
  plot(t, dr_array)
  ylabel('Rudder deflection \delta_r')
  xlabel('Time (s)')
  grid on

Refer to the following documentation pages for more information:

• ode45 solver: https://www.mathworks.com/help/matlab/ref/ode45.html
• odeset options: https://www.mathworks.com/help/matlab/ref/odeset.html
• odeplot real-time plotting: https://www.mathworks.com/help/matlab/ref/odeplot.html
• Basic plotting (plot, figure): https://www.mathworks.com/help/matlab/ref/plot.html
• Subplots: https://www.mathworks.com/help/matlab/ref/subplot.html