Merge pull request #17 from Oxyde2/master

v0.6.1
RimorRes · Feb 17, 2020 · 823db5f · 823db5f
2 parents 662e086 + 8542845
commit 823db5f
Show file tree

Hide file tree

Showing 7 changed files with 208 additions and 73 deletions.
diff --git a/.coverage b/.coverage
diff --git 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
@@ -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
@@ -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
@@ -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,109 +117,202 @@ 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,
                 -0.1425192,
                 13.47530
                 )
         elif b == 8:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 0.000000,
                 3.304895e-5,
                 -0.009062730,
                 0.6516698,
                 -11.03037
                 )
         elif b == 9:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 0.000000,
                 6.693926e-5,
                 -0.01945388,
                 1.719080,
                 -47.75030
                 )
         elif b == 10:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 0.000000,
                 -6.539316e-5,
                 0.02485568,
                 -3.223620,
                 135.9355
                 )
         elif b == 11:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 2.283506e-7,
                 -1.343221e-4,
                 0.02999016,
                 -3.055446,
                 113.5764
                 )
         elif b == 12:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 1.209434e-8,
                 -9.692458e-6,
                 0.003002041,
                 -0.4523015,
                 19.19151
                 )
         elif b == 13:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 8.113942e-10,
                 -9.822568e-7,
                 4.687616e-4,
                 -0.1231710,
                 3.067409
                 )
         elif b == 14:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 9.814674e-11,
                 -1.654439e-7,
                 1.148115e-4,
                 -0.05431334,
                 -2.011365
                 )
         elif b == 15:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 -7.835161e-11,
                 1.964589e-7,
                 -1.657213e-4,
                 0.04305869,
                 -14.77132
                 )
         elif b == 16:
-            return equ(
+            return self.atmo_heterosphere_equ(
+                z,
                 2.813255e-11,
                 -1.120689e-7,
                 1.695568e-4,
                 -0.1188941,
                 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
@@ -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: