how to generate plot

Hi, I trying to plot de, da, dr, and dT but I can not get any line or curve. I would like to know what could be my mistake or what I'm missing. any help will be greatly appreciated. Thanks

clear all
close all

% Lateral-directional nondimensional stability and control derivatives

Cm0 = 0.0134;
Cmw = -0.240;
Cmq = -4.49;
Cmde = -0.364;
thetaEq = 2*(pi/180);	% Initial guess for pitch (rad)
CT0 = 0.197;
rho = 1.225;
S = 0.285;
VEq = 12; % Initial guess for speed (m/s)
phi = 0.5;
Power = 200; % Maximum Power (W)
t = 10;

% Control deflections

deEq = -(Cm0 + Cmw*sin(thetaEq))/Cmde; % Initial guess for elevator (rad)
dTEq = (CT0*rho*S*(VEq^3))/(2*Power); % Nominal throttle setting



de = deEq;
da = 0.5;
%da = - 0.5*phi - 2*pr; % Add some roll stiffness and damping
dr = 0;
dT = dTEq;


plot(t,da)

[James Tursa]

You've got these two lines:

t = 10;
da = 0.5;

So yes, plot(t,da) isn't going to be very interesting.

What are the equations your are trying to plot? Do you have an image of those equations you can post?

[Cesar Cardenas]

@James Tursa In fact I have two files a function and a scrip and I'm trying to plot the control deflections by taking some code lines but I'm stuck with it. Any help will be appreciated. Thanks.

This is the function:

function xdot = GenericFixedWingEOM(t,x)

% GenericFixedWingEOM.m
%
% Nonlinear equations of motion for a rigid, fixed-wing airplane.


global e1 e2 e3 rho m g S b c Inertia Power ...
Cx0 Cxu Cxw Cxw2 ...
Cz0 Czw Czw2 Czq Czde ...
Cm0 Cmw Cmq Cmde ...
Cy0 Cyv Cyp Cyr Cyda Cydr ...
Cl0 Clv Clp Clr Clda Cldr ...
Cn0 Cnv Cnv2 Cnp Cnr Cnda Cndr ...
dTEq deEq ...
VonKarmanFlag Amp Phase Omega Phase2 Omega2 

% Parse out state vector components
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);    

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)];

% if VonKarmanFlag == 1
% 
% %%% WIND VELOCITY (1D VON KARMAN) %%%
% WX = Amp(:,1)'*cos(Omega*X(1)+Phase(:,1));
% WY = Amp(:,2)'*cos(Omega*X(1)+Phase(:,2));
% WZ = Amp(:,3)'*cos(Omega*X(1)+Phase(:,3));
% W = [WX; WY; WZ];
% 
% %%% WIND GRADIENT (1D VON KARMAN) %%%
% WXX = -(Amp(:,1).*Omega)'*sin(Omega*X(1)+Phase(:,1));
% WYX = -(Amp(:,2).*Omega)'*sin(Omega*X(1)+Phase(:,2));
% WZX = -(Amp(:,3).*Omega)'*sin(Omega*X(1)+Phase(:,3));
% GradW = [WXX, 0, 0; WYX, 0, 0; WZX, 0, 0];
% 
% elseif VonKarmanFlag == 2
%     
% %%% WIND VELOCITY (2D VON KARMAN) %%%
% WX = Amp(:,1)'*(cos(Omega*X(1)+Phase(:,1)).*cos(Omega2*X(2)+Phase2(:,1)));
% WY = Amp(:,2)'*(cos(Omega*X(1)+Phase(:,2)).*cos(Omega2*X(2)+Phase2(:,2)));
% WZ = Amp(:,3)'*(cos(Omega*X(1)+Phase(:,3)).*cos(Omega2*X(2)+Phase2(:,3)));;
% W = [WX; WY; WZ];
% 
% %%% WIND GRADIENT (2D VON KARMAN) %%%
% WXX = -(Amp(:,1).*Omega)'*(sin(Omega*X(1)+Phase(:,1)).*cos(Omega2*X(2)+Phase2(:,1)));
% WYX = -(Amp(:,2).*Omega)'*(sin(Omega*X(1)+Phase(:,1)).*cos(Omega2*X(2)+Phase2(:,1)));
% WZX = -(Amp(:,3).*Omega)'*(sin(Omega*X(1)+Phase(:,1)).*cos(Omega2*X(2)+Phase2(:,1)));
% WXY = -(Amp(:,1).*Omega2)'*(cos(Omega*X(1)+Phase(:,1)).*sin(Omega2*X(2)+Phase2(:,1)));
% WYY = -(Amp(:,2).*Omega2)'*(cos(Omega*X(1)+Phase(:,1)).*sin(Omega2*X(2)+Phase2(:,1)));
% WZY = -(Amp(:,3).*Omega2)'*(cos(Omega*X(1)+Phase(:,1)).*sin(Omega2*X(2)+Phase2(:,1)));
% GradW = [WXX, WXY, 0; WYX, WYY, 0; WZX, WZY, 0];
% 
% end

% Wind velocity in body frame, air-relative velocity, and dynamic pressure
% W = RIB'*W; % Convert wind velocity to body frame
W = [0;0;0]; 
Vr = V - W;
    ur = Vr(1);
    vr = Vr(2);
    wr = Vr(3); 
PDyn = (1/2)*rho*norm(Vr)^2;

% Wind gradient in body frame
%Phi = RIB'*GradW'*RIB;

% Effective angular velocity of the fluid (Note: Etkin argues for 
% a different formulation of the "angular velocity of the wind")
%temp = -(Phi-Phi');
%omegaf = [-temp(2,3); temp(1,3); -temp(1,2)]; 
omegaf =[0;0;0];
omegar = omega - omegaf;
    pr = omegar(1);
    qr = omegar(2);
    rr = omegar(3);

%%% AERODYNAMIC FORCES AND MOMENTS %%%

% Aerodynamic angles 
%beta = asin(Vr(2)/norm(Vr));
%alpha = atan2(Vr(3),Vr(1));

% Control deflections
de = deEq;
da = 0;
%da = - 0.5*phi - 2*pr; % Add some roll stiffness and damping
dr = 0;
dT = dTEq;

% Aerodynamic force coefficients
Cx = Cx0 + Cxu*(ur/norm(Vr)) + Cxw*(wr/norm(Vr)) + Cxw2*(wr/norm(Vr))^2; 
Cy = Cy0 + Cyv*(vr/norm(Vr)) + Cyp*((b*pr)/(2*norm(V))) + Cyr*((b*rr)/(2*norm(Vr))) ...
    + Cyda*da + Cydr*dr;
Cz = Cz0 + Czw*(wr/norm(Vr)) + Czw2*(wr/norm(Vr))^2 + Czq*((c*qr)/(2*norm(Vr))) ...
    + Czde*de;

Thrust = dT*Power/norm(Vr);  % Constant power
X = PDyn*S*Cx + Thrust;
Y = PDyn*S*Cy;
Z = PDyn*S*Cz;
Force_Aero = [X; Y; Z];

% Components of aerodynamic moment (modulo "unsteady" terms)
Cl = Cl0 + Clv*(vr/norm(Vr)) + Clp*((b*pr)/(2*norm(Vr))) + Clr*((b*rr)/(2*norm(Vr))) + ...
    Clda*da + Cldr*dr;
Cm = Cm0 + Cmw*(wr/norm(Vr)) + Cmq*((c*qr)/(2*norm(Vr))) + Cmde*de;
Cn = Cn0 + Cnv*(vr/norm(Vr)) + Cnv2*(vr/norm(Vr))^2 + Cnp*((b*pr)/(2*norm(Vr))) + Cnr*((b*rr)/(2*norm(Vr))) + ...
    Cnda*da + Cndr*dr;

L = PDyn*S*b*Cl;
M = PDyn*S*c*Cm;
N = PDyn*S*b*Cn;
Moment_Aero = [L; M; N];

% Sum of forces and moments
Force = RIB'*(m*g*e3) + Force_Aero;
Moment = Moment_Aero;

%%% EQUATIONS OF MOTION %%%

%%% KINEMATIC EQUATIONS %%%
XDot = RIB*V;
ThetaDot = LIB*omega;

%%% DYNAMIC EQUATIONS %%%
VDot = (1/m)*(cross(m*V,omega) + Force);
omegaDot = inv(Inertia)*(cross(Inertia*omega,omega) + Moment);

xdot = [XDot; ThetaDot; VDot; omegaDot];

and this is the 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 Inertia Power ...
Cx0 Cxu Cxw Cxw2 ...
Cz0 Czw Czw2 Czq Czde ...
Cm0 Cmw Cmq Cmde ...
Cy0 Cyv Cyp Cyr Cyda Cydr ...
Cl0 Clv Clp Clr Clda Cldr ...
Cn0 Cnv Cnv2 Cnp Cnr Cnda Cndr ...
dTEq deEq ...
VonKarmanFlag Amp Phase Omega Phase2 Omega2 

% 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)

% Aircraft parameters (Bix3)
m = 1.202;          % Mass (kg)
W = m*g;            % Weight (Newtons)
Ix = 0.095;         % Roll inertia (kg-m^2)
Iy = 215e3;         % Pitch inertia (kg-m^2)
Iz = 447e3;         % Yaw inertia (kg-m^2)
Ixz = 0;            % Roll/Yaw product of inertia (kg-m^2)
Inertia = [Ix, 0, -Ixz; 0, Iy, 0; -Ixz, 0, Iz];

S = 0.285;         % Wing area (m^2)
b = 1.54;           % Wing span (m)
c = 0.188;          % Wing chord (m)
AR = (b^2)/S;       % Aspect ratio

% Longitudinal nondimensional stability and control derivatives
Cx0 = 0; % Define CT0 = 0.197 below. (Cx0 = 0.197 in [Simmons et al, 2019])
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
Cy0 = 0;
Cyv = -0.251;
Cyp = 0.170;
Cyr = 0.350;
Cyda = 0.103;
Cydr = 0.0157;

Cl0 = 0;
Clv = -0.0756;
Clp = -0.319;
Clr = 0.183;
Clda = 0.170;
Cldr = -0.0117;

Cn0 = 0;
Cnv = 0.0408;
Cnv2 = 0.0126;
Cnp = -0.242;
Cnr = -0.166;
Cnda = 0.0416;
Cndr = -0.0618;

% Solve for wings-level, constant-altitude equilibrium flight condition
VEq = 12;               % Initial guess for speed (m/s)
thetaEq = 2*(pi/180);	% Initial guess for pitch (rad)
deEq = -(Cm0 + Cmw*sin(thetaEq))/Cmde; % Initial guess for elevator (rad)
Power = 200;            % Maximum Power (W)
CT0 = 0.197;            % CT0 = Cx0 in [Simmons et al, 2019]
dTEq = (CT0*rho*S*(VEq^3))/(2*Power); % Nominal throttle setting

%Roots = fsolve('WingsLevelEquilibrium',[VEq;thetaEq;deEq]);
%VEq = Roots(1);
%thetaEq = Roots(2);
%deEq = Roots(3);

% Generate amplitude, phase, and frequency parameters for 1D or 2D von Karman turbulence
%VonKarmanFlag = 1

% 1D von Karman turbulence
%if VonKarmanFlag == 1

%[Amp,Phase,Omega] = OneDVonKarmanTurbulence(-4,0,100);
%Omega = Omega/0.3048; % (rad/m)
%disp('One Dimensional Turbulence Spectrum')

% 2D von Karman turbulence
%elseif VonKarmanFlag == 2

%[Amp,Phase,Phase2,Omega,Omega2] = TwoDVonKarmanTurbulence(-4,0,100);
%Omega = Omega/0.3048; % (rad/m)
%Omega2 = Omega2/0.3048; % (rad/m)
%disp('Two Dimensional Turbulence Spectrum')

%end

%Amp = Amp*0.3048; % (m/s)
% Amp = 0*Amp; % Zero out the turbulence.

% Define initial state
X0 = zeros(3,1); % (m)
Theta0 = [0;thetaEq;0]; % (rad)
V0 = [VEq*cos(thetaEq);0;VEq*sin(thetaEq)]; % (m/s)
omega0 = zeros(3,1); % (rad/s)
%deEq = zeros(3,1);

y0 = [X0; Theta0; V0; omega0];
t_final = 10;
[t,y] = ode45('GenericFixedWingEOM',[0:0.1:t_final]',y0);

%%%% CODE TO COMPUTE FLOW RELATIVE VELOCITY, ETC %%%%

figure(1)
subplot(2,1,1)
plot(t,y(:,1:3))
legend('X','Y','Z')
ylabel('Position (m)')
subplot(2,1,2)
plot(t,y(:,4:6)*(180/pi))
legend('\phi','\theta','\psi')
ylabel('Attitude (deg)')

figure(2)
subplot(2,1,1)
plot(t,y(:,7:9))
legend('u','v','w')
ylabel('Velocity (m/s)')
subplot(2,1,2)
plot(t,y(:,10:12)*(180/pi))
ylabel('Angular Velocity (deg/s)')
legend('p','q','r')