Index exceeds the number of array elements (391).

Hi all. That is my function:

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)<=1499)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL1;

elseif (1500<=H_climb(k1)&&H_climb(k1)<=2999)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL2;
        
elseif (3000<=H_climb(k1)&&H_climb(k1)<=3999)
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL3; 
        
elseif (4000<=H_climb(k1)&&H_climb(k1)<=4999)
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL4;
        
elseif (5000<=H_climb(k1)&&H_climb(k1)<=5999)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL5;
        
elseif (6000<=H_climb(k1)&&H_climb(k1)<=9999)
        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]       line 60 %%%%%%%%%%%%%%%%%%%%%%%%%%%


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

Whenever I try to run the code like this:

>> [Vtas_cl] =Velocities(3608.9,5249.3)
Index exceeds the number of array elements (391).

Error in Velocities (line 60)
Vcas_cl(k1) = Vnom_climb_jet(k1)* 0.514;   % [kn] -> [m/s]

Another condition:

>> [Vtas_cl] =Velocities(3608.9,3999.9)
Index exceeds the number of array elements (391).

Error in Velocities (line 60)
Vcas_cl(k1) = Vnom_climb_jet(k1)* 0.514;   % [kn] -> [m/s]

And another one:

>> [Vtas_cl] =Velocities(3608.9,3999.8)

Vtas_cl =

   1.0e+02 *

  Columns 1 through 6

   1.161544478045788   1.161561284482268   1.161578091258861   1.161594898375583   1.161611705832432   1.161628513629424

  Columns 7 through 12

   1.161645321766563   1.161662130243859   1.161678939061322   1.161695748218955   1.161712557716770   1.161729367554776
   

I get the results from column 1 through column 391. What should I change to fix the code?

Thanks.

[Stephen23]

The problem is your IF/ELSEIF logic inside the FOR loop.

Consider what happens if none of those conditions are fulfilled: would any value get added to Vnom_climb_jet? (hint: no)

For example, consider a value like H_climb(k1) = 1499.5, what would your code do? (tip: you should test each condition only once, not twice like you are doing, to avoid this problem).

Preallocating H_climb is probably a good idea too, e.g. using NaNs:

https://www.mathworks.com/help/matlab/matlab_prog/preallocating-arrays.html

[Turgut Ataseven]

Thanks. Here is what I did.

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

[Stephen23]

@Turgut Ataseven: your conditions ares still too complex, which makes your code harder to understand and more liable to bugs (as your question demonstrates). As I wrote earlier, check each conditions just once, not duplicated as you are doing. For example:

if 0<=H_climb(k1) && H_climb(k1)<1500
    ..
elseif H_climb(k1)<3000 % do NOT check 1500<=H_climb(k1) here, the previous condition ensures this!
    ..
elseif H_climb(k1)<4000 % do NOT check 3000<=H_climb(k1) here, the previous condition ensures this!
    ..
    etc
end

Your code is still fragile as you do not have an ELSE or a preallocated array.

Answer 1

[Image Analyst]

This works. I preallocated the array and fixed the if statements, and put in an else to catch unhandled "if" conditions (never encountered any though).

clc;    % Clear the command window.
fprintf('Beginning to run %s.m ...\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clearvars;
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;

% Call the function
Vtas_cl1 = Velocities(3608.9,5249.3)
Vtas_cl2 = Velocities(3608.9,3999.9)

% Define the function
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;
% Preallocate array.
Vnom_climb_jet = zeros(1, length(H_climb));
for k1 = 1:length(H_climb)
    if (0 <= H_climb(k1)&&H_climb(k1) <= 1499)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL1;
    elseif H_climb(k1) <= 2999
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL2;
    elseif H_climb(k1) <= 3999
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL3;
    elseif H_climb(k1) <= 4999
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL4;
    elseif H_climb(k1) <= 5999
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL5;
    elseif H_climb(k1) <= 9999
        Vnom_climb_jet (k1)= min(Vcl_1,250);
    elseif 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;
    else
        uiwait(errordlg('Unhandled if condition'))
    end
    Vcas_cl(k1) = Vnom_climb_jet(k1)* 0.514;   % [kn] -> [m/s]       line 60 %%%%%%%%%%%%%%%%%%%%%%%%%%%
    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