I have an invalid use of operator error
1 view (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
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.
More Answers (0)
See Also
Categories
Find more on String Parsing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!