For loop broken ?

Hi all, I have a for loop in my matlab code which uses forward euler method to calculate ceratain parameters. The first parameter is rho (line 47) which is suppose to be dependent on h(n). However, it seems as the for loop skips the rho as it is not being updated when the loop iterates, althought, below parameter rhotest (line64) with manual h(55) in it gives a correct answer. Any help would be appreciated :) thank you.

clc; clear all; close all


%Aircraft parameters
g=9.80665;
m=41000;
W = m*g; 
Cd0 = 0.02
S=120%m^2
b=34 %m
AR=b^2/S
ef=0.82; % Efficiency Factor
K=0.03 % combined induced drag factor, = 𝑘 ∕ (𝜋 𝐴𝑅)  

Ta = 110000; %Thrust Available
Tp=1; %Thrust Percentage
Tu=Ta*Tp; % (Thrust Used)


Cl=sqrt((pi()/3)*ef*AR*Cd0);
Cd = K+Cd0*Cl^2
%Cd=Cd0+((Cl^2)/(pi()*ef*AR));


tf = 60; %Final time
dt = 1; %Time step
t = 0:dt:tf;

V0 = 85; %Initial Velocity
X0 = 0; %Initial Displacement
h0 = 0; %Initial Height
gamma0 = 0.0872665 % 5deg initial climb. 


%Pre-Fill with zeros in order to avoid buildup in the loop.
V = zeros(1,length(t)) 
V(1) = V0; % Initial velocity in dt frame
X = zeros(1,length(t))
X(1) = X0; % Initial X displacement in dt frame 
h = zeros(1,length(t));
h(1) = h0; %Initial h displacement in dt frame
gamma = zeros(1,length(t));
gamma(1) = gamma0; %Initial gamma in dt frame

for n=2:length(t)

    rho = (20-h(n)/1000)/(20+h(n)/1000)*1.225;

    D = Cd*S*0.5*rho*V(n-1)^2;

    V(n) = V(n-1)+((Ta-D-W*sin(gamma(n-1)))/m)*dt;

    ROC = ((Ta*V(n)-D*V(n))/W);

    gamma(n)= asin(ROC/V(n));

    X(n) = X(n-1)+V(n)*dt*cos(gamma(n));

    h(n) = h(n-1)+V(n)*dt*sin(gamma(n));

end

%% equation checks
rhotest= (20-h(55)/1000)/(20+h(55)/1000)*1.225;
tst2=((Ta-D-W*sin(gamma0))/m)*dt;
%V(n) = V(n-1)+(1/m)*(Ta-D-W*sin(gamma(n-1)));
tst=(1/m)*(Ta-D-W*sin(gamma(1)));

%% Plots 

plot(X,h)

Answer 1

[David Hill]

rho = (20-h(n)/1000)/(20+h(n)/1000)*1.225;%h(n) is always going to be zero
rho = (20-h(n-1)/1000)/(20+h(n-1)/1000)*1.225;%did you want h(n-1)?

[BioZ]

Ah, I was thinking - what is going on there, and it was that simple thing that I missed. Thank you!

Answer 2

[Cris LaPierre]

rho is calculated using the current value of h(n), which your code has defined as zero. Therefore, rho is always going to have the same value.

h = zeros(1,length(t));

Perhaps you meant to use h(n-1) instead?