Utilizing the Euler Method

11 views (last 30 days)
Turgut Ataseven
Turgut Ataseven on 1 Jan 2022
Commented: Turgut Ataseven on 1 Jan 2022
Hi. First of all sorry for asking a huge question.I need a guide to utilize the Euler Method. I am having troubles with the sizes of variables.
This code tries to solve 6 ODEs with 6 state variables [horizontal position (x1 and x2), altitude (x3), the true airspeed (x4), the heading angle (x5) and the mass of the aircraft (x6)] and 3 control inputs [engine thrust (u1), the bank angle (u2) and the flight path angle (u3)] by using Euler's method. Different flight maneuvers are performed for the specified time intervals.
Velocities.m, Cruise_Vel.m, Des_Vel.m, Thr_cl.m, Thr_cr.m, Thr_des.m, fuel_cl.m, fuel_cr.m, fuel_des.m,den.m,drag.m,drag_cr.m,lift.m,lift_cr.m are functions in seperate tabs. Code in the main file (Flight.m) is:
% Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
% Perform cruise flight for t=60 minutes.
% Turn with β=30 bank angle until heading is changed by η=270◦.
% Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
% Complete a 360◦ turn (loiter) at level flight.
% Descent to h3=800 [m] with κ=4.5◦ flight path angle.
% Aircraft Properties
W = .44225E+06; % .44225E+03 tons = .44225E+06 kg
S = .51097E+03; % Surface Area [m^2]
g0 = 9.80665; % Gravitational acceleration [m/s2]
% solving 1st order ODE using numerical methods
t0=0;
tend=3960;
h=0.05;
N=(tend-t0)/h;
t=t0:h:tend;
x = zeros(6,length(t));
% Initial conditions
x(:,1)=[0;0;3608.92;1.0e+02 * 1.161544478045788;0;W];
for i=2:length(t)
if and (t(i-1) > 0,t(i-1)<=60) % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
x3 = 3608.9:5249.3; % Changing altitude [m] -> [ft]
x4 = Velocities(3608.9,5249.3); % Changing speed [m/s]
x5 = 0; % Changing head angle [deg]
f = fuel_cl(3608.9,5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cl(3608.9,5249.3); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 5; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag(3608.9,5249.3,x4); % Changing drag coefficient
Cl = lift(3608.9,5249.3,x4); % Changing lift coefficient
p = den(3608.9,5249.3); % Changing density [kg/m3]
elseif and (t(i-1) >60,t(i-1)<=3660) % Perform cruise flight for t=60 minutes.
x3 = 5249.3; % Changing altitude [m] -> [ft]
x4 = Cruise_Vel(5249.3); % Changing speed [m/s]
x5 = 0; % Changing head angle [deg]
f = fuel_cr(5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cr(5249.3); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(5249.3,5249.3,x4); % Changing drag coefficient
Cl = lift_cr(5249.3,5249.3,x4); % Changing lift coefficient
p = den(5249.3,5249.3); % Changing density [kg/m3]
elseif and (t (i-1)>3660,t(i-1)<=3720) % Turn with β=30 bank angle until heading is changed by η=270◦.
x3 = 5249.3; % Changing altitude [m] -> [ft]
x4 = Cruise_Vel(5249.3); % Changing speed [m/s]
x5 = 0:30:270; % Changing head angle [deg]
f = fuel_cr(5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cr(5249.3); % Changing thrust [N]
u2 = 30; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(5249.3,5249.3,x4); % Changing drag coefficient
Cl = lift_cr(5249.3,5249.3,x4); % Changing lift coefficient
p = den(5249.3,5249.3); % Changing density [kg/m3]
elseif and (t(i-1) >3720,t(i-1)<=3780) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
x3 = 5249.3:3608.9; % Changing altitude [m] -> [ft]
x4 = Des_Vel(5249.3,3608.9); % Changing speed [m/s]
x5 = 270; % Changing head angle [deg]
f = fuel_des(5249.3,3608.9); % Changing fuel flow [kg/min]
u1 = Thr_des(5249.3,3608.9); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 4; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag(5249.3,3608.9,x4); % Changing drag coefficient
Cl = lift(5249.3,3608.9,x4); % Changing lift coefficient
p = den(5249.3,3608.9); % Changing density [kg/m3]
elseif and (t (i-1)>3780,t(i-1)<=3900) % Complete a 360◦ turn (loiter) at level flight.
x3 = 3608.9; % Changing altitude [m] -> [ft]
x4 = Cruise_Vel(3608.9); % Changing speed [m/s]
lon = [270 300 360 60 120 180 240 270];
x5 = wrapTo360(lon); % Changing head angle [deg]
f = fuel_cr(3608.9); % Changing fuel flow [kg/min]
u1 = Thr_cr(3608.9); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(3608.9,3608.9,x4); % Changing drag coefficient
Cl = lift_cr(3608.9,3608.9,x4); % Changing lift coefficient
p = den(3608.9,3608.9); % Changing density [kg/m3]
elseif and (t(i-1) >3900,t(i-1)<=3960) % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
x3 = 3608.9:2624.67; % Changing altitude [m] -> [ft]
x4 = Des_Vel(3608.9,2624.67); % Changing speed [m/s]
x5 = 270; % Changing head angle [deg]
f = fuel_des(3608.9,2624.67); % Changing fuel flow [kg/min]
u1 = Thr_des(3608.9,2624.67); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 4.5; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag(3608.9,2624.67,x4); % Changing drag coefficient
Cl = lift(3608.9,2624.67,x4); % Changing lift coefficient
p = den(3608.9,2624.67); % Changing density [kg/m3]
end
dx1dt = x4 .* cos(x5) .* cos(u3);
dx2dt = x4 .* sin(x5) .* cos(u3);
dx3dt = x4 .* sin(u3);
dx4dt = -C_D.*S.*p.*(x4.^2)./(2.*x6)-g0.*sin(u3)+u1./x6;
dx5dt = -Cl.*S.*p.*x4./(2.*x6).*sin(u2);
dx6dt = -f;
x(1,i)= x(1,i-1) + h * dx1dt; % line 129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x(2,i)= x(2,i-1) + h * dx2dt;
x(3,i)= x(3,i-1) + h * dx3dt;
x(4,i)= x(4,i-1) + h * dx4dt;
x(5,i)= x(5,i-1) + h * dx5dt;
x(6,i)= x(6,i-1) + h * dx6dt;
end
tot=cell2mat(f); % Total fuel consumption during mission [kg/min]
Tot_fuel=sum(tot);
figure(1)
plot3(x1(:),x2(:),x3(:)); % 3D position graph
figure(2)
plot(t,x4(:)); % Vtas − Time graph
figure(3)
plot(t,V_ver(:)); % V_vertical − Time graph
figure(4)
plot(t,x5(:)); % Heading − Time graph
figure(5)
plot(t,x6(:)); % Mass − Time graph
figure(6)
plot(t,u1(:)); % Thrust − Time graph
figure(7)
plot(t,u2(:)); % Bank Angle − Time graph
figure(8)
plot(t,u3(:)); % Flight Path Angle − Time graph
fprintf('Total fuel consumption during mission is %.2f [kg]',Tot_fuel*tend/60);
Aforementioned functions in seperate tabs gets altitude interval as input and outputs the answer. Here is some of them:
function [Vtas_cl] = Velocities(H1,H2)
Vcl_1 = 335; % Standard calibrated airspeed [kt]
Vcl_2 = 172.3; % Standard calibrated airspeed [kt] -> [m/s] (To find the Mach transition altitude)
Vcl_2_in_knots = 335; % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cl = 0.86; % Standard calibrated airspeed [kt]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
g0 = 9.80665; % Gravitational acceleration [m/s2]
a0= 340.294; % Speed of Sound [m/s]
Vd_CL1 = 5; % Climb speed increment below 1500 ft (jet)
Vd_CL2 = 10; % Climb speed increment below 3000 ft (jet)
Vd_CL3 = 30; % Climb speed increment below 4000 ft (jet)
Vd_CL4 = 60; % Climb speed increment below 5000 ft (jet)
Vd_CL5 = 80; % Climb speed increment below 6000 ft (jet)
CV_min = 1.3; % Minimum speed coefficient
Vstall_TO = .14200E+03; % Stall speed at take-off [KCAS]
CAS_climb = Vcl_2;
Mach_climb = M_cl;
delta_trans = (((1+((K-1)/2)*(CAS_climb/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_climb^2)^(K/(K-1)))-1); % Pressure ratio at the transition altitude
teta_trans = delta_trans ^ (-Bt*R/g0); % Temperature ratio at the transition altitude
H_p_trans_climb = (1000/0.348/6.5)*(T0*(1-teta_trans)); % Transition altitude for climb [ft]
H_climb = H1:H2;
for k1 = 1:length(H_climb)
if (0<=H_climb(k1)&&H_climb(k1)<1500)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL1;
elseif (1500<=H_climb(k1)&&H_climb(k1)<3000)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL2;
elseif (3000<=H_climb(k1)&&H_climb(k1)<4000)
Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL3;
elseif (4000<=H_climb(k1)&&H_climb(k1)<5000)
Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL4;
elseif (5000<=H_climb(k1)&&H_climb(k1)<6000)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL5;
elseif (6000<=H_climb(k1)&&H_climb(k1)<10000)
Vnom_climb_jet (k1)= min(Vcl_1,250);
elseif (10000<=H_climb(k1)&&H_climb(k1)<=H_p_trans_climb)
Vnom_climb_jet(k1) = Vcl_2_in_knots;
elseif (H_p_trans_climb<H_climb(k1))
Vnom_climb_jet(k1) = M_cl;
end
Vcas_cl(k1) = Vnom_climb_jet(k1)* 0.514; % [kn] -> [m/s]
H_climb (k1)= H_climb(k1) * 0.3048; % [feet] -> [m]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
deltaT = 0; % Value of the real temperature T in ISA conditions [K]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
P0 = 101325; % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665; % Gravitational acceleration [m/s2]
p0 = 1.225; % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;
T(k1) = T0 + deltaT + Bt * H_climb(k1); % Temperature [K]
P (k1)= P0*((T(k1)-deltaT)/T0).^((-g0)/(Bt*R)); % Pressure [Pa]
p (k1)= P(k1) ./ (R*T(k1)); % Density [kg/m^3]
Vtas_cl(k1) = (2*P(k1)/visc/p(k1)*((1 + P0/P(k1)*((1 + visc*p0*Vcas_cl(k1)*Vcas_cl(k1)/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
end
end
Another one:
function [Thr_jet_climb_ISA] = Thr_cl(H1,H2)
H_climb = H1:H2; % Climb altitude range [m] -> [ft]
for k3 = 1:length(H_climb)
C_Tc_1 = .75979E+06;
C_Tc_2 = .52423E+05;
C_Tc_3 = .40968E-10;
Thr_jet_climb_ISA(k3) = C_Tc_1* (1-(H_climb(k3)/C_Tc_2)+C_Tc_3*H_climb(k3).^2); % Maximum climb and take-off thrust for jet engine [N]
end
end
And the ones with the single altitude input:
function [Vtas_cr] = Cruise_Vel(H_cruise)
Vcr_1 = 320; % Standard calibrated airspeed [kt]
Vcr_2 = 164.62; % Standard calibrated airspeed [kt] -> [m/s] (To find the Mach transition altitude)
Vcr_2_in_knots = 320; % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cr = 0.84; % Standard calibrated airspeed [kt]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
g0 = 9.80665; % Gravitational acceleration [m/s2]
a0= 340.294; % Speed of Sound [m/s]
deltaT = 0; % Value of the real temperature T in ISA conditions [K]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
P0 = 101325; % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665; % Gravitational acceleration [m/s2]
p0 = 1.225; % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;
CAS_cruise = Vcr_2;
Mach_cruise = M_cr;
delta_trans = (((1+((K-1)/2)*(CAS_cruise/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_cruise^2)^(K/(K-1)))-1); % Pressure ratio at the transition altitude
teta_trans = delta_trans ^ (-Bt*R/g0); % Temperature ratio at the transition altitude
H_p_trans_cruise = (1000/0.348/6.5)*(T0*(1-teta_trans)); % Transition altitude for cruise [ft]
if (0<=H_cruise&&H_cruise<=2999)
Vnom_cruise_jet = min(Vcr_1,170);
elseif (3000<=H_cruise&&H_cruise<=5999)
Vnom_cruise_jet = min(Vcr_1,220);
elseif (6000<=H_cruise&&H_cruise<=13999)
Vnom_cruise_jet = min(Vcr_1,250);
elseif (14000<=H_cruise&&H_cruise<=H_p_trans_cruise)
Vnom_cruise_jet = Vcr_2_in_knots;
elseif (H_p_trans_cruise < H_cruise)
Vnom_cruise_jet = M_cr;
end
Vcas_cr = Vnom_cruise_jet * 0.514; % [knots] -> [m/s]
H_cruise = H_cruise * 0.3048; % [feet] -> [m]
T = T0 + deltaT + Bt * H_cruise; % Temperature [K]
P = P0*((T-deltaT)/T0).^((-g0)/(Bt*R)); % Pressure [Pa]
p = P ./ (R*T); % Density [kg/m^3]
Vtas_cr = (2*P/visc/p*((1 + P0/P*((1 + visc*p0*Vcas_cr*Vcas_cr/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
end
Another single-input function:
function [Thr_cra] = Thr_cr(H_cruise)
W = .44225E+06; % .44225E+03 tons = .44225E+06 kg
S = .51097E+03; % Surface Area [m^2]
Vcr_1 = 320; % Standard calibrated airspeed [kt]
Vcr_2 = 164.62; % Standard calibrated airspeed [kt] -> [m/s] (To find the Mach transition altitude)
Vcr_2_in_knots = 320; % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cr = 0.84; % Standard calibrated airspeed [kt]
K = 1.4; % Adiabatic index of air
R = 287.05287; % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065; % ISA temperature gradient with altitude below the tropopause [K/m]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
g0 = 9.80665; % Gravitational acceleration [m/s2]
a0= 340.294; % Speed of Sound [m/s]
deltaT = 0; % Value of the real temperature T in ISA conditions [K]
T0 = 288.15; % Standard atmospheric temperature at MSL [K]
P0 = 101325; % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665; % Gravitational acceleration [m/s2]
p0 = 1.225; % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;
CAS_cruise = Vcr_2;
Mach_cruise = M_cr;
delta_trans = (((1+((K-1)/2)*(CAS_cruise/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_cruise^2)^(K/(K-1)))-1); % Pressure ratio at the transition altitude
teta_trans = delta_trans ^ (-Bt*R/g0); % Temperature ratio at the transition altitude
H_p_trans_cruise = (1000/0.348/6.5)*(T0*(1-teta_trans)); % Transition altitude for cruise [ft]
if (0<=H_cruise&&H_cruise<=2999)
Vnom_cruise_jet = min(Vcr_1,170);
elseif (3000<=H_cruise&&H_cruise<=5999)
Vnom_cruise_jet = min(Vcr_1,220);
elseif (6000<=H_cruise&&H_cruise<=13999)
Vnom_cruise_jet = min(Vcr_1,250);
elseif (14000<=H_cruise&&H_cruise<=H_p_trans_cruise)
Vnom_cruise_jet = Vcr_2_in_knots;
elseif (H_p_trans_cruise < H_cruise)
Vnom_cruise_jet = M_cr;
end
Vcas_cr = Vnom_cruise_jet * 0.514; % [knots] -> [m/s]
H_cruise = H_cruise * 0.3048; % [feet] -> [m]
T = T0 + deltaT + Bt * H_cruise; % Temperature [K]
P = P0*((T-deltaT)/T0).^((-g0)/(Bt*R)); % Pressure [Pa]
p = P ./ (R*T); % Density [kg/m^3]
Vtas_cr = (2*P/visc/p*((1 + P0/P*((1 + visc*p0*Vcas_cr*Vcas_cr/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
Cl = 2*W / (p*Vtas_cr*Vtas_cr*S); % Lift coefficient
C_D0cr = .25669E-01;
C_D2cr = .39082E-01;
C_D = C_D0cr + C_D2cr * (Cl).^2; % Drag Coefficient
D = 1/2 * p * Vtas_cr* Vtas_cr * S * C_D; % Drag Force [N]
Thr_cra = D; % Cruise drag (Thrust) [N]
end
The rest is like the above ones, they get either one or two altitude inputs.
When I run the code, this error message displays:
>> Flight
Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-1641.
Error in Flight (line 129)
x(1,i)= x(1,i-1) + h * dx1dt;
The reason why it says right side is 1-by-1641 is because in the first if-else, x4 is defined as:
x4 = Velocities(3608.9,5249.3); % Changing speed [m/s]
A guide on what to do would be really helpful to me. If you want to analyse the code: https://drive.google.com/drive/folders/1RtM_XMuhQs2voZFXOnoFAgtkxMI3nKFt?usp=sharing
Thanks.
  1 Comment
Stephen23
Stephen23 on 1 Jan 2022
@Turgut Ataseven: so much code makes your question very hard to follow and understand.
Rather than showing all code in the question, consider uploading it by clicking the paperclip button.

Sign in to comment.

Accepted Answer

Torsten
Torsten on 1 Jan 2022
Edited: Torsten on 1 Jan 2022
In the if-statement in Flight.m, set x3 to x(3,i-1) and calculate the other variables x4,x5,f,u1,u2,u3,V_ver,C_D,Cl and p only for this single value, not for the complete range (mostly 3608.9:5249.3).
  1 Comment
Turgut Ataseven
Turgut Ataseven on 1 Jan 2022
Nice approach. Here is what I did.
x = zeros(6,length(t));
x1 = zeros(1,length(t));
x2 = zeros(1,length(t));
x3 = zeros(1,length(t));
x4 = zeros(1,length(t));
x5 = zeros(1,length(t));
x6 = zeros(1,length(t));
u1 = zeros(1,length(t));
u2 = zeros(1,length(t));
u3 = zeros(1,length(t));
C_D= zeros(1,length(t));
p = zeros(1,length(t));
Cl = zeros(1,length(t));
f = zeros(1,length(t));
% Initial conditions
x(:,1)=[0;0;3608.92;1.0e+02 * 1.161544478045788;0;W];
for i=2:length(t)
if and (x(3,i-1) > 3608.92,x(3,i-1)<=5249.3) % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
t =0:60;
x4 = Velocities(3608.9,5249.3); % Changing speed [m/s]
x5 = 0; % Changing head angle [deg]
f = fuel_cl(3608.9,5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cl(3608.9,5249.3); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 5; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag(3608.9,5249.3,x4); % Changing drag coefficient
Cl = lift(3608.9,5249.3,x4); % Changing lift coefficient
p = den(3608.9,5249.3); % Changing density [kg/m3]
elseif and (x(3,i-1) > 5249.3,x(3,i-1)<=5249.3) % Perform cruise flight for t=60 minutes.
for t=60:3660
x4 = Cruise_Vel(5249.3); % Changing speed [m/s]
x5 = 0; % Changing head angle [deg]
f = fuel_cr(5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cr(5249.3); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(5249.3,5249.3,x4); % Changing drag coefficient
Cl = lift_cr(5249.3,5249.3,x4); % Changing lift coefficient
p = den(5249.3,5249.3); % Changing density [kg/m3]
end
% Turn with β=30 bank angle until heading is changed by η=270◦.
for t=3660:3720
x4 = Cruise_Vel(5249.3); % Changing speed [m/s]
x5 = 0:30:270; % Changing head angle [deg]
f = fuel_cr(5249.3); % Changing fuel flow [kg/min]
u1 = Thr_cr(5249.3); % Changing thrust [N]
u2 = 30; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(5249.3,5249.3,x4); % Changing drag coefficient
Cl = lift_cr(5249.3,5249.3,x4); % Changing lift coefficient
p = den(5249.3,5249.3); % Changing density [kg/m3]
end
elseif and (x(3,i-1) < 5249.3,x(3,i-1)>=3608.9) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
t=3720:3780;
x4 = Des_Vel(5249.3,3608.9); % Changing speed [m/s]
x5 = 270; % Changing head angle [deg]
f = fuel_des(5249.3,3608.9); % Changing fuel flow [kg/min]
u1 = Thr_des(5249.3,3608.9); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 4; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_des(5249.3,3608.9,x4); % Changing drag coefficient
Cl = lift_des(5249.3,3608.9,x4); % Changing lift coefficient
p = den_des(5249.3,3608.9); % Changing density [kg/m3]
elseif and (x(3,i-1) < 3608.9,x(3,i-1)>=3608.9) % Complete a 360◦ turn (loiter) at level flight.
t = 3780:3900;
x4 = Cruise_Vel(3608.9); % Changing speed [m/s]
lon = [270 300 360 60 120 180 240 270];
x5 = wrapTo360(lon); % Changing head angle [deg]
f = fuel_cr(3608.9); % Changing fuel flow [kg/min]
u1 = Thr_cr(3608.9); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 0; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_cr(3608.9,3608.9,x4); % Changing drag coefficient
Cl = lift_cr(3608.9,3608.9,x4); % Changing lift coefficient
p = den(3608.9,3608.9); % Changing density [kg/m3]
elseif and (x(3,i-1) < 3608.9,x(3,i-1)>=2624.67) % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
t=3900:3960;
x4 = Des_Vel(3608.9,2624.67); % Changing speed [m/s]
x5 = 270; % Changing head angle [deg]
f = fuel_des(3608.9,2624.67); % Changing fuel flow [kg/min]
u1 = Thr_des(3608.9,2624.67); % Changing thrust [N]
u2 = 0; % Changing bank angle [deg]
u3 = 4.5; % Changing flight path angle [deg]
V_ver = x4*sin(u3); % Changing vertical speed [m/s]
C_D = drag_des(3608.9,2624.67,x4); % Changing drag coefficient
Cl = lift_des(3608.9,2624.67,x4); % Changing lift coefficient
p = den_des(3608.9,2624.67); % Changing density [kg/m3]
else
fprintf("A problem occured.");
end
dx1dt = x4 .* cos(x5) .* cos(u3);
dx2dt = x4 .* sin(x5) .* cos(u3);
dx3dt = x4 .* sin(u3);
dx4dt = -C_D.*S.*p.*(x4.^2)./(2.*x6)-g0.*sin(u3)+u1./x6;
dx5dt = -Cl.*S.*p.*x4./(2.*x6).*sin(u2);
dx6dt = -f;
x(1,i)= x(1,i-1) + h * dx1dt;
x(2,i)= x(2,i-1) + h * dx2dt;
x(3,i)= x(3,i-1) + h * dx3dt;
x(4,i)= x(4,i-1) + h * dx4dt;
x(5,i)= x(5,i-1) + h * dx5dt;
x(6,i)= x(6,i-1) + h * dx6dt;
end
With these statements I've tried to state that altitude is constant. Is there a better way to do it?
elseif and (x(3,i-1) > 5249.3,x(3,i-1)<=5249.3) % Perform cruise flight for t=60 minutes.
and
elseif and (x(3,i-1) < 3608.9,x(3,i-1)>=3608.9) % Complete a 360◦ turn (loiter) at level flight.

Sign in to comment.

More Answers (1)

Burhan Burak AKMAN
Burhan Burak AKMAN on 1 Jan 2022
Change your code like this. May be work.
x(1,i)= x(1,i-1) + h * dx1dt(1,i-1); % line 129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x(2,i)= x(2,i-1) + h * dx2dt(1,i-1);
x(3,i)= x(3,i-1) + h * dx3dt(1,i-1);
x(4,i)= x(4,i-1) + h * dx4dt(1,i-1);
x(5,i)= x(5,i-1) + h * dx5dt(1,i-1);
x(6,i)= x(6,i-1) + h * dx6dt(1,i-1);

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!