diff --git a/.coverage b/.coverage index 87412d1..465e4a3 100644 Binary files a/.coverage and b/.coverage differ diff --git a/CHANGELOG.md b/CHANGELOG.md index 962d6eb..06a6287 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,10 +4,18 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to semantic versioning: [SemVer](https://semver.org/). +## [0.6.1] - 2020-2-17 +### Added +- Unit tests for density heterosphere model + +### Fixed +- Density calculation in `env.py` + ## [0.6.0] - 2019-11-2 ### Added - Unit tests for heterosphere model - Unit tests for `params.py` +- Linting using `flake8` ### Changed - `env.py` for 86km - 1000km support diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index ccc082a..67a6239 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,20 +1,20 @@ # Contributing When contributing to this repository, please first discuss the change you wish to make via issue, -email, or any other method with the owners of this repository before making a change. +email, or any other method with the owners of this repository before making a change. Please note we have a code of conduct, please follow it in all your interactions with the project. ## Pull Request Process -1. Ensure any install or build dependencies are removed before the end of the layer when doing a - build. Also, please properly use the `.gitignore`. +1. Ensure any install or build dependencies are removed before the end of the layer when doing a + build. 2. Commits emojis should follow the [Gitmoji](https://gitmoji.carloscuesta.me/) scheme. -3. Update the README.md and the CHANGELOG.md with details of changes to the interface, this includes new environment +3. Update the README.md and the CHANGELOG.md with details of changes to the interface, this includes new environment variables, exposed ports, useful file locations and container parameters. 4. Increase the version numbers in any examples files and the README.md to the new version that this Pull Request would represent. The versioning scheme we use is [SemVer](http://semver.org/). -5. You may merge the Pull Request in once you have the sign-off of two other developers, or if you +5. You may merge the Pull Request in once you have the sign-off of two other developers, or if you do not have permission to do that, you may request the second reviewer to merge it for you. ## Code of Conduct diff --git a/README.md b/README.md index fdffe28..cfca0ec 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ [![Gitmoji](https://img.shields.io/badge/gitmoji-%20%F0%9F%98%9C%20%F0%9F%98%8D-FFDD67.svg?)](https://gitmoji.carloscuesta.me/) -Python module for **S**imulation and **O**ptimization of **R**ocket **T**rajectories, [**SORT**](https://www.nasa.gov/pdf/140648main_ESAS_17a.pdf) for short. +Python module for **S**imulation and **O**ptimization of **R**ocket **T**rajectories, [**SORT**](https://www.nasa.gov/pdf/140648main_ESAS_17a.pdf#page=19) for short. - Altitude prediction - State prediction - Flight trajectory optimization diff --git a/preflightpy/env.py b/preflightpy/env.py index 440b907..a33d1ab 100644 --- a/preflightpy/env.py +++ b/preflightpy/env.py @@ -27,7 +27,7 @@ class Environment: def __init__(self, vars): # Environmental Constants self.elev, self.t, self.g, self.M_air, self.R, self.gamma, self.Pstatic = vars # noqa - self.hb = [ + self.hb = [ # Layer base altitudes 0, 11000, 20000, @@ -36,7 +36,7 @@ def __init__(self, vars): 51000, 71000 ] - self.Pb = [ + self.Pb = [ # Layer base pressures 101325, 22632.1, 5474.89, @@ -45,7 +45,7 @@ def __init__(self, vars): 66.9389, 3.95642 ] - self.Tb = [ + self.Tb = [ # Layer base temperatures 288.15, 216.65, 216.65, @@ -54,7 +54,7 @@ def __init__(self, vars): 270.65, 214.65 ] - self.Lm = [ + self.Lm = [ # Layer lapse rates -0.0065, 0.0, 0.001, @@ -67,6 +67,10 @@ def __init__(self, vars): def get_geopotential_altitude(self, r: float, z: float) -> float: return r*z / (r+z) + def atmo_heterosphere_equ(self, z: float, a, b, c, d, e): + z_km = z/1000 + return math.exp(a * z_km**4 + b * z_km**3 + c * z_km**2 + d * z_km + e) # noqa + def get_temp(self, z: float, h: float) -> float: if 0 <= h <= 11000: return (288.15 + (self.Lm[0]*(h-0)), 0) @@ -80,10 +84,10 @@ def get_temp(self, z: float, h: float) -> float: return (270.65 + (self.Lm[4]*(h-47000)), 4) elif 51000 < h <= 71000: return (270.65 + (self.Lm[5]*(h-51000)), 5) - elif 71000 < h <= 84856: + elif 71000 < h <= 84852: return (214.65 + (self.Lm[6]*(h-71000)), 6) elif 86000 < z <= 91000: - return (186.67, 7) + return (186.87, 7) elif 91000 < z <= 110000: if 91000 < z <= 100000: layer = 8 @@ -113,17 +117,14 @@ def get_temp(self, z: float, h: float) -> float: def get_pressure(self, z: float, h: float, T: float, b: int) -> float: - def equ(a, b, c, d, e): - zeta = z/1000 - return math.exp(a * zeta**4 + b * zeta**3 + c * zeta**2 + d * zeta + e) # noqa - if b <= 6: if self.Lm[b] != 0: return self.Pb[b] * (self.Tb[b]/T)**(self.g*self.M_air/(self.R*self.Lm[b])) # noqa else: return self.Pb[b] * math.exp(-self.g * self.M_air * (h-self.hb[b]) / (self.R*self.Tb[b])) # noqa elif b == 7: - return equ( + return self.atmo_heterosphere_equ( + z, 0.000000, 2.159582e-6, -4.836957e-4, @@ -131,7 +132,8 @@ def equ(a, b, c, d, e): 13.47530 ) elif b == 8: - return equ( + return self.atmo_heterosphere_equ( + z, 0.000000, 3.304895e-5, -0.009062730, @@ -139,7 +141,8 @@ def equ(a, b, c, d, e): -11.03037 ) elif b == 9: - return equ( + return self.atmo_heterosphere_equ( + z, 0.000000, 6.693926e-5, -0.01945388, @@ -147,7 +150,8 @@ def equ(a, b, c, d, e): -47.75030 ) elif b == 10: - return equ( + return self.atmo_heterosphere_equ( + z, 0.000000, -6.539316e-5, 0.02485568, @@ -155,7 +159,8 @@ def equ(a, b, c, d, e): 135.9355 ) elif b == 11: - return equ( + return self.atmo_heterosphere_equ( + z, 2.283506e-7, -1.343221e-4, 0.02999016, @@ -163,7 +168,8 @@ def equ(a, b, c, d, e): 113.5764 ) elif b == 12: - return equ( + return self.atmo_heterosphere_equ( + z, 1.209434e-8, -9.692458e-6, 0.003002041, @@ -171,7 +177,8 @@ def equ(a, b, c, d, e): 19.19151 ) elif b == 13: - return equ( + return self.atmo_heterosphere_equ( + z, 8.113942e-10, -9.822568e-7, 4.687616e-4, @@ -179,7 +186,8 @@ def equ(a, b, c, d, e): 3.067409 ) elif b == 14: - return equ( + return self.atmo_heterosphere_equ( + z, 9.814674e-11, -1.654439e-7, 1.148115e-4, @@ -187,7 +195,8 @@ def equ(a, b, c, d, e): -2.011365 ) elif b == 15: - return equ( + return self.atmo_heterosphere_equ( + z, -7.835161e-11, 1.964589e-7, -1.657213e-4, @@ -195,7 +204,8 @@ def equ(a, b, c, d, e): -14.77132 ) elif b == 16: - return equ( + return self.atmo_heterosphere_equ( + z, 2.813255e-11, -1.120689e-7, 1.695568e-4, @@ -203,19 +213,106 @@ def equ(a, b, c, d, e): 14.56718 ) - def get_density(self, P: float, T: float, b) -> float: + def get_density(self, z: float, P: float, T: float, b) -> float: if b <= 6: - V = 1 - n = (P*V)/(self.R*T) - m = self.M_air * n - return (P * m)/(n * self.R * T) + return (P * self.M_air)/(self.R * T) + elif b == 7: + return self.atmo_heterosphere_equ( + z, + 0.000000, + -3.322622E-06, + 9.111460E-04, + -0.2609971, + 5.944694 + ) + elif b == 8: + return self.atmo_heterosphere_equ( + z, + 0.000000, + 2.873405e-05, + -0.008492037, + 0.6541179, + -23.62010 + ) + elif b == 9: + return self.atmo_heterosphere_equ( + z, + -1.240774e-05, + 0.005162063, + -0.8048342, + 55.55996, + -1443.338 + ) + elif b == 10: + return self.atmo_heterosphere_equ( + z, + 0.00000, + -8.854164e-05, + 0.03373254, + -4.390837, + 176.5294 + ) + elif b == 11: + return self.atmo_heterosphere_equ( + z, + 3.661771e-07, + -2.154344e-04, + 0.04809214, + -4.884744, + 172.3597 + ) + elif b == 12: + return self.atmo_heterosphere_equ( + z, + 1.906032e-08, + -1.527799E-05, + 0.004724294, + -0.6992340, + 20.50921 + ) + elif b == 13: + return self.atmo_heterosphere_equ( + z, + 1.199282e-09, + -1.451051e-06, + 6.910474e-04, + -0.1736220, + -5.321644 + ) + elif b == 14: + return self.atmo_heterosphere_equ( + z, + 1.140564e-10, + -2.130756e-07, + 1.570762e-04, + -0.07029296, + -12.89844 + ) + elif b == 15: + return self.atmo_heterosphere_equ( + z, + 8.105631e-12, + -2.358417e-09, + -2.635110e-06, + -0.01562608, + -20.02246 + ) + elif b == 16: + return self.atmo_heterosphere_equ( + z, + -3.701195e-12, + -8.608611e-09, + 5.118829e-05, + -0.06600998, + -6.137674 + ) - def get_c(self, T: float): + def get_c(self, T: float) -> float: return math.sqrt((self.gamma * self.R * T) / self.M_air) def get_status(self, z: float): h = round(self.get_geopotential_altitude(6356766, z), 0) self.T, b = self.get_temp(z, h) self.P = self.get_pressure(z, h, self.T, b) - self.Rho = self.get_density(self.P, self.T, b) + self.Rho = self.get_density(z, self.P, self.T, b) self.c = self.get_c(self.T) diff --git a/setup.py b/setup.py index c0e04f1..b01deef 100644 --- a/setup.py +++ b/setup.py @@ -42,7 +42,7 @@ # For a discussion on single-sourcing the version across setup.py and the # project code, see # https://packaging.python.org/en/latest/single_source_version.html - version='0.5.0', # Required + version='0.6.1', # Required # This is a one-line description or tagline of what your project does. This # corresponds to the "Summary" metadata field: diff --git a/tests/env_test.py b/tests/env_test.py index 8a8db25..f48d4a7 100644 --- a/tests/env_test.py +++ b/tests/env_test.py @@ -13,36 +13,39 @@ def test_pressure(self): 1.4, 101325 ]) + # Threshold set to ± 0.1% + environment.get_status(0) + assert abs(100 - round(environment.P, 0) / 101325 * 100) <= 0.1 environment.get_status(5000) - assert round(environment.P, -1) == 54040 + assert abs(100 - round(environment.P, -1) / 54040 * 100) <= 0.1 environment.get_status(11019) - assert round(environment.P, -1) == 22630 + assert abs(100 - round(environment.P, -1) / 22630 * 100) <= 0.1 environment.get_status(20000) - assert round(environment.P, -1) == 5530 + assert abs(100 - round(environment.P, -1) / 5530 * 100) <= 0.1 environment.get_status(70000) - assert round(environment.P, 1) == 5.2 + assert abs(100 - round(environment.P, 1) / 5.2 * 100) <= 0.1 environment.get_status(87000) - assert round(environment.P, 2) == 0.31 + assert abs(100 - round(environment.P, 2) / 0.31 * 100) <= 0.1 environment.get_status(95000) - assert round(environment.P, 3) == 0.076 + assert abs(100 - round(environment.P, 3) / 0.076 * 100) <= 0.1 environment.get_status(105000) - assert round(environment.P, 3) == 0.014 + assert abs(100 - round(environment.P, 3) / 0.014 * 100) <= 0.1 environment.get_status(115000) - assert round(environment.P, 4) == 0.0040 + assert abs(100 - round(environment.P, 4) / 0.0040 * 100) <= 0.1 environment.get_status(140000) - assert round(environment.P, 5) == 0.00072 + assert abs(100 - round(environment.P, 5) / 0.00072 * 100) <= 0.1 environment.get_status(170000) - assert round(environment.P, 5) == 0.00021 + assert abs(100 - round(environment.P, 5) / 0.00021 * 100) <= 0.1 environment.get_status(250000) - assert round(environment.P, 6) == 0.000025 + assert abs(100 - round(environment.P, 6) / 0.000025 * 100) <= 0.1 environment.get_status(400000) - assert round(environment.P, 8) == 1.45e-6 + assert abs(100 - round(environment.P, 8) / 1.45e-6 * 100) <= 0.1 environment.get_status(625000) - assert round(environment.P, 10) == 6.26e-8 + assert abs(100 - round(environment.P, 10) / 6.26e-8 * 100) <= 0.1 environment.get_status(800000) - assert round(environment.P, 10) == 1.70e-8 + assert abs(100 - round(environment.P, 10) / 1.70e-8 * 100) <= 0.1 environment.get_status(1000000) - assert round(environment.P, 11) == 7.51e-9 + assert abs(100 - round(environment.P, 11) / 7.51e-9 * 100) <= 0.1 def test_temperature(self): environment = pre.Environment([ @@ -54,34 +57,35 @@ def test_temperature(self): 1.4, 101325 ]) + # Threshold set to ± 0.1% environment.get_status(5000) - assert round(environment.T, 2) == 255.68 + assert abs(100 - round(environment.T, 2) / 255.68 * 100) <= 0.1 environment.get_status(20063) - assert round(environment.T, 2) == 216.65 + assert abs(100 - round(environment.T, 2) / 216.65 * 100) <= 0.1 environment.get_status(30000) - assert round(environment.T, 2) == 226.51 + assert abs(100 - round(environment.T, 2) / 226.51 * 100) <= 0.1 environment.get_status(47350) - assert round(environment.T, 2) == 270.65 + assert abs(100 - round(environment.T, 2) / 270.65 * 100) <= 0.1 environment.get_status(71802) - assert round(environment.T, 2) == 214.65 + assert abs(100 - round(environment.T, 2) / 214.65 * 100) <= 0.1 environment.get_status(87000) - assert round(environment.T, 2) == 186.67 + assert abs(100 - round(environment.T, 2) / 186.87 * 100) <= 0.1 environment.get_status(100000) - assert round(environment.T, 2) == 195.08 + assert abs(100 - round(environment.T, 2) / 195.08 * 100) <= 0.1 environment.get_status(115000) - assert round(environment.T, 2) == 300.00 + assert abs(100 - round(environment.T, 2) / 300.00 * 100) <= 0.1 environment.get_status(125000) - assert round(environment.T, 2) == 417.23 + assert abs(100 - round(environment.T, 2) / 417.23 * 100) <= 0.1 environment.get_status(160000) - assert round(environment.T, 2) == 696.29 + assert abs(100 - round(environment.T, 2) / 696.29 * 100) <= 0.1 environment.get_status(330000) - assert round(environment.T, 2) == 985.88 + assert abs(100 - round(environment.T, 2) / 985.88 * 100) <= 0.1 environment.get_status(500000) - assert round(environment.T, 2) == 999.24 + assert abs(100 - round(environment.T, 2) / 999.24 * 100) <= 0.1 environment.get_status(750000) - assert round(environment.T, 2) == 999.99 + assert abs(100 - round(environment.T, 2) / 999.99 * 100) <= 0.1 environment.get_status(1000000) - assert round(environment.T, 2) == 1000.00 + assert abs(100 - round(environment.T, 2) / 1000.00 * 100) <= 0.1 def test_density(self): environment = pre.Environment([ @@ -93,14 +97,39 @@ def test_density(self): 1.4, 101325 ]) + # Threshold set to ± 0.1% + environment.get_status(0) + assert abs(100 - round(environment.Rho, 4) / 1.2252 * 100) <= 0.1 environment.get_status(5000) - assert round(environment.Rho, 4) == 0.7365 - environment.get_status(11019) - assert round(environment.Rho, 4) == 0.3639 - environment.get_status(20063) - assert round(environment.Rho, 4) == 0.0880 - environment.get_status(40000.0) - assert round(environment.Rho, 4) == 0.0040 + assert abs(100 - round(environment.Rho, 5) / 7.3643e-1 * 100) <= 0.1 + environment.get_status(11000) + assert abs(100 - round(environment.Rho, 5) / 3.6480e-1 * 100) <= 0.1 + environment.get_status(20000) + assert abs(100 - round(environment.Rho, 6) / 8.8910e-2 * 100) <= 0.1 + environment.get_status(70000) + assert abs(100 - round(environment.Rho, 9) / 8.2829e-5 * 100) <= 0.1 + environment.get_status(87000) + assert abs(100 - round(environment.Rho, 9) / 5.824e-6 * 100) <= 0.1 + environment.get_status(95000) + assert abs(100 - round(environment.Rho, 9) / 1.393e-6 * 100) <= 0.1 + environment.get_status(105000) + assert abs(100 - round(environment.Rho, 10) / 2.325e-7 * 100) <= 0.1 + environment.get_status(115000) + assert abs(100 - round(environment.Rho, 11) / 4.289e-8 * 100) <= 0.1 + environment.get_status(140000) + assert abs(100 - round(environment.Rho, 12) / 3.831e-9 * 100) <= 0.1 + environment.get_status(170000) + assert abs(100 - round(environment.Rho, 13) / 7.815e-10 * 100) <= 0.1 + environment.get_status(250000) + assert abs(100 - round(environment.Rho, 14) / 6.073e-11 * 100) <= 0.1 + environment.get_status(400000) + assert abs(100 - round(environment.Rho, 15) / 2.803e-12 * 100) <= 0.1 + environment.get_status(625000) + assert abs(100 - round(environment.Rho, 17) / 7.998e-14 * 100) <= 0.1 + environment.get_status(800000) + assert abs(100 - round(environment.Rho, 17) / 1.136e-14 * 100) <= 0.1 + environment.get_status(1000000) + assert abs(100 - round(environment.Rho, 18) / 3.561e-15 * 100) <= 0.1 def test_speed_of_sound(self): environment = pre.Environment([ @@ -112,9 +141,10 @@ def test_speed_of_sound(self): 1.4, 101325 ]) + # Threshold set to ± 0.1% environment.get_status(1000) - assert round(environment.c, 1) == 336.4 + assert abs(100 - round(environment.c, 1) / 336.4 * 100) <= 0.1 environment.get_status(10000) - assert round(environment.c, 1) == 299.5 + assert abs(100 - round(environment.c, 1) / 299.5 * 100) <= 0.1 environment.get_status(30000) - assert round(environment.c, 1) == 301.7 + assert abs(100 - round(environment.c, 1) / 301.7 * 100) <= 0.1