I have an invalid use of operator error

2 visualizzazioni (ultimi 30 giorni)
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 Commenti
KSSV
KSSV il 19 Ott 2021
What wrong you have? Does it throw any error? It is running fine in my pc.
per isakson
per isakson il 19 Ott 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));

Accedi per commentare.

Risposta accettata

per isakson
per isakson il 19 Ott 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
  4 Commenti
per isakson
per isakson il 19 Ott 2021
Modificato: per isakson il 19 Ott 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.
Rachel Ong
Rachel Ong il 19 Ott 2021
I see, thank you for your help! :-)

Accedi per commentare.

Più risposte (0)

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by