For loop gives the same output

3 views (last 30 days)
I have a problem with this code that it gives the same output when the loop reaches the 4th column of the inputs, the loop gives the same results (when it should not). For the hell of me I don't know what is wrong! It has to be in the loop not in the numbers! I ran them indvidually for all scenarios and got unique output for every scenario!
%% Defining Scenario Independent Paramters:
%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
%%My 9 scenarios that are defined the text data are:
%Maximum Altitude, h_max= 0.3048*[20000 20000 20000 2000 2000 10000 30000 10000 30000]; %m
%Normalized max acceleration, a_max= [10 10 10 5 20 10 10 5 20]
% Static margin , SM = [1 2 3 2 2 2 2 1 3]
%% Calculating D_max and L_max through the 9 secenarios
Data=readmatrix('Data');% Each col is 1 test
D_max=zeros(1,9);
L_max=zeros(1,9);
for loop=1:9
h_max=Data(1,loop);
a_max=Data(2,loop);
SM=Data(3,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;
%Iteration for best L, D
% Creating Length and Diameter Limiatation for Iteration
L=linspace(0.0001,4,5000); %m
D=linspace(0.0001,2,5000); %m
%Iteration through all the parameters:
for i=1:length(D);
for j=1:length(L)
delta= (P_c/(2*sig_s))*D(i); % thickness of shell m
M_n=delta*rho_s*pi*D(i)*(D(i)+sqrt(D(i)^2+(D(i)^2/4)));
M_f=((D(i)^2)/2)*delta*rho_s;
M_f_b=(pi*D(i)*rho_s*D(i)*delta);
M_s=(pi*D(i)*rho_s*L(j)*delta)+M_n+M_f+M_f_b;%%%%%%%%%%
M_p=(R_max-1)*(M_s+M_L);
L_p=(M_p/(pi*D(i)^2*rho_p/4));
if L_p<L(j)+D(i)
% Center of Pressure for Rocket Nose
X_n=(2/3)*D(i);%m
CN_n=2;
%Center of Pressure for Rocket Fin
a=D(i);%m
s=D(i);%m
b=0;
m=a-b;
fin_hyp=sqrt(2)*D(i);%m
X_f=D(i)+L(j);%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(i))^2)/(1+sqrt(1+((2*fin_hyp)/(a+b))^2));%m
k_fb=1+((D(i)/2)/(s+(D(i)/2)));%m
CN_fb=k_fb*CN_f; %m;
X_cp=((CN_n*X_n)+(CN_fb*X_f_dash))/(CN_n+CN_fb);
%Center of Gravity for Rocket Cone:
X_ncg=2*D(i)/3;
M_n=delta*rho_s*pi*D(i)*(D(i)+sqrt(D(i)^2+(D(i)^2/4)));
M_c=M_L+M_n;
%Center of Gravity for Rocket Fin
X_f_cg=(2*D(i)/3)+L(j)+D(i);
M_f=((D(i)^2)/2)*delta*rho_s;
%Center of Gravity for Rocket Tube 1
L_1=L(j)+D(i)-L_p;
xL_1_cg=(L_1/2)+D(i);
M_L_1=pi*D(i)*L_1*delta*rho_s;
%Center of Gravity for Rocket Tube 2
L_2=L_p;
xL_2_cg=(L_2/2)+D(i)+L_1;
M_L_2=(pi*D(i)*L_2*delta*rho_s)+M_p;
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));
cond(i,j)=X_cp-X_cg-(D(i)*SM);
if cond(i,j)<0 || cond(i,j)>0.00001
continue
end
L_D(i,j)=L(j)/D(i);
lamda(i,j)=M_L/(M_s+M_p);
else
continue
end
end
end
%Output
% Matrix indexing for the D_max and L_max
[M,I] = max(lamda,[],'all','linear');
[row,col] = ind2sub(size(lamda),I);
D_max(1,loop)=D(row)
L_max(1,loop)=L(col)
end

Accepted Answer

Bjarke Skogstad Larsen
Bjarke Skogstad Larsen on 12 May 2020
Edited: Bjarke Skogstad Larsen on 12 May 2020
You need to initialize the matrixes: L_D, lamda and cond.
Always initialize variables before a loop when you only set part of the variable at a time inside the loop, like
lamda(i,j)=M_L/(M_s+M_p);
needs to be initialized like so:
lamda = nan(length(D),length(L));
for ii=1:length(D)
I will include the full working code below (I changed i to ii and j to jj as I don't like using these variables in Matlab due to the possibility of complex values):
clc
clear
%% Defining Scenario Independent Paramters:
%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
%%My 9 scenarios that are defined the text data are:
%Maximum Altitude, h_max= 0.3048*[20000 20000 20000 2000 2000 10000 30000 10000 30000]; %m
%Normalized max acceleration, a_max= [10 10 10 5 20 10 10 5 20]
% Static margin , SM = [1 2 3 2 2 2 2 1 3]
%% Calculating D_max and L_max through the 9 secenarios
Data=readmatrix('Data');% Each col is 1 test
D_max=zeros(1,9);
L_max=zeros(1,9);
for loop=1:9
disp(num2str(loop));
h_max=Data(1,loop);
a_max=Data(2,loop);
SM=Data(3,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;
%Iteration for best L, D
% Creating Length and Diameter Limiatation for Iteration
L=linspace(0.0001,4,5000); %m
D=linspace(0.0001,2,5000); %m
%Iteration through all the parameters:
L_D = nan(length(D),length(L));
lamda = nan(length(D),length(L));
cond = nan(length(D),length(L));
for ii=1:length(D)
for jj=1:length(L)
delta= (P_c/(2*sig_s))*D(ii); % thickness of shell m
M_n=delta*rho_s*pi*D(ii)*(D(ii)+sqrt(D(ii)^2+(D(ii)^2/4)));
M_f=((D(ii)^2)/2)*delta*rho_s;
M_f_b=(pi*D(ii)*rho_s*D(ii)*delta);
M_s=(pi*D(ii)*rho_s*L(jj)*delta)+M_n+M_f+M_f_b;%%%%%%%%%%
M_p=(R_max-1)*(M_s+M_L);
L_p=(M_p/(pi*D(ii)^2*rho_p/4));
if L_p<L(jj)+D(ii)
% Center of Pressure for Rocket Nose
X_n=(2/3)*D(ii);%m
CN_n=2;
%Center of Pressure for Rocket Fin
a=D(ii);%m
s=D(ii);%m
b=0;
m=a-b;
fin_hyp=sqrt(2)*D(ii);%m
X_f=D(ii)+L(jj);%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(ii))^2)/(1+sqrt(1+((2*fin_hyp)/(a+b))^2));%m
k_fb=1+((D(ii)/2)/(s+(D(ii)/2)));%m
CN_fb=k_fb*CN_f; %m;
X_cp=((CN_n*X_n)+(CN_fb*X_f_dash))/(CN_n+CN_fb);
%Center of Gravity for Rocket Cone:
X_ncg=2*D(ii)/3;
M_n=delta*rho_s*pi*D(ii)*(D(ii)+sqrt(D(ii)^2+(D(ii)^2/4)));
M_c=M_L+M_n;
%Center of Gravity for Rocket Fin
X_f_cg=(2*D(ii)/3)+L(jj)+D(ii);
M_f=((D(ii)^2)/2)*delta*rho_s;
%Center of Gravity for Rocket Tube 1
L_1=L(jj)+D(ii)-L_p;
xL_1_cg=(L_1/2)+D(ii);
M_L_1=pi*D(ii)*L_1*delta*rho_s;
%Center of Gravity for Rocket Tube 2
L_2=L_p;
xL_2_cg=(L_2/2)+D(ii)+L_1;
M_L_2=(pi*D(ii)*L_2*delta*rho_s)+M_p;
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));
cond(ii,jj)=X_cp-X_cg-(D(ii)*SM);
if cond(ii,jj)<0 || cond(ii,jj)>0.00001
continue
end
L_D(ii,jj)=L(jj)/D(ii);
lamda(ii,jj)=M_L/(M_s+M_p);
else
continue
end
end
end
%Output
% Matrix indexing for the D_max and L_max
[M,I] = max(lamda,[],'all','linear');
[row,col] = ind2sub(size(lamda),I);
D_max(1,loop)=D(row)
L_max(1,loop)=L(col)
end

More Answers (0)

Categories

Find more on Stress and Strain in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!