For loop gives the same output

조회 수: 6(최근 30일)
Mohsen Nashel
Mohsen Nashel 2020년 5월 12일
댓글: Mohsen Nashel 2020년 5월 13일
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

채택된 답변

Bjarke Skogstad Larsen
Bjarke Skogstad Larsen 2020년 5월 12일
편집: Bjarke Skogstad Larsen 2020년 5월 12일
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
  댓글 수: 1
Mohsen Nashel
Mohsen Nashel 2020년 5월 13일
Thank you so much! It worked!

댓글을 달려면 로그인하십시오.

추가 답변(0개)

태그

제품

Community Treasure Hunt

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

Start Hunting!

Translated by