Clear Filters
Clear Filters

how to include variables for plotting

4 views (last 30 days)
Cesar Cardenas
Cesar Cardenas on 19 Sep 2022
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

Answers (0)

Categories

Find more on Graphics in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!