From c38830d1ac474edad42e0f74ac17104486a195fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Fri, 1 Jan 2016 19:06:10 +0100 Subject: [PATCH 1/3] Added 1984 ISA 1976 Standard atmosphere with tests. These version of ISA is designed to be called with only one value of altitude instead of an array. Performance improvements can be done with numba. --- .../test_environment.py => __init__.py} | 0 src/pyfme/environment/isa.py | 132 ++++++++++++ src/pyfme/environment/tests/test_isa.py | 203 ++++++++++++++++++ 3 files changed, 335 insertions(+) rename src/pyfme/environment/{tests/test_environment.py => __init__.py} (100%) create mode 100644 src/pyfme/environment/isa.py create mode 100644 src/pyfme/environment/tests/test_isa.py diff --git a/src/pyfme/environment/tests/test_environment.py b/src/pyfme/environment/__init__.py similarity index 100% rename from src/pyfme/environment/tests/test_environment.py rename to src/pyfme/environment/__init__.py diff --git a/src/pyfme/environment/isa.py b/src/pyfme/environment/isa.py new file mode 100644 index 0000000..2f622cf --- /dev/null +++ b/src/pyfme/environment/isa.py @@ -0,0 +1,132 @@ +# coding: utf-8 +"""ISA functions. +""" + +from math import exp + +# Constants +R_a = 287.05287 # J/(Kg·K) +g0 = 9.80665 # m/s^2 + +# 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 + +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). + + 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. + + 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 diff --git a/src/pyfme/environment/tests/test_isa.py b/src/pyfme/environment/tests/test_isa.py new file mode 100644 index 0000000..5617a8c --- /dev/null +++ b/src/pyfme/environment/tests/test_isa.py @@ -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) \ No newline at end of file From ad057d3e43b346ba02484e6351084bf0c3135aba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Fri, 1 Jan 2016 20:14:03 +0100 Subject: [PATCH 2/3] Style corrections --- src/pyfme/environment/isa.py | 7 ++++++- src/pyfme/environment/tests/test_isa.py | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/pyfme/environment/isa.py b/src/pyfme/environment/isa.py index 2f622cf..02c781e 100644 --- a/src/pyfme/environment/isa.py +++ b/src/pyfme/environment/isa.py @@ -1,5 +1,9 @@ # coding: utf-8 -"""ISA functions. +"""ISA functions. Implementation based on: + +.. [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 """ from math import exp @@ -29,6 +33,7 @@ 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. diff --git a/src/pyfme/environment/tests/test_isa.py b/src/pyfme/environment/tests/test_isa.py index 5617a8c..532cd5c 100644 --- a/src/pyfme/environment/tests/test_isa.py +++ b/src/pyfme/environment/tests/test_isa.py @@ -200,4 +200,4 @@ def test_results_under_84km(): 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) \ No newline at end of file + assert_almost_equal(rho, expected_rho[ii], decimal=4) From 093f6ffb31213bb0b46350c8bbb9d29ec8ff9d08 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=81lex=20S=C3=A1ez?= Date: Fri, 1 Jan 2016 20:28:09 +0100 Subject: [PATCH 3/3] Moved table with reference values to module docstring and added a notes section and the available altitude range in the function. --- src/pyfme/environment/isa.py | 41 ++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/src/pyfme/environment/isa.py b/src/pyfme/environment/isa.py index 02c781e..64eb2cb 100644 --- a/src/pyfme/environment/isa.py +++ b/src/pyfme/environment/isa.py @@ -3,7 +3,23 @@ .. [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 + + 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 @@ -12,22 +28,7 @@ R_a = 287.05287 # J/(Kg·K) g0 = 9.80665 # m/s^2 -# 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 - +# 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 @@ -40,7 +41,7 @@ def atm(h): Parameters ---------- h : float - Geopotential altitude (m). + Geopotential altitude (m). h values must range from 0 to 84500 m. Returns ------- @@ -56,6 +57,10 @@ def atm(h): 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,