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.radius : float
+        Length of the non-unique semi-axis (radius) of the inflated hemispheroid
+        parachute in meters.
+    Parachute.height : float, None
+        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),
+        radius=1.5,
+        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.
+        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.
+        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,15 @@ def __init__(
         self.clean_pressure_signal_function = Function(0)
         self.noisy_pressure_signal_function = Function(0)
         self.noise_signal_function = Function(0)
+        self.radius = radius
+        self.height = height or radius
+        self.porosity = porosity
+        self.ka = 1.068 * (
+            1
+            - 1.465 * self.porosity
+            - 0.25975 * self.porosity**2
+            + 1.2626 * self.porosity**3
+        )
 
         alpha, beta = self.noise_corr
         self.noise_function = lambda: alpha * self.noise_signal[-1][
@@ -264,6 +296,9 @@ def to_dict(self, include_outputs=False):
             "sampling_rate": self.sampling_rate,
             "lag": self.lag,
             "noise": self.noise,
+            "radius": self.radius,
+            "height": self.height,
+            "porosity": self.porosity,
         }
 
         if include_outputs:
@@ -290,6 +325,9 @@ def from_dict(cls, data):
             sampling_rate=data["sampling_rate"],
             lag=data["lag"],
             noise=data["noise"],
+            radius=data.get("radius", 1.5),
+            height=data.get("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),
+        radius=1.5,
+        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.
+        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.
+        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,
+            radius,
+            height,
+            porosity,
+        )
         self.parachutes.append(parachute)
         return self.parachutes[-1]
 

rocketpy/simulation/flight.py

@@ -765,7 +765,21 @@ 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
+                            ),
+                            lambda self, ka=parachute.ka: setattr(
+                                self, "parachute_ka", ka
+                            ),
                         ]
                         self.flight_phases.add_phase(
                             node.t + parachute.lag,
@@ -1011,7 +1025,28 @@ 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.radius: setattr(
+                                                self,
+                                                "parachute_radius",
+                                                parachute_radius,
+                                            ),
+                                            lambda self,
+                                            parachute_height=parachute.height: setattr(
+                                                self,
+                                                "parachute_height",
+                                                parachute_height,
+                                            ),
+                                            lambda self,
+                                            parachute_porosity=parachute.porosity: setattr(
+                                                self,
+                                                "parachute_porosity",
+                                                parachute_porosity,
+                                            ),
+                                            lambda self, ka=parachute.ka: setattr(
+                                                self, "ka", ka
+                                            ),
                                         ]
                                         self.flight_phases.add_phase(
                                             overshootable_node.t + parachute.lag,
@@ -1963,24 +1998,27 @@ def u_dot_parachute(self, t, u, post_processing=False):
         wind_velocity_x = self.env.wind_velocity_x.get_value_opt(z)
         wind_velocity_y = self.env.wind_velocity_y.get_value_opt(z)
 
-        # Get Parachute data
-        cd_s = self.parachute_cd_s
-
         # 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
         # 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 = (
+            self.ka
+            * rho
+            * (4 / 3)
+            * np.pi
+            * self.parachute_radius**2
+            * self.parachute_height
+        )
 
         # Calculate freestream speed
         freestream_x = vx - wind_velocity_x
@@ -1989,14 +2027,14 @@ def u_dot_parachute(self, t, u, post_processing=False):
         free_stream_speed = (freestream_x**2 + freestream_y**2 + freestream_z**2) ** 0.5
 
         # Determine drag force
-        pseudo_drag = -0.5 * rho * cd_s * free_stream_speed
+        pseudo_drag = -0.5 * rho * self.parachute_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

rocketpy/stochastic/stochastic_parachute.py

@@ -29,6 +29,12 @@ class StochasticParachute(StochasticModel):
         time-correlation).
     name : list[str]
         List with the name of the parachute object. This cannot be randomized.
+    radius : tuple, list, int, float
+        Radius of the parachute in meters.
+    height : tuple, list, int, float
+        Height of the parachute in meters.
+    porosity : tuple, list, int, float
+        Porosity of the parachute.
     """
 
     def __init__(
@@ -39,6 +45,9 @@ def __init__(
         sampling_rate=None,
         lag=None,
         noise=None,
+        radius=None,
+        height=None,
+        porosity=None,
     ):
         """Initializes the Stochastic Parachute class.
 
@@ -63,13 +72,22 @@ def __init__(
         noise : list
             List of tuples in the form of (mean, standard deviation,
             time-correlation).
+        radius : tuple, list, int, float
+            Radius of the parachute in meters.
+        height : tuple, list, int, float
+            Height of the parachute in meters.
+        porosity : tuple, list, int, float
+            Porosity of the parachute.
         """
         self.parachute = parachute
         self.cd_s = cd_s
         self.trigger = trigger
         self.sampling_rate = sampling_rate
         self.lag = lag
         self.noise = noise
+        self.radius = radius
+        self.height = height
+        self.porosity = porosity
 
         self._validate_trigger(trigger)
         self._validate_noise(noise)
@@ -81,6 +99,9 @@ def __init__(
             lag=lag,
             noise=noise,
             name=None,
+            radius=radius,
+            height=height,
+            porosity=porosity,
         )
 
     def _validate_trigger(self, trigger):

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):