//http://www.basicairdata.eu/calculation-routines.html
//Barometric altitude calculation according to US standard atmosphere 1962/72 and ISA
//(c) JLJ www.basicairdata.eu
//Standard temperature 288.15 °K and 101325 Pa
//Height < 11000 m
clear;
graf=0; //0,1 If set to 1 a plot is shown
//ISA atmosphere
p0isa=101325 //Pa ISA standard
T0=288.15//°K ISA STANDARD
g=9.80665;
R=8.314462175//J/K/mol
Mwa=28.9644;//kg/1000/mol
Rair=R*1000/Mwa;
if graf==0 then
Ps=input('Input static pressure reading Pa : ')
Ps=Ps/133.3223684211/25.4//Pa to inHg Conversion
//Using Goodrich 4081 Air data handbook formula 
h=(29.92126^0.190255-Ps^0.190255)/0.000013125214; //US atmosphere 1962
//Back to SI
h=h*0.3048;
printf("Pressure height %.2f m \n",h);
pisa=p0isa*(1-0.0065*h/T0)^(g/Rair/0.0065) //ISA definition
printf("Pressure inverse reading according ISA %.2f Pa\n",pisa);
i=1;
end
//Standard atmosphere calculation
altitude=linspace(0,11000,11000);//Altitude values from 0 m to 11000m 
Tair=T0-altitude*0.0065; //Calculate the temperature at each altitude
//Perfect gas law, dry air
//pV=n*R*T=m/Mw*R*T=Rair*T  ->p=rhoair*Rair*T  rhoair=p/Rair/T
//p(altitude)
pisap=p0isa*(1-0.0065*altitude./T0)^(g/Rair/0.0065) //ISA definition
//rhoair(altitude)
rhoair= pisap./Tair./Rair //kg/m^3
//For relative density rhoair(altitude)/rhobase*100
rhoairrel=rhoair./1.225*100;
i=1;
for  temp=Tair
viscosity(i)=18.27e-3*(291.15+120)/(temp+120)*(temp/291.15)^(3/2)
i=i+1;
end

//Plot section
if graf==1 then
subplot(211)
plot(rhoair,altitude);
a=gca(); 
a.title.font_size=4;
a.x_label.font_size=3;
a.y_label.font_size=3;
a.labels_font_size=3;
pollo=a.children.children(1)
pollo.thickness=4;
xtitle('Density profile with altitude')
xgrid
xlabel('Density kg/m3')
ylabel('Altitude m')
subplot(212)
plot(rhoairrel,altitude);
xgrid
a=gca(); 
a.x_label.font_size=3;
a.y_label.font_size=3;
a.labels_font_size=3;
pollo=a.children.children(1)
pollo.thickness=4;
xlabel('Percent of Standard sea level density')
ylabel('Altitude m')
//new figure
scf();
a=gca();
xgrid
plot(viscosity,altitude);
a.title.font_size=4;
a.x_label.font_size=3;
a.y_label.font_size=3;
a.labels_font_size=3;
pollo=a.children.children(1)
pollo.thickness=4;
xtitle('Viscosity profile with altitude')
xlabel('Viscosity mPas or Cp')
ylabel('Altitude m')
end

//Function definition for pressure correction example
function F=isapressure(x)
  F=x-p0isa*(1-0.0065*h/T0)^(g/Rair/0.0065);
endfunction

pressure_reading=101800;
p0isa=101800;
function F=isaaltitude(x)
  F=pressure_reading-p0isa*(1-0.0065*x/T0)^(g/Rair/0.0065);
endfunction