Merge pull request #511 from RocketPy-Team/enh/precalculate-barometric-height

ENH Precalculate Barometric Height

Author: phmbressan

.pylintrc

@@ -187,6 +187,7 @@ function-naming-style=snake_case
 # Good variable names which should always be accepted, separated by a comma.
 good-names=FlightPhases,
            WindroseAxes,
+           barometric_height_ISA,
 
 
 # Good variable names regexes, separated by a comma. If names match any regex,

CHANGELOG.md

@@ -40,9 +40,11 @@ straightforward as possible.
 - MNT: Add repr method to Parachute class [#490](https://github.com/RocketPy-Team/RocketPy/pull/490)
 - ENH: Function Reverse Arithmetic Priority [#488](https://github.com/RocketPy-Team/RocketPy/pull/488)
 - DOC: Update Header Related Docs
+- ENH Precalculate Barometric Height [#511](https://github.com/RocketPy-Team/RocketPy/pull/511)
 - MNT: Encapsulate quaternion conversions [#537](https://github.com/RocketPy-Team/RocketPy/pull/537)
 - MNT: improve the low pass filter and document an example [#538](https://github.com/RocketPy-Team/RocketPy/pull/538)
 
+
 ### Fixed
 
 - ENH: Parachute trigger doesn't work if "Apogee" is used instead of "apogee" [#489](https://github.com/RocketPy-Team/RocketPy/pull/489)

rocketpy/environment/environment.py

@@ -111,6 +111,10 @@ class Environment:
     Environment.pressure : Function
         Air pressure in Pa as a function of altitude. Can be accessed as regular
         array, or called as a Function. See Function for more information.
+    Environment.barometric_height : Function
+        Geometric height above sea level in m as a function of pressure. Can be
+        accessed as regular array, or called as a Function. See Function for
+        more information.
     Environment.temperature : Function
         Air temperature in K as a function of altitude. Can be accessed as
         regular array, or called as a Function. See Function for more
@@ -1315,6 +1319,8 @@ def process_standard_atmosphere(self):
 
         # Save temperature, pressure and wind profiles
         self.pressure = self.pressure_ISA
+        self.barometric_height = self.barometric_height_ISA
+
         self.temperature = self.temperature_ISA
         self.wind_direction = Function(
             0,
@@ -1433,6 +1439,7 @@ def process_custom_atmosphere(
         if pressure is None:
             # Use standard atmosphere
             self.pressure = self.pressure_ISA
+            self.barometric_height = self.barometric_height_ISA
         else:
             # Use custom input
             self.pressure = Function(
@@ -1441,6 +1448,11 @@ def process_custom_atmosphere(
                 outputs="Pressure (Pa)",
                 interpolation="linear",
             )
+            self.barometric_height = self.pressure.inverse_function().set_discrete(
+                0, max_expected_height, 100, extrapolation="constant"
+            )
+            self.barometric_height.set_inputs("Pressure (Pa)")
+            self.barometric_height.set_outputs("Height Above Sea Level (m)")
             # Check maximum height of custom pressure input
             if not callable(self.pressure.source):
                 max_expected_height = max(self.pressure[-1, 0], max_expected_height)
@@ -1605,6 +1617,15 @@ def process_windy_atmosphere(self, model="ECMWF"):
             outputs="Pressure (Pa)",
             interpolation="linear",
         )
+        # Linearly extrapolate pressure to ground level
+        bar_height = data_array[:, (0, 1)]
+        self.barometric_height = Function(
+            bar_height,
+            inputs="Pressure (Pa)",
+            outputs="Height Above Sea Level (m)",
+            interpolation="linear",
+            extrapolation="natural",
+        )
         self.temperature = Function(
             data_array[:, (1, 2)],
             inputs="Height Above Sea Level (m)",
@@ -1732,6 +1753,15 @@ def process_wyoming_sounding(self, file):
             outputs="Pressure (Pa)",
             interpolation="linear",
         )
+        # Linearly extrapolate pressure to ground level
+        bar_height = data_array[:, (0, 1)]
+        self.barometric_height = Function(
+            bar_height,
+            inputs="Pressure (Pa)",
+            outputs="Height Above Sea Level (m)",
+            interpolation="linear",
+            extrapolation="natural",
+        )
 
         # Retrieve temperature from data array
         data_array[:, 2] = data_array[:, 2] + 273.15  # Converts C to K
@@ -1845,6 +1875,7 @@ def process_noaaruc_sounding(self, file):
 
         # Extract pressure as a function of height
         pressure_array = []
+        barometric_height_array = []
         for line in lines:
             # Split line into columns
             columns = re.split(" +", line)[1:]
@@ -1858,7 +1889,9 @@ def process_noaaruc_sounding(self, file):
                     if max(columns) != 99999:
                         # Save value
                         pressure_array.append(columns)
+                        barometric_height_array.append([columns[1], columns[0]])
         pressure_array = np.array(pressure_array)
+        barometric_height_array = np.array(barometric_height_array)
 
         # Extract temperature as a function of height
         temperature_array = []
@@ -1905,6 +1938,15 @@ def process_noaaruc_sounding(self, file):
             outputs="Pressure (Pa)",
             interpolation="linear",
         )
+        # Converts 10*hPa to Pa and save values
+        barometric_height_array[:, 0] = 10 * barometric_height_array[:, 0]
+        self.barometric_height = Function(
+            barometric_height_array,
+            inputs="Pressure (Pa)",
+            outputs="Height Above Sea Level (m)",
+            interpolation="linear",
+            extrapolation="natural",
+        )
 
         # Convert 10*C to K and save values
         temperature_array[:, 1] = (
@@ -2274,6 +2316,15 @@ def process_forecast_reanalysis(self, file, dictionary):
             outputs="Pressure (Pa)",
             interpolation="linear",
         )
+        # Linearly extrapolate pressure to ground level
+        bar_height = data_array[:, (0, 1)]
+        self.barometric_height = Function(
+            bar_height,
+            inputs="Pressure (Pa)",
+            outputs="Height Above Sea Level (m)",
+            interpolation="linear",
+            extrapolation="natural",
+        )
         self.temperature = Function(
             data_array[:, (1, 2)],
             inputs="Height Above Sea Level (m)",
@@ -2803,6 +2854,15 @@ def select_ensemble_member(self, member=0):
             outputs="Pressure (Pa)",
             interpolation="linear",
         )
+        # Linearly extrapolate pressure to ground level
+        bar_height = data_array[:, (0, 1)]
+        self.barometric_height = Function(
+            bar_height,
+            inputs="Pressure (Pa)",
+            outputs="Height Above Sea Level (m)",
+            interpolation="linear",
+            extrapolation="natural",
+        )
         self.temperature = Function(
             data_array[:, (1, 2)],
             inputs="Height Above Sea Level (m)",
@@ -2965,7 +3025,12 @@ def pressure_function(h):
             outputs="Pressure (Pa)",
         )
 
-        return None
+        # Discretize Function to speed up the trajectory simulation.
+        self.barometric_height_ISA = self.pressure_ISA.inverse_function().set_discrete(
+            pressure[-1], pressure[0], 100, extrapolation="constant"
+        )
+        self.barometric_height_ISA.set_inputs("Pressure (Pa)")
+        self.barometric_height_ISA.set_outputs("Height Above Sea Level (m)")
 
     def calculate_density_profile(self):
         """Compute the density of the atmosphere as a function of

rocketpy/simulation/flight.py

@@ -702,10 +702,7 @@ def __init__(
                     parachute.noisy_pressure_signal.append([node.t, pressure + noise])
                     # Gets height above ground level considering noise
                     hAGL = (
-                        self.env.pressure.find_input(
-                            pressure + noise,
-                            self.y_sol[2],
-                        )
+                        self.env.barometric_height(pressure + noise)
                         - self.env.elevation
                     )
                     if parachute.triggerfunc(pressure + noise, hAGL, self.y_sol):
@@ -1011,10 +1008,7 @@ def __init__(
                                     )
                                     # Gets height above ground level considering noise
                                     hAGL = (
-                                        self.env.pressure.find_input(
-                                            pressure + noise,
-                                            overshootable_node.y[2],
-                                        )
+                                        self.env.barometric_height(pressure + noise)
                                         - self.env.elevation
                                     )
 

tests/test_environment.py

@@ -20,6 +20,7 @@ def test_standard_atmosphere(mock_show, example_env):
     assert example_env.info() == None
     assert example_env.all_info() == None
     assert abs(example_env.pressure(0) - 101325.0) < 1e-8
+    assert abs(example_env.barometric_height(101325.0)) < 1e-2
     assert example_env.prints.print_earth_details() == None
 
 
@@ -43,6 +44,7 @@ def test_custom_atmosphere(mock_show, example_env):
     )
     assert example_env.all_info() == None
     assert abs(example_env.pressure(0) - 101325.0) < 1e-8
+    assert abs(example_env.barometric_height(101325.0)) < 1e-2
     assert abs(example_env.wind_velocity_x(0) - 5) < 1e-8
     assert abs(example_env.temperature(100) - 300) < 1e-8
 
@@ -71,6 +73,7 @@ def test_wyoming_sounding_atmosphere(mock_show, example_env):
             pass
     assert example_env.all_info() == None
     assert abs(example_env.pressure(0) - 93600.0) < 1e-8
+    assert abs(example_env.barometric_height(example_env.pressure(0)) - 722.0) < 1e-8
     assert abs(example_env.wind_velocity_x(0) - -2.9005178894925043) < 1e-8
     assert abs(example_env.temperature(100) - 291.75) < 1e-8