how to obtain plots from a given function
8 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Torsten
am 17 Sep. 2022
Plotting should be done when the ODE integrator returns a solution, not from the ODE function used to get the derivatives.
Antworten (0)
Siehe auch
Kategorien
Mehr zu Environmental Models finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!