# I have an invalid use of operator error

1 view (last 30 days)
Rachel Ong on 19 Oct 2021
Commented: Rachel Ong on 19 Oct 2021
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
##### 2 CommentsShowHide 1 older comment
per isakson on 19 Oct 2021
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));

per isakson on 19 Oct 2021
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 on 19 Oct 2021
I see, thank you for your help! :-)