# For loop broken ?

2 views (last 30 days)
Ali Khabibulayev on 12 Oct 2021
Commented: Ali Khabibulayev on 12 Oct 2021
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)

David Hill on 12 Oct 2021
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)?
Ali Khabibulayev on 12 Oct 2021
Ah, I was thinking - what is going on there, and it was that simple thing that I missed. Thank you!

Cris LaPierre on 12 Oct 2021
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?

R2021b

### Community Treasure Hunt

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

Start Hunting!