how to obtain plots from a given function

2 次查看（过去 30 天）
Cesar Cardenas
Cesar Cardenas 2022-9-17
评论： Cesar Cardenas 2022-9-19
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 个评论
Torsten
Torsten 2022-9-17
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 2022-9-19
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')

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Aerodynamics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by