I have an invalid use of operator error

function [T p rho] = atmos(h)

h1 = 17;    % Height of tropopause
h2 = 20;    % End height of table
g = 9.81;
R = 287; 
c = 6.51e-3;       % temperature lapse dt/dh = - c = -6.51 degcelcius/km

T0 = 30+273.15;       % Temperature sea level
p0 = 101325;       % pressure sealevel
rho0 = 101325/R/T0;      % density sealevel = pressure / R*T, R=287, T = 30 degcelcius       

T1 = T0 - c*h1*1000;    % Temperature at 17km
p1 = p0 * (T1/T0)^5.2506;   % Pressure at 17km
rho1 = rho0 * (T1/T0)^4.2506;   % Density at 17km
T2 = T1;    % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1)*1000);  % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1)*1000);  % Density at 20km

if h <= h1
    disp('Troposphere');
    T = T0 - c*h*1000;
    p = p0 * (T/T0)^5.2506;
    rho = rho0 * (T/T0)^4.2506;
elseif h <= h2
    disp('Tropopause');
    T = T1;
    p = p1 * exp(-g/(R*T)*(h-h1)*1000);
    rho = rho1 * exp(-g/(R*T)*(h-h1)*1000);
end

h = 0:0.2:20;
for i=1:size(h,2)
    [T(i), p(i), rho(i)] = atmos(h(i));
end

figure;
plot(T-273.15, h);
xlabel('Temperature (^{o}C)');
ylabel('Altitude (km)');
grid on;

figure;
plot (p,h);
xlabel('Pressure (N/m^2)');
ylabel('Altitude (km)');
grid on;

figure;
plot (rho,h);
xlabel('Density (kg/m^3)');
ylabel('Altitude (km)');
grid on;

end

I dont know what went wrong but there is a problem with the if h <= h1 part

[KSSV]

What wrong you have? Does it throw any error? It is running fine in my pc.

[per isakson]

I fail to reproduce the problem with if h <= h1 . How did you call atmos() ?

However, I get

>> atmos(12)
...
Troposphere
Troposphere
Out of memory. The likely cause is an infinite recursion within
the program.
Error in atmos (line 34)
[T(i), p(i), rho(i)] = atmos(h(i)); 

Answer 1

[per isakson]

I have modified your function to avoid the recursive call of the function.

atmos

ans = 101×1

303.1500 301.8480 300.5460 299.2440 297.9420 296.6400 295.3380 294.0360 292.7340 291.4320

function [T,p,rho] = atmos()
h = 0:0.2:20;
len = size(h,2);
T = nan(len,1);
p = nan(len,1);
rho = nan(len,1);
for i=1:len
    [T(i), p(i), rho(i)] = atmos_(h(i));
end

figure;
plot(T-273.15, h);
xlabel('Temperature (^{o}C)');
ylabel('Altitude (km)');
grid on;

figure;
plot (p,h);
xlabel('Pressure (N/m^2)');
ylabel('Altitude (km)');
grid on;

figure;
plot (rho,h);
xlabel('Density (kg/m^3)');
ylabel('Altitude (km)');
grid on;

end
function [T,p,rho] = atmos_(h)

h1 = 17;    % Height of tropopause
h2 = 20;    % End height of table
g = 9.81;
R = 287;
c = 6.51e-3;       % temperature lapse dt/dh = - c = -6.51 degcelcius/km

T0 = 30+273.15;       % Temperature sea level
p0 = 101325;       % pressure sealevel
rho0 = 101325/R/T0;      % density sealevel = pressure / R*T, R=287, T = 30 degcelcius

T1 = T0 - c*h1*1000;    % Temperature at 17km
p1 = p0 * (T1/T0)^5.2506;   % Pressure at 17km
rho1 = rho0 * (T1/T0)^4.2506;   % Density at 17km
T2 = T1;    % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1)*1000);      %#ok<NASGU> % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1)*1000);  %#ok<NASGU> % Density at 20km

if h <= h1
    % disp('Troposphere');
    T = T0 - c*h*1000;
    p = p0 * (T/T0)^5.2506;
    rho = rho0 * (T/T0)^4.2506;
elseif h <= h2
    % disp('Tropopause');
    T = T1;
    p = p1 * exp(-g/(R*T)*(h-h1)*1000);
    rho = rho1 * exp(-g/(R*T)*(h-h1)*1000);
end
end

[Rachel Ong]

Hello, thank you for your help. I'm a bit new to this so can I ask what this part means?

len = size(h,2);
T = nan(len,1);
p = nan(len,1);
rho = nan(len,1);

[Walter Roberson]

size(h,2) is the number of columns in the variable h. You then store the number of columns into the variable named len

nan(len,1) means to created an array of Not A Value entries that has len rows and 1 columns. len was defined as the number of columns in h, so the result will be an array that has as many rows as h had columns, and each entry has the value NaN (Not A Number)

[per isakson]

These lines allocates memory for the vectors, T, p and rho. With large arrays it's important for performance. I use nan, because if nan is not over-written in the loop this mostly becomes obvious. See Preallocation.

[Rachel Ong]

I see, thank you for your help! :-)