Index exceeds the number of array elements for a For loop

I'm trying to find the best mass fuel burn rate for a rocket with my code using the formula:

acceleration = Thrust/mass - Weight - Drag.

The code is designed so that if the fuel runs out: acceleration = -weight + drag (taking upwards the positive direction). And if drag > (thrust - weight) the acceleration is 0.

I keep getting this error however:

Index exceeds the number of array elements. Index must not exceed 1

for the line: ay_2(i) = T(i-1) - W_2(i-1) - D_2(i-1);

What is going wrong?

%%

clear

clc

%time

timestep = 1;

t_end = 120;

t_start = 0;

t = t_start:timestep:t_end;

G = 6.67408e-11; % Universal Gravitational Constant

M = 5.9722e24; % Mass of the earth

R = 6371e3; % Radius of the earth(m)

A = sqrt((1)^2+(0.25/2)^2)*0.25/2*pi; % Area of the rocket head

Cd = 0.4; % Drag Coefficient for the rocket

g_0 = 9.81;

v_e = 4.5E3;

%mass

m_engine = 0.82;

m_payload = 0.18;

m_fuel_initial = 300;

m_to = m_payload + m_fuel_initial + m_engine;

m(1) = m_to;

%intial values

r(1) = R;

m(1)= m_to;

rho_2(1) = 1.225;

g(1) = 9.81;

vy_2(1) = 0;

W_2(1) = m_to*g(1);

%dmdt = 2

dmdt_2 = 2; %fuel burn rate

T_max_2(1) = v_e*dmdt_2;

ay_2(1) = T_max_2/m_to - 9.81;

vy_2(1) = 0;

y_2(1) = 0;

r_2(1) = R;

m_fuel_2(1) = 300;

D_2(1) = 0;

for i = 2:length(t)

if m_fuel_2 > 0

m(i) = m(i-1) - timestep*dmdt_2;

m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;

T(i) = T_max_2;

else m_fuel <= 0

m(i) = m(i-1) - timestep*dmdt_2;

m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;

T(i) = 0;

ay_2 = -W_2 + D(i-1);

end

if D_2 < T_max_2 - W_2

ay_2(i) = T(i-1) - W_2(i-1) - D_2(i-1);

else D_2 > T_max_2 + W_2;

ay_2(i) = 0;

end

D(i) = 0.5*A*Cd*rho_2(i-1)*vy_2(i-1);

W(i) = G*M/(r(i-1))^2;

vy_2(i) = vy_2(i-1) + ay_2(i)*timestep;

y_2(i) = y_2(i-1) + vy_2(i)*timestep;

r_2(i) = y_2(i) + R;

end

Answer 1

[Mathieu NOE]

hello Asit

there was some obvious bugs in your recursion / variables used

and also I corrected the "else" lines which do not need any mirrored statement - as explained below (how to use if / elseif / else )

if expression

statements

elseif expression

statements

else

statements

end

last but not least , the recursion for rho_2 is missing , so I simply used the initial (fix) value

corrected code :

%%
clear 
clc
%time
timestep = 1;
t_end = 120;
t_start = 0;
t = t_start:timestep:t_end;
G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = sqrt((1)^2+(0.25/2)^2)*0.25/2*pi;                 % Area of the rocket head
Cd = 0.4;               % Drag Coefficient for the rocket
g_0 = 9.81;
v_e = 4.5E3;
%mass
m_engine = 0.82;
m_payload = 0.18;
m_fuel_initial = 300;
m_to = m_payload + m_fuel_initial + m_engine;
m(1) = m_to;
%intial values
r(1) = R;
m(1)= m_to;
rho_2(1) = 1.225;
g(1) = 9.81;
vy_2(1) = 0;
W_2(1) = m_to*g(1);
%dmdt = 2
dmdt_2 = 2; %fuel burn rate
T_max_2(1) = v_e*dmdt_2;
ay_2(1) = T_max_2/m_to - 9.81;
vy_2(1) = 0; 
y_2(1) = 0;
r_2(1) = R;
m_fuel_2(1) = 300; 
D_2(1) = 0;
for i = 2:length(t)
    
   if m_fuel_2 > 0
    m(i) = m(i-1) - timestep*dmdt_2;
    m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;
    T(i) = T_max_2;
%    else m_fuel <= 0;
   else 
    m(i) = m(i-1) - timestep*dmdt_2;
    m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;
    T(i) = 0;
    ay_2 = -W_2 + D(i-1);
   end
   
    if D_2 < T_max_2 - W_2
        ay_2(i) = T(i-1) - W_2(i-1) - D_2(i-1);
%     else D_2 > T_max_2 + W_2;
    else
        ay_2(i) = 0;
    end
%     D(i) = 0.5*A*Cd*rho_2(i-1)*vy_2(i-1); % replace D (not defined) by D_2
%     D_2(i) = 0.5*A*Cd*rho_2(i-1)*vy_2(i-1); % no equation of recursion  for rho2 ? 
    D_2(i) = 0.5*A*Cd*rho_2*vy_2(i-1); % fix rho2 value
%     W(i) = G*M/(r(i-1))^2;  % replace W (not defined) by W_2

%     W_2(i) = G*M/(r(i-1))^2;
    W_2(i) = G*M/(r_2(i-1))^2;  % replace r (not defined) by r_2
    vy_2(i) = vy_2(i-1) + ay_2(i)*timestep;
    y_2(i) = y_2(i-1) + vy_2(i)*timestep;
    r_2(i) = y_2(i) + R;
end

[Asit Rahman]

It almost works. The only problem is that the altitude continues to increase, when it should decrease after running out of fuel

[Mathieu NOE]

ok

i tried to first find out more bugs..like missing index in recursion , and also errors in equations betwen acceleration , thrust (is a force => must be divided by mass) , drag force (is a force => must be divided by mass)

also a drag force is proportionnal to the squared velocity

I also introduced the equation between the density vs altitude which is crucial to have a more accurate drag force

unfortunately my look up table is valide only up to 10,000 m , so I tried to extrapolate that for higher altitudes , but this can create some strange results once the latitude is above 10,000 m

updated code is :

%%
clear 
clc
%time
timestep = 1;
t_end = 120;
t_start = 0;
t = t_start:timestep:t_end;
G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = sqrt((1)^2+(0.25/2)^2)*0.25/2*pi;                 % Area of the rocket head
Cd = 0.4;               % Drag Coefficient for the rocket
g_0 = 9.81;
v_e = 4.5E3;
%mass
m_engine = 0.82;
m_payload = 0.18;
m_fuel_initial = 300;
m_to(1) = m_payload + m_fuel_initial + m_engine;
m(1) = m_to;
%intial values
r(1) = R;
m(1)= m_to;
rho_2(1) = 1.225;
g(1) = 9.81;
vy_2(1) = 0;
% W_2(1) = m_to*g(1);
W_2(1) = g(1);
%dmdt = 2
dmdt_2 = 3; %fuel burn rate
T_max_2(1) = v_e*dmdt_2;
ay_2(1) = T_max_2/m_to - 9.81;
vy_2(1) = 0; 
y_2(1) = 0;
r_2(1) = R;
m_fuel_2(1) = 300; 
D_2(1) = 0;

%% air density (kgs/m3) vs altitude (m) look up table
% this table is know to be true up to alt = 10,000 m (extrapolated above more or less accurately)
alt = [0 1 2 3 4 5 6 7 8 10 20 50]*1e3;		%altitude in m
rho = 1.225*[1 0.907 0.822 0.742 0.669 0.601 0.538 0.481 0.428 0.337 0.01 0]; % air density (kgs/m3)


% for debugg only
cond_m_fuel(1) = 1;
cond_D_2(1) = 1;

for i = 2:length(t)
    
    if m_fuel_2 > 0
        cond_m_fuel(i) = 1;
        m(i) = m(i-1) - timestep*dmdt_2;
        m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;
        T(i) = T_max_2;
        
        %    else m_fuel <= 0;
    else 
        cond_m_fuel(i) = 0;
        %     m(i) = m(i-1) - timestep*dmdt_2;
        m(i) = m(i-1); % no more fuel burn => total mass constant
        %     m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2; % I think this is not correct
        m_fuel_2(i) = m_fuel_2(i-1); % no more fuel burn = constant mass
        T(i) = 0; % no fuel burn => no thrust
        ay_2(i) = -W_2(i-1) - D_2(i-1)./m(i);% here : indexing was missing
        
    end
    
    if D_2(i-1) < T_max_2 - W_2(i-1).*m(i) % here : indexing was missing, product W2 (gravity) by mass was missing (to be a force)
        cond_D_2(i) = 1;
        ay_2(i) = T(i)./m(i) - W_2(i-1) - D_2(i-1)./m(i); % acceleration = thrust force / mass - drag force / mass - gravity (W2)(division by mass was missing)
        %     else D_2 > T_max_2 + W_2;
    else
        cond_D_2(i) = 0;
        ay_2(i) = 0;
    end
    %     D(i) = 0.5*A*Cd*rho_2(i-1)*vy_2(i-1); % replace D (not defined) by D_2
    %     D_2(i) = 0.5*A*Cd*rho_2(i-1)*vy_2(i-1); % no equation of recursion  for rho2 ? 
    rho_2 = interp1(alt,rho,y_2(i-1)); % update air density vs altitude (from ground)
    D_2(i) = 0.5*A*Cd*rho_2*(vy_2(i-1)).^2; % fix rho2 value, velocity must be squared (Formula to calculate Drag Force is given by d = 1/2 * ρ * u² * A * Cd)
    %     W(i) = G*M/(r(i-1))^2;  % replace W (not defined) by W_2
    %     W_2(i) = G*M/(r(i-1))^2; % done
    %     W_2(i) = G*M/(r_2(i-1))^2;  % replace r (not defined) by r_2
    W_2(i) = G*M/(r_2(i-1))^2;  % acceleration (gravity) vs distance to earth center
    vy_2(i) = vy_2(i-1) + ay_2(i)*timestep;
    y_2(i) = y_2(i-1) + vy_2(i)*timestep;
    r_2(i) = y_2(i) + R;
end
subplot(421),plot(t,ay_2);legend('ay 2');
subplot(422),plot(t,T);legend('T');
subplot(423),plot(t,m_fuel_2);legend('m fuel 2');
subplot(424),plot(t,D_2);legend('D 2');
subplot(425),plot(t,D_2 - (T_max_2 - W_2));legend('D 2 - (T max 2 - W 2)');
subplot(426),plot(t,vy_2);legend('vy 2');
subplot(427),plot(t,cond_m_fuel,t,cond_D_2);legend('cond m fuel','cond D 2');
subplot(428),plot(t,y_2);legend('y 2');

[Mathieu NOE]

now , I was skeptical about what we are doing here (the second if loop)

    if D_2(i-1) < T_max_2 - W_2(i-1).*m(i) % here : indexing was missing, product W2 (gravity) by mass was missing (to be a force)
        cond_D_2(i) = 1;
        ay_2(i) = T(i)./m(i) - W_2(i-1) - D_2(i-1)./m(i); % acceleration = thrust force / mass - drag force / mass - gravity (W2)(division by mass was missing)
        %     else D_2 > T_max_2 + W_2;
    else
        cond_D_2(i) = 0;
        ay_2(i) = 0;
    end

to me that sounds like we force the acceleration to zero because the drag force exceeds thrust minus gravity

what is the physical reason behind ?

[Mathieu NOE]

according to my previous comment I tried to make the code different :

the "strange"(IMHO) second if loop has been replaced by the fact that the acceleration can never be more negative that -W2 - so because we hav not an accurate enough model of the drag, as soon as the thrust falls to zero, the remaining acceleration is only -W2 (and the rocket would fall back towards earth). The max altitude is not bigh enough vs earth radius to make W2 significantly lower than 9.81 (only a few % drop)

%%
clear 
clc
%time
timestep = 1;
t_end = 300;
t_start = 0;
t = t_start:timestep:t_end;
G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = sqrt((1)^2+(0.25/2)^2)*0.25/2*pi;                 % Area of the rocket head
Cd = 0.4;               % Drag Coefficient for the rocket
g_0 = 9.81;
v_e = 4.5E3;
%mass
m_engine = 0.82;
m_payload = 0.18;
m_fuel_initial = 300;
m_to(1) = m_payload + m_fuel_initial + m_engine;
m(1) = m_to;
%intial values
r(1) = R;
m(1)= m_to;
rho_2(1) = 1.225;
g(1) = 9.81;
vy_2(1) = 0;
% W_2(1) = m_to*g(1);
W_2(1) = g(1);
%dmdt = 2
dmdt_2 = 3; %fuel burn rate
T_max_2(1) = v_e*dmdt_2;
ay_2(1) = T_max_2/m_to - 9.81;
vy_2(1) = 0; 
y_2(1) = 0;
r_2(1) = R;
m_fuel_2(1) = 300; 
D_2(1) = 0;

%% air density (kgs/m3) vs altitude (m) look up table
% this table is know to be true up to alt = 10,000 m (extrapolated above more or less accurately)
alt = [0 1 2 3 4 5 6 7 8 10 20 100 1000]*1e3;		%altitude en m
rho = 1.225*[1 0.907 0.822 0.742 0.669 0.601 0.538 0.481 0.428 0.337 0.2 0.03 0.02];


% for debugg only
cond_m_fuel(1) = 1;
cond_D_2(1) = 1;

for i = 2:length(t)

   if m_fuel_2 > 0
    cond_m_fuel(i) = 1;
    m(i) = m(i-1) - timestep*dmdt_2;
    m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;
    T(i) = T_max_2;

%    else m_fuel <= 0;
   else 
    cond_m_fuel(i) = 0;
%     m(i) = m(i-1) - timestep*dmdt_2;
    m(i) = m(i-1); % no more fuel burn => total mass constant
%     m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2; % I think this is not correct
    m_fuel_2(i) = m_fuel_2(i-1); % no more fuel burn = constant mass
    T(i) = 0; % no fuel burn => no thrust
%     ay_2(i) = -W_2(i-1) - D_2(i-1)./m(i);% here : indexing was missing

   end
   
   ay_2(i) = T(i)./m(i) - W_2(i-1) - D_2(i-1)./m(i); % acceleration = thrust force / mass - drag force / mass - gravity (W2)(division by mass was missing)
   
   if ay_2(i)<- W_2(i-1)
       ay_2(i) = - W_2(i-1);
   end
   
%     if D_2(i-1) < T_max_2 - W_2(i-1).*m(i) % here : indexing was missing, product W2 (gravity) by mass was missing (to be a force)
%         cond_D_2(i) = 1;
%         ay_2(i) = T(i)./m(i) - W_2(i-1) - D_2(i-1)./m(i); % acceleration = thrust force / mass - drag force / mass - gravity (W2)(division by mass was missing)
% %     else D_2 > T_max_2 + W_2;
%     else
%         cond_D_2(i) = 0;
%         ay_2(i) = 0;
%     end

%     D(i) = 0.5*A*Cd*rho_2(i-1)*vy_2(i-1); % replace D (not defined) by D_2
%     D_2(i) = 0.5*A*Cd*rho_2(i-1)*vy_2(i-1); % no equation of recursion  for rho2 ? 
    rho_2 = interp1(alt,rho,y_2(i-1)); % update air density vs altitude (from ground)
    D_2(i) = 0.5*A*Cd*rho_2*(vy_2(i-1)).^2; % fix rho2 value, velocity must be squared (Formula to calculate Drag Force is given by d = 1/2 * ρ * u² * A * Cd)
%     W(i) = G*M/(r(i-1))^2;  % replace W (not defined) by W_2
%     W_2(i) = G*M/(r(i-1))^2; % done
%     W_2(i) = G*M/(r_2(i-1))^2;  % replace r (not defined) by r_2
    W_2(i) = G*M/(r_2(i-1))^2;  % acceleration (gravity) vs distance to earth center
    vy_2(i) = vy_2(i-1) + ay_2(i)*timestep;
    y_2(i) = y_2(i-1) + vy_2(i)*timestep;
    r_2(i) = y_2(i) + R;
end
figure(1)
subplot(421),plot(t,ay_2,t,-W_2);legend('ay 2','gravity W2');
subplot(423),plot(t,m_fuel_2);legend('m fuel 2');
subplot(425),plot(t,T);legend('T');
subplot(427),plot(t,m);legend('m');
subplot(422),plot(t,D_2);legend('D 2');
subplot(424),plot(t,vy_2);legend('vy 2');
subplot(426),plot(t,y_2,t,1e4.*ones(size(t)),'--r');legend('y 2','limit of accurate rho model');

[Asit Rahman]

I see the problem with the previous code now. The code was supposed to show that when the rocket ran out of fuel, it would be subjected to just weight and drag and fall back to ground. Unfortunately when I take away that part of the code, I'm left with another index error:

clear

clc

%time

timestep = 1;

t_end = 60;

t_start = 0;

t = t_start:timestep:t_end;

G = 6.67408e-11; % Universal Gravitational Constant

M = 5.9722e24; % Mass of the earth

R = 6371e3; % Radius of the earth(m)

A = sqrt((1)^2+(0.25/2)^2)*0.25/2*pi; % Area of the rocket head

Cd = 0.4; % Drag Coefficient for the rocket

g_0 = 9.81;

v_e = 4.5E3;

%mass

m_engine = 0.82;

m_payload = 0.18;

m_fuel_initial = 300;

m_to = m_payload + m_fuel_initial + m_engine;

m(1) = m_to;

%intial values

r(1) = R;

m(1)= m_to;

rho_2(1) = 1.225;

g(1) = 9.81;

vy_2(1) = 0;

W_2(1) = m_to*g(1);

%dmdt = 2

dmdt_2 = 10; %fuel burn rate

T_max_2(1) = v_e*dmdt_2;

ay_2(1) = T_max_2/m_to - 9.81;

vy_2(1) = 0;

y_2(1) = 0;

r_2(1) = R;

m_fuel_2(1) = 300;

D_2(1) = 0;

for i = 2:length(t)

if m_fuel_2 > 0

m(i) = m(i-1) - timestep*dmdt_2;

m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;

T(i) = T_max_2;

ay_2 = T_max_2 - g(i-1) + D_2(i-1);

else

m(i) = m(i-1) - timestep*dmdt_2;

m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;

T(i) = 0;

ay_2 = -g(i-1) + D_2(i-1);

end

rho(i)=1.225*10.^((-3*y_2(i-1)/50000));

g_2(i) = G*M/(r_2(i-1))^2;

D_2(i) = 0.5*A*Cd*rho_2(i-1)*vy_2(i-1);

vy_2(i) = vy_2(i-1) + ay_2(i)*timestep;

y_2(i) = y_2(i-1) + vy_2(i)*timestep;

r_2(i) = y_2(i) + R;

end

subplot(3,2,1)

plot(t,ay_2)

xlabel('time s')

ylabel('a_y m/s^2')

subplot(3,2,2)

plot(t,D_2)

xlabel('time s')

ylabel('Drag kN')

subplot(3,2,3)

plot(t,vy_2)

xlabel('time s')

ylabel('v_y km/s')

subplot(3,2,4)

plot(t,m_fuel_2)

xlabel('time s')

ylabel('fuel mass kg')

subplot(3,2,5)

plot(t,y_2)

xlabel('time s')

ylabel('alt km')

subplot(3,2,6)

plot(t,g)

xlabel('time s')

ylabel('g m/s^2')

[Mathieu NOE]

hello

here some corrections to your code (see remarks)

clear 
clc
%time
timestep = 1;
t_end = 260;
t_start = 0;
t = t_start:timestep:t_end;
G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = sqrt((1)^2+(0.25/2)^2)*0.25/2*pi;                 % Area of the rocket head
Cd = 0.4;               % Drag Coefficient for the rocket
g_0 = 9.81;
v_e = 4.5E3;
%mass
m_engine = 0.82;
m_payload = 0.18;
m_fuel_initial = 300;
m_to = m_payload + m_fuel_initial + m_engine;
m(1) = m_to;
%intial values
r(1) = R;
m(1)= m_to;
rho_2(1) = 1.225;
g_2(1) = 9.81;
vy_2(1) = 0;
% W_2(1) = m_to*g_2(1); % not used
%dmdt = 2
dmdt_2 = 2; %fuel burn rate
T_max_2(1) = v_e*dmdt_2;
ay_2(1) = T_max_2/m_to - 9.81;
vy_2(1) = 0; 
y_2(1) = 0;
r_2(1) = R;
m_fuel_2(1) = 300; 
D_2(1) = 0;
for i = 2:length(t)
    
   if m_fuel_2(i-1) > 0   % indexing was missing
    m(i) = m(i-1) - timestep*dmdt_2;
    m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;
    T(i) = T_max_2;
    ay_2(i) = T_max_2./m(i) - g_2(i-1) - D_2(i-1)./m(i); % here indexing was missing : ay_2(i) instead of ay_2 ;
                                        % drag force is negative =>changed sign 
                                        % drag and thrust forces must be divided by mass to give acceleration terms
   else 
%     m(i) = m(i-1) - timestep*dmdt_2;
    m(i) = m(i-1);
%     m_fuel_2(i) = m_fuel_2(i-1) - timestep*dmdt_2;
    m_fuel_2(i) = 0;
    T(i) = 0;
    ay_2(i) = -g_2(i-1) - D_2(i-1)./m(i); % same comments as above
   end
  
   % restrict acceleration to structural limits (10 G)  (reduce progressively thrust and fuel burn)
  if ay_2(i)>100
      ay_2(i)=100;
  end
  
    vy_2(i) = vy_2(i-1) + ay_2(i)*timestep;
    y_2(i) = y_2(i-1) + vy_2(i)*timestep;
    r_2(i) = y_2(i) + R;
    % now we can update rho, g_2, D_2 once vy_2, y_2, r_2
    rho(i)=1.225*10.^((-3*y_2(i)/50000)); 
    g_2(i) = G*M/(r_2(i))^2;
    D_2(i) = 0.5*A*Cd*rho(i)*(vy_2(i))^2; % replaced rho_2 (why ??) by rho ( 2 lines above); also darg is normally proportionnal to squared velocity

end
subplot(3,2,1)
plot(t,ay_2)
xlabel('time s')
ylabel('a_y m/s^2')
subplot(3,2,2) 
plot(t,D_2/1000)
xlabel('time s')
ylabel('Drag kN')
subplot(3,2,3) 
plot(t,vy_2/1000)
xlabel('time s')
ylabel('v_y km/s')
subplot(3,2,4) 
plot(t,m_fuel_2)
xlabel('time s')
ylabel('fuel mass kg')
subplot(3,2,5) 
plot(t,y_2/1000)
xlabel('time s')
ylabel('alt km')
subplot(3,2,6) 
plot(t,g_2)
xlabel('time s')
ylabel('g m/s^2')  

[Asit Rahman]

This works perfectly! Thank you so much!

[Mathieu NOE]

My pleasure !