Skip to content
This repository has been archived by the owner on Aug 13, 2020. It is now read-only.

Commit

Permalink
Merge pull request #17 from AlexS12/isa
Browse files Browse the repository at this point in the history 
Add 1984 ISA 1976 Standard atmosphere with tests
  • Loading branch information
@astrojuanlu
astrojuanlu committed Jan 4, 2016
2 parents f323ab4 + 093f6ff commit c2bed14
Show file tree
Hide file tree
Showing 3 changed files with 345 additions and 0 deletions.
0 ...fme/environment/tests/test_environment.py → src/pyfme/environment/__init__.py
View file
File renamed without changes.
142 changes: 142 additions & 0 deletions src/pyfme/environment/isa.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,142 @@
# coding: utf-8
"""ISA functions. Implementation based on:
.. [1] U.S. Standard Atmosphere, 1976, U.S. Government Printing Office,
Washington, D.C., 1976
From: https://en.wikipedia.org/wiki/U.S._Standard_Atmosphere
Geopotential
altitude Static Standard Temperature
Subscript above MSL Pressure Temperature Lapse Rate
-------------------------------------------------------------
(m) (pascals) (K) (K/m)
0 0 101325 288.15 -0.0065
1 11000 22632.1 216.65 0
2 20000 5474.89 216.65 0.001
3 32000 868.019 228.65 0.0028
4 47000 110.906 270.65 0
5 51000 66.9389 270.65 -0.0028
6 71000 3.95642 214.65 -0.002
"""

from math import exp

# Constants
R_a = 287.05287 # J/(Kg·K)
g0 = 9.80665 # m/s^2

# Layer constants
h0 = (0, 11000, 20000, 32000, 47000, 51000, 71000, 84500) # m
T0_layers = (288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65) # K
p0_layers = (101325.0, 22632.1, 5474.89, 868.019, 110.906, 66.9389, 3.95642) # Pa
alpha_layers = (-0.0065, 0, 0.001, 0.0028, 0, -0.0028, -0.002) # K / m


def atm(h):
"""ISA 1976 Standard atmosphere temperature, pressure and density.
Parameters
----------
h : float
Geopotential altitude (m). h values must range from 0 to 84500 m.
Returns
-------
T : float
Temperature (K).
p : float
Pressure (Pa).
rho : float
Density (kg/m³)
Raises
------
ValueError
If the value of the altitude is outside the defined layers.
Notes
-----
Check layers and reference values in [2].
References
----------
.. [1] U.S. Standard Atmosphere, 1976, U.S. Government Printing Office,
Washington, D.C., 1976
.. [2] https://en.wikipedia.org/wiki/U.S._Standard_Atmosphere
"""

if h < 0.0:
raise ValueError("Altitude cannot be less than 0 m.")

elif h0[0] <= h < h0[1]: # Troposphere
T0 = T0_layers[0]
p0 = p0_layers[0]
alpha = alpha_layers[0]

T = T0 + alpha * h
p = p0 * (T0 / (T0 + alpha * h)) ** (g0 / (R_a * alpha))

elif h0[1] <= h < h0[2]: # Tropopause
T0 = T0_layers[1]
p0 = p0_layers[1]
alpha = alpha_layers[1]
h_diff = h - h0[1]

T = T0
p = p0 * exp(-g0 * h_diff / (R_a * T0))

elif h0[2] <= h < h0[3]: # Stratosphere 1
T0 = T0_layers[2]
p0 = p0_layers[2]
alpha = alpha_layers[2]
h_diff = h - h0[2]

T = T0 + alpha * h_diff
p = p0 * (T0 / (T0 + alpha * h_diff)) ** (g0 / (R_a * alpha))

elif h0[3] <= h < h0[4]: # Stratosphere 2
T0 = T0_layers[3]
p0 = p0_layers[3]
alpha = alpha_layers[3]
h_diff = h - h0[3]

T = T0 + alpha * h_diff
p = p0 * (T0 / (T0 + alpha * h_diff)) ** (g0 / (R_a * alpha))

elif h0[4] <= h < h0[5]: # Stratopause
T0 = T0_layers[4]
p0 = p0_layers[4]
alpha = alpha_layers[4]
h_diff = h - h0[4]

T = T0
p = p0 * exp(-g0 * h_diff / (R_a * T0))

elif h0[5] <= h < h0[6]: # Mesosphere 1
T0 = T0_layers[5]
p0 = p0_layers[5]
alpha = alpha_layers[5]
h_diff = h - h0[5]

T = T0 + alpha * h_diff
p = p0 * (T0 / (T0 + alpha * h_diff)) ** (g0 / (R_a * alpha))

elif h0[6] <= h <= h0[7]: # Mesosphere 2
T0 = T0_layers[6]
p0 = p0_layers[6]
alpha = alpha_layers[6]
h_diff = h - h0[6]

T = T0 + alpha * h_diff
p = p0 * (T0 / (T0 + alpha * h_diff)) ** (g0 / (R_a * alpha))

else:
raise ValueError("Altitude cannot be greater than {} m.".format(h0[7]))

rho = p / (R_a * T)

return T, p, rho
203 changes: 203 additions & 0 deletions src/pyfme/environment/tests/test_isa.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,203 @@
# coding: utf-8
"""Tests of the ISA functions.
All numerical results are validated against the `COESA`_ standard.
.. _`COESA`: http://hdl.handle.net/2060/19770009539
Based on scikit-aero (c) 2012 scikit-aero authors.
"""
import numpy as np
from numpy.testing import (assert_equal, assert_almost_equal)

import pytest

from pyfme.environment.isa import atm


def test_sea_level():
h = 0.0 # m
expected_T = 288.15 # K
expected_p = 101325.0 # Pa
expected_rho = 1.2250 # kg / m3

T, p, rho = atm(h)

# Reads: "Assert if T equals expected_T"
assert_equal(T, expected_T)
assert_equal(p, expected_p)
assert_almost_equal(rho, expected_rho, decimal=4)


def test_altitude_is_out_of_range(recwarn):
wrong_h = (-1.0, 84501) # m

for h in wrong_h:
with pytest.raises(ValueError):
atm(h)


def test_results_under_11km():
h = np.array([0.0,
50.0,
550.0,
6500.0,
10000.0,
11000.0
]) # m
expected_T = np.array([288.150,
287.825,
284.575,
245.900,
223.150,
216.650
]) # K

expected_rho = np.array([1.2250,
1.2191,
1.1616,
0.62384,
0.41271,
0.36392
]) # kg / m3

for ii, h_ in enumerate(h):
T, p, rho = atm(h_)
assert_almost_equal(T, expected_T[ii], decimal=3)
assert_almost_equal(rho, expected_rho[ii], decimal=4)


def test_results_under_20km():
h = np.array([12000,
14200,
17500,
20000
]) # m
expected_T = np.array([216.650,
216.650,
216.650,
216.650,
]) # K

expected_rho = np.array([0.31083,
0.21971,
0.13058,
0.088035
]) # kg / m3

for ii, h_ in enumerate(h):
T, p, rho = atm(h_)
assert_almost_equal(T, expected_T[ii], decimal=3)
assert_almost_equal(rho, expected_rho[ii], decimal=4)


def test_results_under_32km():
h = np.array([22100,
24000,
28800,
32000
]) # m
expected_T = np.array([218.750,
220.650,
225.450,
228.650
]) # K

expected_rho = np.array([0.062711,
0.046267,
0.021708,
0.013225
]) # kg / m3

for ii, h_ in enumerate(h):
T, p, rho = atm(h_)
assert_almost_equal(T, expected_T[ii], decimal=3)
assert_almost_equal(rho, expected_rho[ii], decimal=4)


def test_results_under_47km():
h = np.array([32200,
36000,
42000,
47000
]) # m
expected_T = np.array([229.210,
239.850,
256.650,
270.650
]) # K

expected_rho = np.array([0.012805,
0.0070344,
0.0028780,
0.0014275
]) # kg / m3

for ii, h_ in enumerate(h):
T, p, rho = atm(h_)
assert_almost_equal(T, expected_T[ii], decimal=3)
assert_almost_equal(rho, expected_rho[ii], decimal=4)


def test_results_under_51km():
h = np.array([47200,
49000,
51000
]) # m
expected_T = np.array([270.650,
270.650,
270.650
]) # K

expected_rho = np.array([0.0013919,
0.0011090,
0.00086160
]) # kg / m3

for ii, h_ in enumerate(h):
T, p, rho = atm(h_)
assert_almost_equal(T, expected_T[ii], decimal=3)
assert_almost_equal(rho, expected_rho[ii], decimal=4)


def test_results_under_71km():
h = np.array([51500,
60000,
71000
]) # m
expected_T = np.array([269.250,
245.450,
214.650
]) # K

expected_rho = np.array([0.00081298,
2.8832e-4,
6.4211e-5
]) # kg / m3

for ii, h_ in enumerate(h):
T, p, rho = atm(h_)
assert_almost_equal(T, expected_T[ii], decimal=3)
assert_almost_equal(rho, expected_rho[ii], decimal=4)


def test_results_under_84km():
h = np.array([52000,
60000,
84500
]) # m
expected_T = np.array([267.850,
245.450,
187.650
]) # K

expected_rho = np.array([7.6687e-4,
2.8832e-4,
7.3914e-6
]) # kg / m3

for ii, h_ in enumerate(h):
T, p, rho = atm(h_)
assert_almost_equal(T, expected_T[ii], decimal=3)
assert_almost_equal(rho, expected_rho[ii], decimal=4)

0 comments on commit c2bed14

Please sign in to comment.