Calculating Air Density from altitude and then storing in a matrix

Hi, everyone. Basically I have a variable, h (the altitude) which is changing and is thus changing the air density, rho. For altitudes from 0 to 65000 ft, I want to calculate rho and then store the values of rho in a vector. Next, I want to use that vector to plot a graph of W/S and W/P. The details can be seen in my code below:

for h =1:65000
  if h > 36152
  T=-70;
  p=473.1 * exp(1.73 - .000048*h);
  rho(h)=p/ (1718 * (T+459.7));
  elseif h < 36152
  T=59 - .00356*h;
  p=2116* ( (T+459.7)/518.6 )^5.256;
  rho(h)=p/ (1718 * (T+459.7));
  end   
end

%

%
sigma= rho./0.00237717;
RC= 300;
x_10= 16.789; 
eta_p= .68;
for i = 1:length(sigma)
  % density equation
  ws = linspace(1,4,length(sigma));
  wp(i) = eta_p./((RC./33000)+(ws(i)^0.5./(19.*x_10.*sigma(i)^.5)));
end
% Graph W/S and W/P
figure,
plot(ws,wp,'r-')
xlabel('W/S')
ylabel('W/P')

When I run this code, my graph comes out with a random vertical line. I'm a MATLAB beginner, so I'm sure I'm making some dumb mistake (maybe its a mathematical error, though I doubt it, and my code seems fine?), so I would really appreciate it if someone could point why my code isn't running the way I intend, or maybe explain a more efficient method that I could use.

Answer 1

[KSSV]

You need to modify your if else condition.....

Repalce:

   if h >36152

with

   if h >= 36152

actually when h=36152....value of rho is zero..so the problem. And also you have to initialize the values which you are filling in loop.

N = 65000 ;
rho = zeros(1,N) ;
for h =1:N
    if h >= 36152
        T=-70;
        p=473.1 * exp(1.73 - .000048*h);
        rho(h)=p/ (1718 * (T+459.7));
    elseif h < 36152
        T=59 - .00356*h;
        p=2116* ( (T+459.7)/518.6 )^5.256;
        rho(h)=p/ (1718 * (T+459.7));
    end
end
%
%
sigma= rho./0.00237717;
RC= 300;
x_10= 16.789;
eta_p= .68;
wp = zeros(1,N) ;
for i = 1:length(sigma)
    % density equation
    ws = linspace(1,4,length(sigma));
    wp(i) = eta_p./((RC./33000)+(ws(i)^0.5./(19.*x_10.*sigma(i)^.5)));
end
% Graph W/S and W/P
figure,
plot(ws,wp,'r-')
xlabel('W/S')
ylabel('W/P')

Note that above code can be vctorised. As you are a beginner, it is better to start with loops first and then vectorise.