# Utilizing the Euler Method

2 views (last 30 days)
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.
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.

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

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

### Categories

Find more on Roads, Lanes, and Markings in Help Center and File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!