Index exceeds the number of array elements (1).

Hi. My 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.

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,drag.m,lift.m,den.m are functions in seperate files. Code in the main file (FlightPlan.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]

% Equations

% dx1 = x4*cos(x5/180*pi)*cos(u3/180*pi);       
% dx2 = x4*sin(x5/180*pi)*cos(u3/180*pi);       
% dx3 = x4*sin(u3/180*pi);               
% dx4 = -C_D*S*p*x4.*x4/(2*x6)-g0*sin(u3/180*pi)+u1/x6; 
% dx5 = -Cl*S*p*x4/(2*x6)/sin(u2/180*pi);
% dx6 = -f;       

% Initial conditions

x1(1) = 0;               % Initial position [m]
x2(1) = 0;               % Initial position [m]
x3(1) = 3608.92;         % Initial altitude [ft]
x4(1) = 1.0e+02 * 1.161544478045788;   % Initial speed [m/s]
x5(1) = 0;               % Assuming aircraft headed to North initially.
x6(1) = W;               % Initial mass [kg]

u1(1) = 1.0e+05 * 7.078900012026454;   % Initial thrust [N]
u2(1) = 0;               % Initial bank angle [deg]
u3(1) = 5;               % Initial flight path angle [deg]

C_D(1) = 0.026278212364444;            % Initial drag coefficient
Cl(1) = 0.124852132431623;             % Initial lift coefficient
p(1) = 1.027625283894403;              % Initial density [kg/m3]
f(1) = 1.0e+02 * 4.497203403070063;    % Initial fuel flow [kg/min]

dx1dt = @(x4,x5,u3)            x4.*cos(x5*pi/180).*cos(u3*pi/180);       
dx2dt = @(x4,x5,u3)            x4*sin(x5*pi/180)*cos(u3*pi/180);       
dx3dt = @(x4,u3)               x4*sin(u3*pi/180);               
dx4dt = @(C_D,p,x6,x4,u3,u1)   -C_D*S*p*x4.*x4/(2*x6)-g0*sin(u3)+u1/x6; 
dx5dt = @(Cl,p,x6,x4,u2)       -Cl*S*p*x4/(2*x6)/sin(u2);
dx6dt = @(f)                   -f;       


% solving 1st order ODE using numerical methods 

t0=0;           
tend=66;
y0=0;
h=0.2; 
N=(tend-t0)/h;
t=t0:h:tend;

if and (t >= 0,t<=1) % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
   x3(t) = linspace(3608.9,5249.3,2);                   % Changing altitude [m] -> [ft]
   x4(t) = Velocities(3608.9,5249.3);                   % Changing speed [m/s]
   x5(t) = 0;                                           % Changing head angle [deg]
   f(t) = fuel_cl(3608.9,5249.3);                       % Changing fuel flow [kg/min] 
   
   u1(t) = Thr_cl(3608.9,5249.3);                       % Changing thrust [N]
   u2(t) = 0;                                           % Changing bank angle [deg]
   u3(t) = 5;                                           % Changing flight path angle [deg]
   V_ver(t) = x4(t)*sin(u3(t)*pi/180);                  % Changing vertical speed [m/s]
   C_D = drag(3608.9,5249.3,x4(t));                     % Changing drag coefficient
   Cl = lift(3608.9,5249.3,x4(t));                      % Changing lift coefficient
   p = den(3608.9,5249.3);                              % Changing density [kg/m3]
end

if and (t >1,t<=61) % Perform cruise flight for t=60 minutes.
   x3(t) = 5249.3;                                      % Changing altitude [m] -> [ft]
   x4(t) = Cruise_Vel(5249.3);                          % Changing speed [m/s]
   x5(t) = 0;                                           % Changing head angle [deg]
   f(t) = fuel_cr(5249.3);                              % Changing fuel flow [kg/min] 
   
   u1(t) = Thr_cr(5249.3);                              % Changing thrust [N]
   u2(t) = 0;                                           % Changing bank angle [deg]
   u3(t) = 0;                                           % Changing flight path angle [deg]
   V_ver(t) = x4(t)*sin(u3(t)*pi/180);                  % Changing vertical speed [m/s]
   C_D = drag(5249.3,5249.3,x4(t));                     % Changing drag coefficient
   Cl = lift(5249.3,5249.3,x4(t));                      % Changing lift coefficient
   p = den(5249.3,5249.3);                              % Changing density [kg/m3]
end

if and (t >61,t<=62) % Turn with β=30 bank angle until heading is changed by η=270◦.
   x3(t) = 5249.3;                                      % Changing altitude [m] -> [ft]
   x4(t) = Cruise_Vel(5249.3);                          % Changing speed [m/s]
   x5(t) = 0:30:270;                                    % Changing head angle [deg]
   f(t) = fuel_cr(5249.3);                              % Changing fuel flow [kg/min] 
   
   u1(t) = Thr_cr(5249.3);                              % Changing thrust [N]
   u2(t) = 30;                                          % Changing bank angle [deg]
   u3(t) = 0;                                           % Changing flight path angle [deg]
   V_ver(t) = x4(t)*sin(u3(t)*pi/180);                  % Changing vertical speed [m/s]
   C_D = drag(5249.3,5249.3,x4(t));                     % Changing drag coefficient
   Cl = lift(5249.3,5249.3,x4(t));                      % Changing lift coefficient
   p = den(5249.3,5249.3);                              % Changing density [kg/m3]
end

if and (t >62,t<=63) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
   x3(t) = linspace(5249.3,3608.9,2);                   % Changing altitude [m] -> [ft]
   x4(t) = Des_Vel(5249.3,3608.9);                      % Changing speed [m/s]
   x5(t) = 270;                                         % Changing head angle [deg]
   f(t) = fuel_des(5249.3,3608.9);                      % Changing fuel flow [kg/min] 
   
   u1(t) = Thr_des(5249.3,3608.9);                      % Changing thrust [N]
   u2(t) = 0;                                           % Changing bank angle [deg]
   u3(t) = 4;                                           % Changing flight path angle [deg]
   V_ver(t) = x4(t)*sin(u3(t)*pi/180);                  % Changing vertical speed [m/s]
   C_D = drag(5249.3,3608.9,x4(t));                     % Changing drag coefficient
   Cl = lift(5249.3,3608.9,x4(t));                      % Changing lift coefficient
   p = den(5249.3,3608.9);                              % Changing density [kg/m3]
end

if and (t >63,t<=65) % Complete a 360◦ turn (loiter) at level flight.
   x3(t) = 3608.9;                                      % Changing altitude [m] -> [ft]
   x4(t) = Cruise_Vel(3608.9);                          % Changing speed [m/s]
   lon = [270 300 360 60 120 180 240 270];
   x5(t) = wrapTo360(lon);                              % Changing head angle [deg]
   f(t) = fuel_cr(3608.9);                              % Changing fuel flow [kg/min] 
   
   u1(t) = Thr_cr(3608.9);                              % Changing thrust [N] 
   u2(t) = 0;                                           % Changing bank angle [deg]
   u3(t) = 0;                                           % Changing flight path angle [deg]
   V_ver(t) = x4(t)*sin(u3(t)*pi/180);                  % Changing vertical speed [m/s]
   C_D = drag(3608.9,3608.9,x4(t));                     % Changing drag coefficient
   Cl = lift(3608.9,3608.9,x4(t));                      % Changing lift coefficient
   p = den(3608.9,3608.9);                              % Changing density [kg/m3]
end

if and (t >65,t<=66) % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
   x3(t) = linspace(3608.9,2624.67,2);                  % Changing altitude [m] -> [ft]
   x4(t) = Des_Vel(3608.9,2624.67);                     % Changing speed [m/s] 
   x5(t) = 270;                                         % Changing head angle [deg]
   f(t) = fuel_des(3608.9,2624.67);                     % Changing fuel flow [kg/min] 
   
   u1(t) = Thr_des(3608.9,2624.67);                     % Changing thrust [N] 
   u2(t) = 0;                                           % Changing bank angle [deg]
   u3(t) = 4.5;                                         % Changing flight path angle [deg]
   V_ver(t) = x4(t)*sin(u3(t)*pi/180);                  % Changing vertical speed [m/s]
   C_D = drag(3608.9,2624.67,x4(t));                    % Changing drag coefficient
   Cl = lift(3608.9,2624.67,x4(t));                     % Changing lift coefficient
   p = den(3608.9,2624.67);                             % Changing density [kg/m3]
end


for i=1:N

   x1(i+1)=x1(i)+h*dx1dt(x4(i),x5(i),u3(i));            % line 152 %%%%%%%%%%%%%%%%%%%%%
   x2(i+1)=x2(i)+h*dx2dt(x4(i),x5(i),u3(i));
   x3(i+1)=x3(i)+h*dx3dt(x4(i),u3(i));
   x4(i+1)=x4(i)+h*dx4dt(C_D(i),p(i),x6(i),x4(i),u3(i),u1(i));
   x5(i+1)=x5(i)+h*dx5dt(Cl(i),p(i),x6(i),x4(i),u2(i));
   x6(i+1)=x6(i)+h*dx6dt(f(i));
  
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);
legend('Euler');

When I run the code, this error displays:

Index exceeds the number of array elements (1).

Error in FlightPlan (line 152)
   x1(i+1)=x1(i)+h*dx1dt(x4(i),x5(i),u3(i));
   

After checking the Workspace, I saw that state variables (x1,x2,x3,x4,x5,x6) were updated as 1x2 matrices (first columns are initial conditions).

I want your help to utilize the Euler's Method on my code. How should I update the loops to plot Variable vs Time graphs? Is using (t) in loops is problematic? I need your help.

Thanks.

Answer 1

[Walter Roberson]

t=t0:h:tend;

if and (t >= 0,t<=1) % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.

t is a vector. and (t >= 0,t<=1) is a vector. if of a vector is considered true only if all of the values in the vector are non-zero; in the context of a logical test, that means that the condition had to be met for all of the values.

You need to switch from if to using logical indexing.

https://blogs.mathworks.com/loren/2013/02/20/logical-indexing-multiple-conditions/