Clear Filters
Clear Filters

how to obtain plots from a given function

2 views (last 30 days)
Hi, not sure how to plot xd1, xd2,,,,,xd3,,,any help will be greatly appreciated. Thanks
function xdot = EoMv2(t,x)
global mUser Ixx Iyy Izz Ixz S b cBar SMI CONHIS u tuHis deluHis ...
ACtype MODEL TrimHist RUNNING
D2R = pi/180;
R2D = 180/pi;
% Terminate if Altitude becomes Negative
[value,isterminal,direction] = event(t,x);
% Earth-to-Body-Axis Transformation Matrix
HEB = DCM(x(10),x(11),x(12));
% Atmospheric State
x(6) = min(x(6),0); % Limit x(6) <= 0 m
[airDens,airPres,temp,soundSpeed] = Atmos(-x(6));
% Body-Axis Wind Field
windb = WindField(-x(6),x(10),x(11),x(12));
% Body-Axis Gravity Components
gb = HEB * [0;0;9.80665];
% Air-Relative Velocity Vector
x(1) = max(x(1),0); % Limit axial velocity to >= 0 m/s
Va = [x(1);x(2);x(3)] + windb;
V = sqrt(Va'*Va);
alphar = atan(Va(3)/abs(Va(1)));
alpha = R2D * alphar;
betar = asin(Va(2)/V);
beta = R2D*betar;
Mach = V/soundSpeed;
qbar = 0.5*airDens*V^2;
% Incremental Flight Control Effects
if CONHIS >=1 && RUNNING == 1
[uInc] = interp1(tuHis,deluHis,t);
uInc = (uInc)';
uTotal = u + uInc;
else
uTotal = u;
end
% Force and Moment Coefficients
if MODEL == 'Alph'
[CD,CL,CY,Cl,Cm,Cn,mRef,S,Ixx,Iyy,Izz,Ixz,cBar,b,Thrust] = AlphaModel(x,uTotal,alphar,betar,V);
end
% if MODEL == 'Mach'
% [CD,CL,CY,Cl,Cm,Cn,mRef,S,Ixx,Iyy,Izz,Ixz,cBar,b,Thrust] = MachModel(x,uTotal,Mach,alphar,betar,V);
% end
qbarS = qbar*S;
CX = -CD*cos(alphar) + CL*sin(alphar); % Body-axis X coefficient
CZ = -CD*sin(alphar) - CL*cos(alphar); % Body-axis Z coefficient
% State Accelerations
Xb = (CX*qbarS + Thrust)/mRef;
Yb = CY*qbarS/mRef;
Zb = CZ*qbarS/mRef;
Lb = Cl*qbarS*b;
Mb = Cm*qbarS*cBar;
Nb = Cn*qbarS*b;
nz = -Zb/9.80665; % Normal load factor
% Dynamic Equations
xd1 = Xb + gb(1) + x(9)*x(2) - x(8)*x(3);
xd2 = Yb + gb(2) - x(9)*x(1) + x(7)*x(3);
xd3 = Zb + gb(3) + x(8)*x(1) - x(7)*x(2);
y = HEB' * [x(1);x(2);x(3)];
xd4 = y(1);
xd5 = y(2);
xd6 = y(3);
xd7 = (Izz * Lb + Ixz*Nb - (Ixz*(Iyy - Ixx - Izz)*x(7) + ...
(Ixz^2 + Izz*(Izz - Iyy))*x(9)) * x(8))/(Ixx * Izz - Ixz^2);
xd8 = (Mb - (Ixx - Izz)*x(7)*x(9) - Ixz*(x(7)^2 - x(9)^2))/Iyy;
xd9 = (Ixz*Lb + Ixx*Nb + (Ixz*(Iyy - Ixx - Izz)*x(9) + ...
(Ixz^2 + Ixx*(Ixx - Iyy))*x(7))*x(8))/(Ixx*Izz - Ixz^2);
cosPitch = cos(x(11));
if abs(cosPitch) <= 0.00001
cosPitch = 0.00001 * sign(cosPitch);
end
tanPitch = sin(x(11)) / cosPitch;
xd10 = x(7) + (sin(x(10))*x(8) + cos(x(10))*x(9))*tanPitch;
xd11 = cos(x(10))*x(8) - sin(x(10))*x(9);
xd12 = (sin(x(10))*x(8) + cos(x(10))*x(9))/cosPitch;
xdot = [xd1;xd2;xd3;xd4;xd5;xd6;xd7;xd8;xd9;xd10;xd11;xd12];
end
  2 Comments
Torsten
Torsten on 17 Sep 2022
Plotting should be done when the ODE integrator returns a solution, not from the ODE function used to get the derivatives.
Cesar Cardenas
Cesar Cardenas on 19 Sep 2022
Thanks for your reply. Yes, thanks I solved it. Now, I need to apply a square doublet with 1 second of duration. I do not really know how to do it. Also, I need to incorporate this dynamic model:
this is my attempt but not sure,? Any help will be greatly appreciated. Thanks
%Dynamic Actuator Model
d = (e^-(t*s)/t1*s + 1)*dc;
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')

Sign in to comment.

Answers (0)

Categories

Find more on Environmental Models 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!