ENH: Improve parachute geometric parametrization

README.md

@@ -261,6 +261,9 @@ main = calisto.add_parachute(
     sampling_rate=105,
     lag=1.5,
     noise=(0, 8.3, 0.5),
+    parachute_radius=1.5,
+    parachute_height=1.5,
+    porosity=0.0432,
 )
 
 drogue = calisto.add_parachute(
@@ -270,6 +273,9 @@ drogue = calisto.add_parachute(
     sampling_rate=105,
     lag=1.5,
     noise=(0, 8.3, 0.5),
+    parachute_radius=1.5,
+    parachute_height=1.5,
+    porosity=0.0432,
 )
 ```
 

docs/user/first_simulation.rst

@@ -276,6 +276,9 @@ Finally, we can add any number of Parachutes to the ``Rocket`` object.
         sampling_rate=105,
         lag=1.5,
         noise=(0, 8.3, 0.5),
+        parachute_radius=1.5,
+        parachute_height=1.5,
+        porosity=0.0432,
     )
 
     drogue = calisto.add_parachute(
@@ -285,6 +288,9 @@ Finally, we can add any number of Parachutes to the ``Rocket`` object.
         sampling_rate=105,
         lag=1.5,
         noise=(0, 8.3, 0.5),
+        parachute_radius=1.5,
+        parachute_height=1.5,
+        porosity=0.0432,
     )
 
 We can then see if the rocket is stable by plotting the static margin:

docs/user/rocket/rocket_usage.rst

@@ -302,6 +302,9 @@ apogee and another that will be deployed at 800 meters above ground level:
         sampling_rate=105,
         lag=1.5,
         noise=(0, 8.3, 0.5),
+        parachute_radius=1.5,
+        parachute_height=1.5,
+        porosity=0.0432,
     )
 
     drogue = calisto.add_parachute(
@@ -311,6 +314,9 @@ apogee and another that will be deployed at 800 meters above ground level:
         sampling_rate=105,
         lag=1.5,
         noise=(0, 8.3, 0.5),
+        parachute_radius=1.5,
+        parachute_height=1.5,
+        porosity=0.0432,
     )
 
 .. seealso::

rocketpy/rocket/parachute.py

@@ -9,7 +9,7 @@
 
 
 class Parachute:
-    """Keeps parachute information.
+    """Keeps information of the parachute, which is modeled as a hemispheroid.
 
     Attributes
     ----------
@@ -92,6 +92,15 @@ class Parachute:
         Function of noisy_pressure_signal.
     Parachute.clean_pressure_signal_function : Function
         Function of clean_pressure_signal.
+    Parachute.parachute_radius : float
+        Length of the non-unique semi-axis (radius) of the inflated hemispheroid
+        parachute in meters.
+    Parachute.parachute_height : float
+        Length of the unique semi-axis (height) of the inflated hemispheroid
+        parachute in meters.
+    Parachute.porosity : float
+        Porosity of the parachute is the ratio of open space in the canopy to total
+        canopy area.
     """
 
     def __init__(
@@ -102,6 +111,9 @@ def __init__(
         sampling_rate,
         lag=0,
         noise=(0, 0, 0),
+        parachute_radius=1.5,
+        parachute_height=None,
+        porosity=0.0432,
     ):
         """Initializes Parachute class.
 
@@ -154,6 +166,17 @@ def __init__(
             The values are used to add noise to the pressure signal which is
             passed to the trigger function. Default value is ``(0, 0, 0)``.
             Units are in Pa.
+        parachute_radius : float, optional
+            Length of the non-unique semi-axis (radius) of the inflated hemispheroid
+            parachute. Default value is 1.5.
+            Units are in meters.
+        parachute_height : float, optional
+            Length of the unique semi-axis (height) of the inflated hemispheroid
+            parachute. Default value is the radius of the parachute.
+            Units are in meters.
+        porosity : float, optional
+            Porosity of the parachute is the ratio of open space in the canopy to total
+            canopy area. Default value is 0.0432 (for consistency with previous versions).
         """
         self.name = name
         self.cd_s = cd_s
@@ -170,6 +193,11 @@ def __init__(
         self.clean_pressure_signal_function = Function(0)
         self.noisy_pressure_signal_function = Function(0)
         self.noise_signal_function = Function(0)
+        self.parachute_radius = parachute_radius
+        if parachute_height is None:
+            parachute_height = parachute_radius
+        self.parachute_height = parachute_height
+        self.porosity = porosity
 
         alpha, beta = self.noise_corr
         self.noise_function = lambda: alpha * self.noise_signal[-1][
@@ -264,6 +292,9 @@ def to_dict(self, include_outputs=False):
             "sampling_rate": self.sampling_rate,
             "lag": self.lag,
             "noise": self.noise,
+            "parachute_radius": self.parachute_radius,
+            "parachute_height": self.parachute_height,
+            "porosity": self.porosity,
         }
 
         if include_outputs:
@@ -290,6 +321,9 @@ def from_dict(cls, data):
             sampling_rate=data["sampling_rate"],
             lag=data["lag"],
             noise=data["noise"],
+            parachute_radius=data.get("parachute_radius", 1.5),
+            parachute_height=data.get("parachute_height", None),
+            porosity=data.get("porosity", 0.0432),
         )
 
         return parachute

rocketpy/rocket/rocket.py

@@ -1433,7 +1433,16 @@ def add_free_form_fins(
         return fin_set
 
     def add_parachute(
-        self, name, cd_s, trigger, sampling_rate=100, lag=0, noise=(0, 0, 0)
+        self,
+        name,
+        cd_s,
+        trigger,
+        sampling_rate=100,
+        lag=0,
+        noise=(0, 0, 0),
+        parachute_radius=1.5,
+        parachute_height=None,
+        porosity=0.0432,
     ):
         """Creates a new parachute, storing its parameters such as
         opening delay, drag coefficients and trigger function.
@@ -1492,16 +1501,37 @@ def add_parachute(
             The values are used to add noise to the pressure signal which is
             passed to the trigger function. Default value is (0, 0, 0). Units
             are in pascal.
+        parachute_radius : float, optional
+            Length of the non-unique semi-axis (radius) of the inflated hemispheroid
+            parachute. Default value is 1.5.
+            Units are in meters.
+        parachute_height : float, optional
+            Length of the unique semi-axis (height) of the inflated hemispheroid
+            parachute. Default value is the radius of the parachute.
+            Units are in meters.
+        porosity : float, optional
+            Porosity of the parachute is the ratio of open space in the canopy to total
+            canopy area. Default value is 0.0432 (for consistency with previous versions).
 
         Returns
         -------
         parachute : Parachute
-            Parachute  containing trigger, sampling_rate, lag, cd_s, noise
-            and name. Furthermore, it stores clean_pressure_signal,
+            Parachute containing trigger, sampling_rate, lag, cd_s, noise, radius,
+            height, porosity and name. Furthermore, it stores clean_pressure_signal,
             noise_signal and noisyPressureSignal which are filled in during
             Flight simulation.
         """
-        parachute = Parachute(name, cd_s, trigger, sampling_rate, lag, noise)
+        parachute = Parachute(
+            name,
+            cd_s,
+            trigger,
+            sampling_rate,
+            lag,
+            noise,
+            parachute_radius,
+            parachute_height,
+            porosity,
+        )
         self.parachutes.append(parachute)
         return self.parachutes[-1]
 

rocketpy/simulation/flight.py

@@ -765,7 +765,18 @@ def __simulate(self, verbose):
                         callbacks = [
                             lambda self, parachute_cd_s=parachute.cd_s: setattr(
                                 self, "parachute_cd_s", parachute_cd_s
-                            )
+                            ),
+                            lambda self,
+                            parachute_radius=parachute.parachute_radius: setattr(
+                                self, "parachute_radius", parachute_radius
+                            ),
+                            lambda self,
+                            parachute_height=parachute.parachute_height: setattr(
+                                self, "parachute_height", parachute_height
+                            ),
+                            lambda self, parachute_porosity=parachute.porosity: setattr(
+                                self, "parachute_porosity", parachute_porosity
+                            ),
                         ]
                         self.flight_phases.add_phase(
                             node.t + parachute.lag,
@@ -1011,7 +1022,25 @@ def __simulate(self, verbose):
                                             lambda self,
                                             parachute_cd_s=parachute.cd_s: setattr(
                                                 self, "parachute_cd_s", parachute_cd_s
-                                            )
+                                            ),
+                                            lambda self,
+                                            parachute_radius=parachute.parachute_radius: setattr(
+                                                self,
+                                                "parachute_radius",
+                                                parachute_radius,
+                                            ),
+                                            lambda self,
+                                            parachute_height=parachute.parachute_height: setattr(
+                                                self,
+                                                "parachute_height",
+                                                parachute_height,
+                                            ),
+                                            lambda self,
+                                            parachute_porosity=parachute.porosity: setattr(
+                                                self,
+                                                "parachute_porosity",
+                                                parachute_porosity,
+                                            ),
                                         ]
                                         self.flight_phases.add_phase(
                                             overshootable_node.t + parachute.lag,
@@ -1965,22 +1994,28 @@ def u_dot_parachute(self, t, u, post_processing=False):
 
         # Get Parachute data
         cd_s = self.parachute_cd_s
+        parachute_radius = self.parachute_radius
+        parachute_height = self.parachute_height
+        porosity = self.parachute_porosity
 
         # Get the mass of the rocket
         mp = self.rocket.dry_mass
 
         # Define constants
-        ka = 1  # Added mass coefficient (depends on parachute's porosity)
-        R = 1.5  # Parachute radius
+        ka = 1.068 * (
+            1 - 1.465 * porosity - 0.25975 * porosity**2 + 1.2626 * porosity**3
+        )
         # to = 1.2
         # eta = 1
         # Rdot = (6 * R * (1 - eta) / (1.2**6)) * (
         #     (1 - eta) * t**5 + eta * (to**3) * (t**2)
         # )
         # Rdot = 0
 
+        # tf = 8 * nominal diameter / velocity at line stretch
+
         # Calculate added mass
-        ma = ka * rho * (4 / 3) * np.pi * R**3
+        ma = ka * rho * (4 / 3) * np.pi * parachute_radius**2 * parachute_height
 
         # Calculate freestream speed
         freestream_x = vx - wind_velocity_x
@@ -1991,12 +2026,12 @@ def u_dot_parachute(self, t, u, post_processing=False):
         # Determine drag force
         pseudo_drag = -0.5 * rho * cd_s * free_stream_speed
         # pseudo_drag = pseudo_drag - ka * rho * 4 * np.pi * (R**2) * Rdot
-        Dx = pseudo_drag * freestream_x
+        Dx = pseudo_drag * freestream_x  # add eta efficiency for wake
         Dy = pseudo_drag * freestream_y
         Dz = pseudo_drag * freestream_z
         ax = Dx / (mp + ma)
         ay = Dy / (mp + ma)
-        az = (Dz - 9.8 * mp) / (mp + ma)
+        az = (Dz - mp * self.env.gravity.get_value_opt(z)) / (mp + ma)
 
         # Add coriolis acceleration
         _, w_earth_y, w_earth_z = self.env.earth_rotation_vector

tests/unit/test_flight.py

@@ -252,7 +252,7 @@ def test_aerodynamic_moments(flight_calisto_custom_wind, flight_time, expected_v
         ("t_initial", (1.654150, 0.659142, -0.067103)),
         ("out_of_rail_time", (5.052628, 2.013361, -1.75370)),
         ("apogee_time", (2.321838, -1.613641, -0.962108)),
-        ("t_final", (-0.025792, 0.012030, 159.202481)),
+        ("t_final", (-0.025792, 0.012030, 159.061515)),
     ],
 )
 def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_values):

tests/unit/test_utilities.py

@@ -114,7 +114,7 @@ def test_fin_flutter_analysis(flight_calisto_custom_wind):
     assert np.isclose(flutter_mach(np.inf), 1.0048188594647927, atol=5e-3)
     assert np.isclose(safety_factor(0), 64.78797, atol=5e-3)
     assert np.isclose(safety_factor(10), 2.1948620401502072, atol=5e-3)
-    assert np.isclose(safety_factor(np.inf), 61.64222220697017, atol=5e-3)
+    assert np.isclose(safety_factor(np.inf), 61.669562809629035, atol=5e-3)
 
 
 def test_flutter_prints(flight_calisto_custom_wind):