For loop giving me wrong values

I am trying to find the output for different parameters. It gives me zeros for masses that can not be zero. I think I have a problem with the for loop. Please help!

Intro: I have 9 values for maximum hight, h_max, 9 values for maximum accelaration a_max, 9 values for structure margin sm_max, uploaded in the document "Data".

I also have 9 values for maximum length, L_max, 9 values for maximum diameter , D_max.

The code is functional in returning the values for

R_out
Weq_out
t_b
p_0
delta_ratio
X_CG
X_CP

But fails when it gets to

M_p
M_0
M_s

clear all
%% Calculating the outputs for all the 9 D_max, L_max 
%For Rocket
%Payload mass:
M_L=1;%kg
%Number of fins:
N=3;
%Shell density:
rho_s=2700; %kg/m^3
%Propellant density:
rho_p=1772; %kg/m^3
%Shell working stress:
sig_s=60*10^6; %pa
%Gravity:
g=9.81; % m/s^2
%For air:
%Atmorspheric TeM_perature at sea level:
T_atm=298;
%Specific heat ratio of air:
gamma=1.4;
%Density of air at sea level:
rho=1.225;
%Pressure at sea level:
P_a=101325; %pa

% Here I start to run the tests for the 9 values 
Data=readmatrix('Data');% Each col is 1 test
L_max = [1.0755    1.0651    1.2883    0.7842    1.4316    1.7396    1.8300    0.6058    1.7892];
D_max = [0.0929    0.0853    0.0769    0.0709    0.1065    0.0753    0.0749    0.0825    0.0957];

%Here I set up the output matrices for the loop
R_out=zeros(1,9);
Weq_out=zeros(1,9);
t_b=zeros(1,9);
p_0=zeros(1,9);
delta_ratio=zeros(1,9)
D_out=D_max;
L_out=L_max;
X_CG=zeros(1,9);
X_CP=zeros(1,9);
M_p=zeros(1,9);
M_0=zeros(1,9);
M_s = zeros(1,9);
epsilon_c=zeros(1,9);

for loop=1:9
    disp(num2str(loop));
    h_max=Data(1,loop);
    a_max=Data(2,loop);
    SM=Data(3,loop);
    L=L_max(1,loop);
    D=D_max(1,loop);
    R_max=1+a_max;
    
    
    W_eq=sqrt((h_max*g)/(((log(R_max)/2)*(log(R_max)-2))+((R_max-1)/R_max)));
    t_bmax=(R_max-1)*W_eq/(g*R_max);
    M_eq=W_eq/sqrt(gamma*287*T_atm);
    P_c=P_a*(1+(((gamma-1)/2)*M_eq^2))^(gamma/(gamma-1));
    P_0_a=P_c/P_a;
            delta= (P_c/(2*sig_s))*D; % thickness of shell m
            M_n=delta*rho_s*pi*D*(D+sqrt(D^2+(D^2/4)));
            M_f=((D^2)/2)*delta*rho_s;
            M_f_b=(pi*D*rho_s*D*delta);
            M_s=(pi*D*rho_s*L*delta)+M_n+M_f+M_f_b;%%%%%%%%%%
            M_p=(R_max-1)*(M_s+M_L);
            L_p=(M_p/(pi*D^2*rho_p/4));
            
            %% Ignore this part! 
                % Center of Pressure for Rocket Nose
                X_n=(2/3)*D;%m
                CN_n=2;
                %Center of Pressure for Rocket Fin
                a=D;%m
                s=D;%m
                b=0;
                m=a-b;
                fin_hyp=sqrt(2)*D;%m
                X_f=D+L;%m
                delta_X_f=((m*(a+2*b))/(3*(a+b)))+((1/6)*(a+b-((a*b)/(a+b))));%m
                X_f_dash=X_f+delta_X_f;%m
                CN_f=(4*N*(s/D)^2)/(1+sqrt(1+((2*fin_hyp)/(a+b))^2));%m
                k_fb=1+((D/2)/(s+(D/2)));%m
                CN_fb=k_fb*CN_f; %m;
                X_cp=((CN_n*X_n)+(CN_fb*X_f_dash))/(CN_n+CN_fb);
                X_ncg=2*D/3;
                
                %% Porblem begins here ! 
                
                M_n=delta*rho_s*pi*D*(D+sqrt(D^2+(D^2/4)));
                M_c=M_L+M_n;
                
                %Center of Gravity for Rocket Fin
                X_f_cg=(2*D/3)+L+D;
                M_f=((D^2)/2)*delta*rho_s;
                
                %Center of Gravity for Rocket Tube 1
                L_1=L+D-L_p;
                xL_1_cg=(L_1/2)+D;
                M_L_1=pi*D*L_1*delta*rho_s;
                
                %Center of Gravity for Rocket Tube 2
                L_2=L_p;
                xL_2_cg=(L_2/2)+D+L_1;
                M_L_2=(pi*D*L_2*delta*rho_s)+M_p;
                % Center of gravity final
                M_f_b = (pi*D*rho_s*D*delta)
                X_cg=((M_c*X_ncg)+(M_L_1*xL_1_cg)+(xL_2_cg*M_L_2)+(X_f_cg*M_f*3))/(M_c+M_L_1+M_L_2+(M_f*3));
                M_s=(pi*D*rho_s*delta)+M_n+M_f+M_f_b
                M_p=(R_max-1)*(M_s+M_L);
                M_0=M_s+M_p+M_L 
                % PROBLEM: You see here M_0 is the summation of M_s M_p and M_L 
                % but M_L is defined previously = 1 
                % so the output here can not be zero yet it returns zero
               
                
                epsilon=M_s/(M_s+M_p);
                lamda=M_L/(M_s+M_p);
                
                
R_out(1,loop)=R_max
Weq_out(1,loop)=W_eq
t_b(1,loop)=t_bmax
p_0(1,loop)=P_0_a
delta_ratio(1,loop)=delta
D_out=D_max
L_out=L_max
X_CG(1,loop)=X_cg
X_CP(1,loop)=X_cp
M_p(1,loop)=M_p
M_0(1,loop)=M_0
epsilon_c(1,loop)= epsilon
M_s(1,loop)=M_s
end

[Rik]

The code in your loop is not clearly explained by the comments. It is fairly complex, so it is not obvious where the issue is.

[Bjarke Skogstad Larsen]

I think you have some problems with (old) workspace values interfering. Try starting out with clear command to clear the workspace and work through the problem from there.

[Mohsen Nashel]

Rik, the loop is to find the solution to all the equations listed in it through 9 scenarios listed in L_max, D_max and the data matrix, which is obvious in the loop. There is nothing complicated about it, it's just a bunch of equation iterating through the 9 given values of 5 input paramters.

Bjarke Skogstad Larsen

I think you are right but I used clear and correct the code and ran it it is till giving the same values. Which is why I think it might be a problem of outputs of the loop itself

i.e :

Look at the equation

M_0=M_s+M_p+M_L

M_L is defined as 1 at the beginging of the code, which means M_0 can never be zero. Does that mean it might be getting old values going though the loop? Indexing that using (ii) maybe ? I don't know !

[Rik]

I agree, nothing complicated:

W_eq=sqrt((h_max*g)/(((log(R_max)/2)*(log(R_max)-2))+((R_max-1)/R_max)));
X_cg=((M_c*X_ncg)+(M_L_1*xL_1_cg)+(xL_2_cg*M_L_2)+(X_f_cg*M_f*3))/(M_c+M_L_1+M_L_2+(M_f*3));

Obvious, right?

If the results don't match what you think they should be, there are two options

1. You made a mistake in implementing your algorithm. (which we can't check due to the lack of comments explaining your code)
2. You thought something was true, but it turns out to be incorrect. (which we can't check because you aren't giving any context)

Have a read here and here. It will greatly improve your chances of getting an answer.

[Mohsen Nashel]

My bad! Check the notations

The probelm where I shown in the code. It has to be something I am not getting about the for loop because it might be retaining old values when going through the loop. So the equations are not the problem, nor the alrgrithm, just he loop because it runs individually fine.

[Rik]

Have you tried replacing clear all with clear and using breakpoints to step through your code line by line?

Answer 1

[John D'Errico]

I tried to look through the code for a few minutes now, looking for something obvious. So then I downloaded your data, put your code into MATLAB. Nothing seemed clear there either. Then I ran the code.

Sadly, it is dumping a lot of crap to the screen that means nothing to me, since I have no idea what is reasonable and what is not.

Learn to use the debugger. Step through the code. THINK ABOUT WHAT YOU ARE DOING.

I tried running the code. I do see something that looks very strange in what you are doing. On the first pass through the loop, I see this get dumped out:

M_p =
       11.029
M_0 =
       12.132

Then on the third pass through, I see this get dumped out:

M_p =
       10.677            0       10.677
M_0 =
       11.744            0       11.744

At the end, we have:

M_0 =
       23.211            0            0            0            0            0            0            0       23.211

So M_p and M_0 are growing in size. Is that the goal here? We are given no clue as to what you expect.

Yet I don't see any line where you even assign some result into a vector element inside the loop. You are treating these variables as if they are all scalars inside the loop. You use the * and / operators, as if the variables are all scalars.

That in turn suggests that you probably are misusing vectors in some place. I look to see how these variables are created:

M_p=(R_max-1)*(M_s+M_L);
M_0=M_s+M_p+M_L 

So they are still treated as scalars. How strange.

But then, at the VERY end of your loop, I see this block of code:

R_out(1,loop)=R_max
Weq_out(1,loop)=W_eq
t_b(1,loop)=t_bmax
p_0(1,loop)=P_0_a
delta_ratio(1,loop)=delta
D_out=D_max
L_out=L_max
X_CG(1,loop)=X_cg
X_CP(1,loop)=X_cp
M_p(1,loop)=M_p
M_0(1,loop)=M_0
epsilon_c(1,loop)= epsilon
M_s(1,loop)=M_s

SIGH. SERIOUSLY?

Go back up, just to check that before the loop, you created VECTORS M_p and M_0.

M_p=zeros(1,9);
M_0=zeros(1,9);

But then we return to the code, and inside the main part of the loop, remember you treat everything as if it is a SCALAR variable.

M_p=(R_max-1)*(M_s+M_L);
M_0=M_s+M_p+M_L 

Do you see the problem?

When you do this:

M_p=(R_max-1)*(M_s+M_L);

MATLAB OVERWRITES THE VECTOR M_p. Yet, you are also using the variable M_p to store the result of this entire operation, for each pass through the loop. The same thing is done with each of those vectors.