PLEASE HELP!!
So i have an assignment i need to do for school, and i have hit an intriguing obstacle. It happens to be the first matlab obstacle i havent been able to surpass.
INFORMATION-So, when the you look at the earth's atmosphere under standard temperature and pressure, you can find out the values and ratios of temperature, pressure, and density values with a lone input of altitude.
Pressure and density are determined by calculating temperature which is calculated by altitude. Temperature goes through a series of gradient (sloped temperature change with altitude) and isothermal (no temperature change with altitude) sections.
MY PROBLEM-When switching from a gradient section (0m<altitude<11000m) to an isothermal section (11000<=Altitude<25000) within an group of if-statements, my matlab presents an error of "??? Reference to a cleared variable Prat." where Prat (Pressure Ratio) is a calculated function in "if 0m<altitude<11000" and "elseif (11000<=Altitude<25000)"
OTHER TESTS- the program runs flawlessly when the altitude stays within one single if/elseif statement, but any transition between altitude ranges presents this error. So, i'm pretty confident that there isn't a simple error in coding.
QUESTION- how do i transition between altitude ranges?? for example, entering altitude=[0:10000] or altitude=[11000:24000] works fine, but entering altitude=[10999:11000] wont work
Code:%Goal: Calculate Pressure, Pressure Ratio, Density, Density Ratio, %Temperature, and Temperature Ratio based on an Altitude input. %This program operates using air as a medium. The Gas Constant, density, %and atmospheric pressure are constant and refer to standard atmospheric temperature %and pressure. Gravitational pull is also considered constant. clear,clc alt=input('Enter your altitude: '); go=9.8; %m/s^2 R=287.05; %J/K mol P1=101325; %N/,^2 rho1=1.225; %kg/m^3 T1=288.16; %K h1=0; a=[-6.5*10^-3 3*10^-3 -4.5*10^-3 4*10^3]; aeng=a.*.54864; i=1; while i<=length(alt) if alt==0 T(i)=288.6 %K Prat(i)=1 rhorat(i)=1 elseif alt>=0 & alt<11000 a=a(1); %K/m T(i)=T1+a.*(alt(i)-h1) Prat(i)=(T(i)./T1).^-(go./(a.*R)) rhorat(i)=exp(1).^-(go./(a.*T(i))+1) elseif alt>=11000 & alt<25000 T(i)=216.66 %K Prat(i)=exp(1).^-(go./(R.*T(i)).*(alt(i)-h1)) rhorat(i)=exp(1).^-(go./(R.*T(i)).*(alt(i)-h1)) elseif alt>=25000 & alt<47000 a=a(2); %K/m T(i)=T1+a.*(alt(i)-h1) Prat(i)=(T(i)./T1).^-(go./(a.*R)) rhorat(i)=exp(1).^-(go./(a.*T(i))+1) elseif alt>=47000 &alt<53000 T(i)=282.66 %K Prat(i)=exp(1).^-(go./(R.*T(i)).*(alt(i)-h1)) rhorat(i)=Prat(i+1) elseif alt>=53000 & alt<79000 a=a(3); %K/m T(i)=T1+a.*(alt(i+1)-h1) Prat(i)=(T(i)./T1).^-(go./(a.*R)) rhorat(i)=exp(1).^-(go./(a.*T(i))+1) elseif alt>=79000 & alt<90000 T(i)=165.66 %K Prat(i)=exp(1).^-(go./(R.*T(i)).*(alt(i)-h1)) rhorat(i)=Prat(i) elseif alt>=90000 & alt<105000 a=a(1); %K/m T(i)=T1+a.*(alt(i)-h1) Prat(i)=(T(i)./T1).^-(go./(a.*R)) rhorat(i)=exp(1).^-(go./(a.*T(i))+1) elseif alt==105000 T(i)=225.66; %K Prat(i)=exp(1).^-(go./(R.*T(i)).*(alt(i)-h1)) rhorat(i)=Prat(i+1) elseif alt<0 | alt>105000 fprintf('Error: This code only evaluates outputs within a range of 0-105000 meters\n') end P=P1.*Prat rho=rho1.*rhorat sigma=rho./rho1 i=i+1; end
If alt=[100,2000] what do you think
does?Code:if alt==0 x=x+1 elseif (alt>0) & (alt<1000) x=x+2 end
Also the line:
rhorat(i)=Prat(i+1)
appears to use Prat(i+1) before its calculated.
Other comments: These are not show stoppers, but they are good practice:
1. Do not use a while loop when you know the trip count around a loop. use a for loop instead.
2. Declare arrays of known dimensions, do not let Matlab rearrange memory every time you access a different element.
I guess I've been using some offbeat logic on this code then, I guess
I was under the impression I could use two separate equations for one variable so long as they were indexed correctly (which, somehow my matlab wasn't doing outside of a loop in an earlier version of my code [obviously due to some other error])
I assume you are implying that my if/elseif setup will just cause my Prat from the first elseif to try and overwrite the Prat from the second elseif? Sorry I need the clarification; it seems logical it just want intuitive for me :/
also, good call on the Prat(i+1)--it was an earlier (and incorrect) indexing I was using before.... I guess I'm tripping up all over the place :/. I hope I'm just rusty
what advantages does a for loop have over a while loop?
And, if you don't mind, do you know of any way I could successfully set up what I was trying to do? Since the idea is to plot the outputs at the end, I suppose I could break up the different pressure and temperature sections into different equations and stick them together at the end
Thanks a lot for the advice so far!! I really appreciate it!
When you apply a logic test to a vector or array you get a vector of results. If it is in an if structure it appears to be interpreted as:
if (any(alt==0))
You should not be applying a logical test to a vector in an if statement, you should be applying it to a specific element like:
if alt(i)==0
CB
You, sir, are a gentleman and a scholar.
I can't express how awesome i felt when the code started working right.
I completely understand where you're coming from now, so thanks a ton for the information. It broadened my insight.
Anyway, there is still one roadblock I'm hitting, that I can't seem to figure out. Its an indexing error whenever i run between 25000 and 47000 meters. I've seen the error before, but i cant seem to get rid of it this time :/
"??? Index exceeds matrix dimensions"
Here's my code as of now (and no, i didn't ignore you before; but for the interest of time i haven't changed it to a for loop yet)
%Goal: Calculate Pressure, Pressure Ratio, Density, Density Ratio,
%Temperature, and Temperature Ratio based on an Altitude input.
%This program operates using air as a medium. The Gas Constant, density,
%and atmospheric pressure are constant and refer to standard atmospheric temperature
%and pressure. Gravitational pull is also considered constant.
clear,clc
alt=input('Enter your altitude: ');
go=9.8; %m/s^2
R=287.05; %J/K mol
P1=101325; %N/,^2
rho1=1.225; %kg/m^3
T1=288.16; %K
h1=0;
a=[-6.5*10^-3 3*10^-3 -4.5*10^-3 4*10^3];
aeng=a.*.54864;
ao=0;
i=1;
while i<=length(alt)
if alt(i)==0
T(i)=288.6 %K
Prat(i)=1
rhorat(i)=1
elseif alt(i)>=0 & alt(i)<11000
a=a(1); %K/m
T(i)=T1+a.*(alt(i)-h1)
Prat(i)=(T(i)./T1).^-(go./(a.*R))
rhorat(i)=exp(1).^-(go./(a.*T(i))+1)
elseif alt(i)>=11000 & alt(i)<25000
T(i)=216.66+ao.*(alt(i)-h1) %K
Prat(i)=exp(1).^-(go./(R.*T(i)).*(alt(i)-h1))
rhorat(i)=exp(1).^-(go./(R.*T(i)).*(alt(i)-h1))
elseif alt(i)>=25000 & alt(i)<=47000
a=a(2); %K/m
T(i)=T1+a.*(alt(i)-h1)
Prat(i)=(T(i)./T1).^-(go./(a.*R))
rhorat(i)=exp(1).^-(go./(a.*T(i))+1)
end
P(i)=P1.*Prat(i)
rho(i)=rho1.*rhorat(i)
sigma(i)=rho(i)./rho1
i=i+1
end