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..64eb2cb --- /dev/null +++ b/src/pyfme/environment/isa.py @@ -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 diff --git a/src/pyfme/environment/tests/test_isa.py b/src/pyfme/environment/tests/test_isa.py new file mode 100644 index 0000000..532cd5c --- /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)