Merge pull request #57 from PaulCreusy/40-add-a-function-to-plot-the-…

…evolution-of-quantities-in-an-atmosphere-classe

40 add a function to plot the evolution of quantities in an atmosphere classe

docs/atmosphere_module.rst

@@ -69,11 +69,19 @@ Here is the list of hypotheses of the model:
 2. The air is dry.
 3. The gravity force is uniform.
 4. The atmosphere is at equilibrium. Therefore :math:`dp = -\rho g_0 dh`
-5. The temperature changes with altitude following these relations:
-
-    1. Between 0 and 11km: :math:`T(h) = 288.15 + h (216.65 - 288.15)/11000`
-    2. Between 11 and 20km: :math:`T(h) = 216.65`
-    3. Between 20 and 32km: :math:`T(h) = 216.65 + (h - 20000) (228.65 - 216.65)/12000`
-    4. Between 32 and 47km: :math:`T(h) = 228.65 + (h - 32000) (270.65 - 228.65)/12000`
+5. The temperature changes with altitude following a table defined in the :download:`US Standard Atmosphere NASA definition from 1976 <_static/US_Standard_Atmosphere_1976.pdf>`.
 
 The formulas of the other quantities directly derive from theses hypotheses.
+
+Graphs
+^^^^^^
+
+Here is an overview of the temperature and pressure evolutions.
+
+.. image:: ./figures/standard_temperature.png
+  :width: 400
+  :alt: Temperature Evolution
+
+.. image:: ./figures/standard_pressure.png
+  :width: 400
+  :alt: Pressure Evolution

flight_mech/atmosphere.py

@@ -6,11 +6,32 @@
 # Imports #
 ###########
 
+# Python imports #
+
+from typing import Literal
+from abc import ABC, abstractmethod
+
 # Dependencies #
 
 import numpy as np
 from scipy.optimize import minimize, Bounds
-from abc import ABC, abstractmethod
+
+# Local imports #
+
+from flight_mech._common import plot_graph
+
+#############
+# Constants #
+#############
+
+VARIABLE_TO_UNIT = {
+    "pressure": "Pa",
+    "temperature": "K",
+    "density": "kg.m-3",
+    "sound_speed": "m.s-1",
+    "dynamic_viscosity": "kg.m-1.s-1",
+    "kinematic_viscosity": "m2.s-1"
+}
 
 #############
 # Functions #
@@ -101,6 +122,49 @@ def compute_altitude_from_sigma(sigma: float) -> float:
         """
         pass
 
+    def plot_graph(
+            self,
+            variable: Literal["pressure", "temperature", "density", "sound_speed", "dynamic_viscosity"],
+            min_altitude: float = 0.,
+            max_altitude: float = 20000.,
+            nb_points: int = 100,
+            **kwargs):
+        """
+        Plot the graph of evolution of the given variable in the atmosphere.
+
+        Parameters
+        ----------
+        variable : Literal["pressure", "temperature", "density", "sound_speed", "dynamic_viscosity"]
+            Name of the variable to plot.
+        min_altitude : float, optional
+            Minimum altitude on the plot in meters, by default 0.
+        max_altitude : float, optional
+            Maximum altitude on the plot in meters, by default 20000.
+        nb_points : int, optional
+            Number of points in the plot.
+
+        Note
+        ----
+        For more details on the optional arguments, please check flight_mech._common.plot_graph.
+        """
+
+        # Define an altitude array
+        altitude_array = np.linspace(min_altitude, max_altitude, nb_points)
+
+        # Compute the variable array
+        variable_array = getattr(
+            self, f"compute_{variable}_from_altitude")(altitude_array)
+
+        # Plot
+        plot_graph(
+            x_array=variable_array,
+            y_array=altitude_array,
+            title=f"{variable.capitalize()} graph",
+            y_label="Altitude in meters",
+            x_label=f"{variable.capitalize()} [{VARIABLE_TO_UNIT[variable]}]",
+            **kwargs
+        )
+
 class ConstantAtmosphere(AtmosphereModel):
     """
     A constant atmosphere model. Used for test purposes only.
@@ -183,7 +247,7 @@ def compute_temperature_from_altitude(z: float):
             (z >= 32000) * (z < 47000) + \
             270.15 * (z >= 47000) * (z < 51000) + \
             (270.15 + (z - 51000) * -2.8 / 1000) * (z >= 51000) * (z < 71000) + \
-            (214.15 + (z - 71000) * -2. / 1000) * (z >= 71000) * (z < 86000)
+            (214.15 + (z - 71000) * -2. / 1000) * (z >= 71000) * (z <= 86000)
 
         return temperature
 
@@ -203,16 +267,16 @@ def compute_pressure_from_altitude(z: float):
             Air pressure in Pa.
         """
 
-        pressure = (101325 * np.power(1 - 22.588e-6 * z, 5.256)) * (z < 11000) + \
-            (22632 * np.exp(-157.77e-6 * (z - 11000))) * (z >= 11000) * (z < 20000) + \
-            (5474.9 * np.power(1 + 4.615e-6 * (z - 20000), 34.163)) * (z >= 20000) * (z < 32000) + \
+        pressure = (z < 11000) * (101325 * np.power(np.abs(1 - 22.588e-6 * z), 5.256)) + \
+            (z >= 11000) * (z < 20000) * (22632 * np.exp(-157.77e-6 * (z - 11000))) + \
+            (5474.9 * np.power(1 - 4.615e-6 * (z - 20000), 34.163)) * (z >= 20000) * (z < 32000) + \
             (868.014 * np.power(1 + 12.245e-6 * (z - 32000), -12.2)) * \
             (z >= 32000) * (z < 47000) + \
             (110.906 * np.exp(-126.293e-6 * (z - 47000))) * (z >= 47000) * (z < 51000) + \
-            (66.939 * np.power((270.65) / (270.65 - 2.8 * (z - 51000) / 1000), 12.2)) * \
+            (66.939 * np.power((270.65) / (270.65 - 2.8 * (z - 51000) / 1000), -12.2)) * \
             (z >= 51000) * (z < 71000) + \
             (3.9564 * np.power((214.65) / (214.65 - 2. * (z - 71000) / 1000), -17.09)) * \
-            (z >= 71000) * (z < 86000)
+            (z >= 71000) * (z <= 86000)
 
         return pressure
 

pyproject.toml

@@ -36,6 +36,9 @@ keywords = [
     "wing",
     "turbojet",
     "airfoil-analysis",
+    "atmosphere",
+    "gas",
+    "motor",
 ]
 
 [project.urls]