#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""**ICAO Standard Atmosphere (ISA) - functions and routines**
| © Dietmar Thaler 2014-2016
| E-Mail: info@foehnwall.at
| Web : http://www.foehnwall.at/pythonmodules.html
.. warning:: **adisalib** calculates ISA only up to 32000 gpm!
.. note:: This module makes use of the **numpy** library. NumPy is the
fundamental package for scientific computing with Python. Make sure
you have it installed. See http://www.numpy.org/ and
http://sourceforge.net/projects/numpy/
Both **Python 2.7** and **Python >= 3.4** compatibility assured.
The core routines of the library are:
* isa_pressure(height)
* isa_density(height)
* isa_temperature(height)
* isa_height(pressure)
* isa_QNH(barometer_altitude, pressure)
* isa_QFE(barometer_altitude, arp_altitude, pressure)
* isa_densityaltitude(density)
* isa_table(h1,h2,dh)
Added are a couple of constants, variables and conversion routines
(see below). A simple demo program that produces a matplotlib-graphic
can be found at
`www.foehnwall.at/isa.html. <http://www.foehnwall.at/isa.html>`_
--------------------------------------------------------------
*This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.*
*This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.*
*You should have received a copy of the GNU General Public License.
If not, see* <http://www.gnu.org/licenses/gpl.html> .
--------------------------------------------------------------
"""
# Make python2 compatible to python3 and vice versa
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals # maybe problems with module csv
import sys
err_NUMPY = """ERROR: numpy missing!
This module requires the **numpy** library. NumPy is the fundamental
package for scientific computing with Python. Make sure you have it
installed. See http://www.numpy.org/ and
http://sourceforge.net/projects/numpy/"""
try:
import numpy as np
except:
print(err_NUMPY)
sys.exit(-3)
__version__ = "0.2.0-2016-12-20" # : module version
__copyright__ = "D. Thaler"
_FOOT2METER_ = 0.3048 # m
Tice = 273.15 # : freezing point water under standard conditions in K
isa_g0 = 9.80665 # : ISA standard gravity m/s^2
isa_M0 = 28.9664 # : ISA molar mass (kg / kmol)
isa_R = 8314.32 # : ISA universal ideal gas constant [J/(kmol K)]
isa_Rl = isa_R / isa_M0 # : ISA gas constant of dry air J/(kg K)
isa_p0 = 1013.25 # : ISA surface pressur hPa
isa_T0 = 288.15 # : ISA surface temperature (at z=0)
isa_rho0 = isa_p0 * 100. / isa_T0 / isa_Rl # : ISA surface density kg/m^3
isa_GAMMA0 = 0.0065 # : ISA temperature gradient K/m from h=0 to h=11000 gpm
isa_h1 = 11000. # : ISA tropopause geopot. altitude
isa_p1 = isa_p0 * ((1.0 - isa_h1 / isa_T0 * isa_GAMMA0) ** (isa_g0 / isa_GAMMA0 / isa_Rl)) # : ISA pressure at 11000 gpm (tropopause)
isa_T1 = isa_T0 - isa_GAMMA0 * isa_h1 # : ISA temperature between h = 11000 and h = 20000 gpm (tropopause)
isa_rho1 = isa_p1 * 100. / isa_T1 / isa_Rl # : ISA tropopause density kg/m^3
isa_h2 = 20000. # : ISA begin stratosphere altitude in gpm
isa_T2 = isa_T1
isa_p2 = isa_p1 * np.exp(-isa_g0 * (isa_h2 - isa_h1) / isa_Rl / isa_T1) # : ISA pressure at 20000 gpm (begin stratosphere)
isa_rho2 = isa_p2 * 100. / isa_T2 / isa_Rl # : ISA surface density kg/m^3
isa_GAMMA2 = -0.0010 # : ISA temperature gradient K/m from geopotential h=20000 to h=32000 m
isa_h3 = 32000. # : ISA-stratosphere altitude in gpm
isa_p3 = isa_p2 * ((1.0 - (isa_h3 - isa_h2) / isa_T2 * isa_GAMMA2) ** (isa_g0 / isa_GAMMA2 / isa_Rl)) # : ISA pressure at 32000 gpm (begin stratosphere)
isa_T3 = isa_T1 - isa_GAMMA2 * (isa_h3 - isa_h2) # : ISA temperature at h = 32000 gpm
isa_rho3 = isa_p3 * 100. / isa_T3 / isa_Rl # : ISA density kg/m^3 at h = 32000 m
isa_T_error = -999. # : error value for ISA temperature
isa_p_error = -999. # : error value for ISA pressure
isa_rho_error = -999. # : error value for ISA density
isa_h_error = -9999. # : error value for ISA altitude
def __makearray__(obj):
"""Converts a scalar or a list etc ... to a numpy array and returns
a numpy array unchanged"""
if np.isscalar(obj):
return np.array([obj])
else:
return np.array(obj)
[docs]def f2m(foot):
"""Converts foot to meter"""
return foot2meter(foot)
[docs]def m2f(meter):
"""Converts meter to foot"""
return meter2foot(meter)
[docs]def kelvin2celsius(tk):
"""Converts temperature in Kelvin to temperature in Ceslius"""
return tk - Tice
[docs]def c2k(tk):
"""Converts temperature in Celsius to temperature in Kelvin"""
return celsius2kelvin(tk)
[docs]def celsius2kelvin(tc):
"""Converts temperature in Celsius to temperature in Kelvin"""
return tc + Tice
[docs]def k2c(tk):
"""Converts temperature in Kelvin to temperature in Ceslius"""
return kelvin2celsius(tk)
[docs]def isa_temperature(hh):
"""Temperature of ISA-Standardatmosphere as function of altitude
| z ... altitude in gpm
| isa_temperature ... ISA temperature in K
"""
h = __makearray__(hh)
T = np.empty_like(h)
idx = h < isa_h1
T[idx] = isa_T0 - isa_GAMMA0 * h[idx]
idx = h >= isa_h1
T[idx] = isa_T1
idx = h >= isa_h2
T[idx] = isa_T1 - isa_GAMMA2 * (h[idx] - isa_h2)
idx = h > isa_h3
T[idx] = isa_T_error
return T
[docs]def isa_pressure(height):
"""Pressure in ISA below 32000 gpm as function of height
| hh ... altiude in gpm
| isa_pressure .. pressure at level hh in hPa
"""
h = __makearray__(height)
p = np.empty_like(h)
idx = h < isa_h1
p[idx] = isa_p0 * ((1.0 - h[idx] / isa_T0 * isa_GAMMA0) ** (isa_g0 / isa_GAMMA0 / isa_Rl))
idx = h >= isa_h1
dh = h - isa_h1
p[idx] = isa_p1 * np.exp(-isa_g0 * dh[idx] / isa_Rl / isa_T1)
idx = h >= isa_h2
dh = h - isa_h2
p[idx] = isa_p2 * ((1.0 - dh[idx] / isa_T2 * isa_GAMMA2) ** (isa_g0 / isa_GAMMA2 / isa_Rl))
idx = h > isa_h3
p[idx] = isa_p_error
return p
[docs]def isa_height(pressure):
"""Height in ISA below 32000 gpm as function of pressure
| pp ... pressure in hPa
| isa_height .. height in gpm
"""
p = __makearray__(pressure) # make array
h = np.empty_like(p)
idx = p > isa_p1
h[idx] = isa_T0 / isa_GAMMA0 * (1.0 - ((p[idx] / isa_p0) ** (isa_GAMMA0 * isa_Rl / isa_g0)))
idx = p <= isa_p1
h[idx] = isa_h1 - isa_Rl * isa_T1 / isa_g0 * np.log(p[idx] / isa_p1)
idx = p <= isa_p2
h[idx] = isa_h2 + isa_T2 / isa_GAMMA2 * (1.0 - (p[idx] / isa_p2) ** (isa_GAMMA2 * isa_Rl / isa_g0))
idx = p < isa_p3
h[idx] = isa_h_error
return h
[docs]def isa_QNH(barometer_altitude, pressure):
"""Returns QNH in hPa below 3200 gpm aas numpy array
| barometer_altitude in gpm
| pressure ... barometer reading in hPa at barometer altitude as numpy array
"""
hh0 = isa_height(pressure)
dhh = hh0 - barometer_altitude
return isa_pressure(dhh)
[docs]def isa_QFE(barometer_altitude, arp_altitude, pressure):
"""Returns QFE (pressure at aerodrome reference point arp)
as function of:
| barometer_altitude in gpm
| arp_altitude in gpm
| pressure ... pressure reading at baromter_altitude in hPa
"""
dh = barometer_altitude - arp_altitude
h_qnh = isa_height(pressure)
h_arp = h_qnh - dh
return isa_pressure(h_arp)
[docs]def isa_density(hh):
"""Density in ISA below 32000 gpm
| hh ... altitude in gpm
| isa_density ... density kg/m**3 at level hh
"""
h = __makearray__(hh)
rho = np.empty_like(h)
idx = (h <= isa_h1)
rho[idx] = isa_rho0 * ((isa_T0 - isa_GAMMA0 * h[idx]) / isa_T0) ** (isa_g0 / isa_Rl / isa_GAMMA0 - 1.)
idx = (h > isa_h1)
rho[idx] = isa_rho1 * np.exp(-isa_g0 * (h[idx] - isa_h1) / isa_Rl / isa_T1)
idx = (h > isa_h2)
rho[idx] = isa_rho2 * ((isa_T2 - isa_GAMMA2 * (h[idx] - isa_h2)) / isa_T2) ** (isa_g0 / isa_Rl / isa_GAMMA2 - 1.)
idx = (h > isa_h3)
rho[idx] = isa_rho_error
return rho
[docs]def isa_densityaltitude(density):
"""
ISA-density altitude for a given air density
| rho ... actual air density
| isa_densityaltitude ... altitude where ISA-density equals rho (in gpam)"""
lam1 = -1.0 / (1.0 - isa_g0 / isa_Rl / isa_GAMMA0)
lam2 = -1.0 / (1.0 - isa_g0 / isa_Rl / isa_GAMMA2)
rho = __makearray__(density)
h = np.empty_like(rho)
idx = (rho >= isa_density(isa_h1))
h[idx] = isa_T0 / isa_GAMMA0 * (1.0 - (rho[idx] / isa_rho0) ** lam1)
idx = (rho < isa_density(isa_h1))
h[idx] = isa_h1 + isa_Rl * isa_T1 / isa_g0 * np.log(isa_rho1 / rho[idx])
idx = (rho < isa_density(isa_h2))
h[idx] = isa_h2 + isa_T2 / isa_GAMMA2 * (1.0 - (rho[idx] / isa_rho2) ** lam2)
idx = (rho < isa_density(isa_h3))
h[idx] = isa_h_error
return h
[docs]def isa_table(h1=0, h2=32000, dh=500, csv=False):
"""Prints a table with ISA-Data.
| h1 and h2 ... lower and upper vertical boundary in Meter
| dh ... vertical increment in Meter
| csv ... flag (locical)
| True ... the output is a CSV-formatted file with a semicolon ';'
| as field separator and a dot '.' as decimal separator.
| False ... Textfile
"""
FMTDATTXT = "%5.0f %6.2F %7.2f %7.4f"
FMTDATCSV = "%5.0f; %6.2F; %7.2f; %7.4f"
if csv:
fmtdat = FMTDATCSV
else:
fmtdat = FMTDATTXT
h1 = float(h1)
h2 = float(h2)
__ = float(dh)
isa_hh = np.arange(h1, h2 + dh, dh)
isa_temp = isa_temperature(isa_hh)
isa_pppp = isa_pressure(isa_hh)
isa_dens = isa_density(isa_hh)
if csv:
print(" h; T; p; rho")
else:
print(" gpm Celsius hPa kg/m³")
ii = 0
for h in isa_hh:
print(fmtdat % (h, isa_temp[ii] - Tice, isa_pppp[ii], isa_dens[ii]))
ii = ii + 1
return None
def _testadisalib_():
"""A simple testing routine for module adisalib"""
print((_testadisalib_.__doc__))
print(("Version " + __version__))
isa_hhhh = np.array([-300, 100, 0, 500, 1000, 1500, 2000, 2500,
3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000,
11000, 12000, 14000, 16000, 18000, 20000,
24000, 26000, 28000, 30000, 32000]) * 1.0
isa_temp = isa_temperature(isa_hhhh)
isa_pppp = isa_pressure(isa_hhhh)
isa_hhhp = isa_height(isa_pppp)
isa_dens = isa_density(isa_hhhh)
isa_densh = isa_densityaltitude(isa_dens)
print(29 * "=" + " ISA Temp " + 29 * '=')
print()
ii = 0
LSP = "====== ====== ========= ========= ============= ==========="
print(LSP)
print("h[gpm] T[C] h->p[hPa] p->h[gpm] h->rho[kg/m3] rho->h[gpm]")
print(LSP)
fformat = "%5.0f %6.2f %7.2f %7.1f %7.4f %7.1f"
for h in isa_hhhh:
print(fformat % (h, isa_temp[ii] - Tice, isa_pppp[ii], isa_hhhp[ii], isa_dens[ii], isa_densh[ii]))
ii = ii + 1
print(LSP)
def _testqnh_():
LINE0 = 54 * '='
LINE1 = len(LINE0) * '-'
# Testdata
baro_alt = 650.0 # 649.0 # Barometer altitude m
arp_alt = 638.0 # ARP altitude m
# Barometer readings hPa
press = np.array([966.3, 954.2, 955.2, 948.8, 930.9, 924.6])
qnh_baros = isa_QNH(baro_alt, press)
qfes = isa_QFE(baro_alt, arp_alt, press)
qnh_arps = isa_QNH(arp_alt, qfes)
dqfes = qfes - press
print()
print(LINE0)
print('Test QNH and QFE calculations')
print('-----------------------------')
print()
fmtalt = 'Barometer Altitude: {0:6.1f} m ARP altitude: {1:6.1f} m'
fmthead = 'Baro Read | QNH QFE | QFE-> QNH | dQFE'
fmtdat = ' {0:7.2f} | {1:7.2f} {2:7.2f} | {3:7.2f} | {4:5.3f}'
print(fmtalt.format(baro_alt, arp_alt))
print()
print(fmthead)
print(LINE1)
ii = 0
for p in press:
print(fmtdat.format(p, qnh_baros[ii], qfes[ii], qnh_arps[ii], dqfes[ii]))
ii = ii + 1
return None
def _test_plot_isa_():
ZMIN = -500. # gpm
ZMAX = 32000. # gpm
DZ = 20.0 # gpm
try:
import matplotlib.pyplot as plt
except:
print("WARNING: library matplotlib is not installed")
print("This is not vital for the module **admeteopy**")
print()
return -5
hhh = np.arange(ZMIN, ZMAX + DZ, DZ)
plt.subplot(1, 4, 1)
plt.plot(k2c(isa_temperature(hhh)), hhh, 'b-')
plt.grid()
plt.title('ISA-Temperature')
plt.xlabel(" Celsius")
plt.ylabel("gpm")
plt.ylim(ZMIN, ZMAX)
plt.subplot(1, 4, 2)
plt.plot(isa_pressure(hhh), hhh, 'b-')
plt.grid()
plt.title('ISA-Pressure')
plt.ylim(ZMIN, ZMAX)
plt.xlabel("hPa")
dens = isa_density(hhh)
densalt = isa_densityaltitude(dens)
plt.subplot(1, 4, 3)
plt.plot(isa_density(hhh), hhh, 'b-')
plt.grid()
plt.xlabel("$kg/m^3$")
plt.title('ISA-Density')
plt.ylim(ZMIN, ZMAX)
plt.subplot(1, 4, 4)
plt.plot(dens, densalt, 'g-')
plt.grid()
plt.xlabel("$kg/m^3$")
plt.title('ISA-Density Altitude')
plt.ylim(ZMIN, ZMAX)
plt.show()
return 0
if __name__ == "__main__":
_testadisalib_()
_testqnh_()
print(" --- ISA Table ---")
isa_table(csv=True)
_test_plot_isa_()