Graph Not Plotting for the For Loop Code

I'm trying to plot a graph of air density against altitude using the following density altitude eq on the NASA website:

https://www.grc.nasa.gov/www/k-12/airplane/atmosmet.html

This is the code I've written:

global rho


for h = 0:10:40000
    
    if (h < 11000)
        T = 25.05 - h.*0.00649;
        p = 101.29* ((T + 273.15)./288.08)^5.256;
        
        
    else
        if (11000 < h) & (h< 250000)
            T = -56.46;
            p = 22.65* exp(1.73 - h.*0.000157);
            
        else
            if (h > 25000)
                T = -131.21 + h.*0.00299;
                p = 2.488*((T + 273.15)/216.6)^-11.3888;
            end
        end
    end
    
    rho = p / (0.2869*(T + 273.15));
    
end

plot (h,rho)
xlabel('height')
ylabel('Density')

But nothing is showing up on the graph plot. What am I missing?

Answer 1

[Stephen23]

"What am I missing? "

Assigning the values to an output array via indexing.

Because you did not use indexing to assign the values to an array your code simply overwrites the previous iteration's values until the final iteration, leaving you with one data value (which is not visible when plotted as the default line has no marker).

This works and plots all of the loop results (only you can check if it is correct):

G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = 21*pi;                 % Area of the rocket 
Cd = 0.4;               % Drag Coefficient for the rocket
t = 0:1:4000;
h_v = 0:100:5000;
rho = nan(size(h_v)); % <--- preallocate the output array

for k = 1:numel(h_v) % <--- loop over indices, not data values
    h = h_v(k);
    if (h < 11000)
        T = 25.05 - h.*0.00649;
        p = 101.29* ((T + 273.15)./288.08)^5.256;
    else
        if (11000 < h) & (h< 250000)
            T = -56.46;
            p = 22.65* exp(1.73 - h.*0.000157);
        else
            if (h > 25000)
                T = -131.21 + h.*0.00299;
                p = 2.488*((T + 273.15)/216.6)^-11.3888;
            end
        end
    end
    rho(k) = p / (0.2869*(T + 273.15));
    %   ^ use indexing to assign output to an array
end
plot (h_v,rho)
xlabel('height')
ylabel('Density')

Logical indexing may be a simpler approach than using a loop.

[Asit Rahman]

This worked, thank you so much!

[Torsten]

I think you will have to replace

if (11000 < h) & (h< 250000)

by

if (11000 < h) & (h< 25000)

Answer 2

[Torsten]

H=linspace(0,40000,4001)
rho=zeros(1,numel(H))
for i=1:numel(H)
    h = H(i);
    if h > 25000
        T = -131.21 + h.*0.00299;
        p = 2.488*((T + 273.15)/216.6)^-11.3888;
    elseif  h<= 25000 && h>11000 
        T = -56.46;
        p = 22.65* exp(1.73 - h.*0.000157);
    else
        T = 25.05 - h.*0.00649;
        p = 101.29* ((T + 273.15)./288.08)^5.256;
    end
    rho(i) = p / (0.2869*(T + 273.15));
end
plot(H,rho)

[Asit Rahman]

This also works, thank you!

Answer 3

[Arif Hoq]

your output is a single value. if your code is going well then try this with a marker

plot (h,rho, 'o')

[Arif Hoq]

G = 6.67408e-11;        % Universal Gravitational Constant
M = 5.9722e24;          % Mass of the earth
R = 6371e3;             % Radius of the earth(m)
A = 21*pi;                 % Area of the rocket 
Cd = 0.4;               % Drag Coefficient for the rocket
t = 0:1:4000;
%h = 9.8*t.^2;
%h = 0:100:5000;
for h = 0:10:4000
    
    if (h < 11000)
    T = 25.05 - h.*0.00649;
    p = 101.29* ((T + 273.15)./288.08)^5.256;
    
   
    else
        if (11000 < h) & (h< 250000)
        T = -56.46;
        p = 22.65* exp(1.73 - h.*0.000157);
        
        else
            if (h > 25000)
            T = -131.21 + h.*0.00299;
            p = 2.488*((T + 273.15)/216.6)^-11.3888;
            end
        end
    end
           
rho = p / (0.2869*(T + 273.15));
end
plot (h,rho, 'o')
xlabel('height')
ylabel('Density')

[Asit Rahman]

It's not supposed to be a single point, it's supposed to plot how density changeswith altitude. That's what I put in the for loop to check the density in iterations of 10m