# I have an invalid use of operator error

2 views (last 30 days)

Show older comments

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 Comments

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));

### Accepted Answer

per isakson
on 19 Oct 2021

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

atmos

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

##### 4 Comments

per isakson
on 19 Oct 2021

Edited: per isakson
on 19 Oct 2021

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!