BUG: Fix hard-coded radius value for parachute added mass calculation

CHANGELOG.md

@@ -58,6 +58,7 @@ Attention: The newest changes should be on top -->
 
 ### Fixed
 
+- BUG: Fix hard-coded radius value for parachute added mass calculation [#889](https://github.com/RocketPy-Team/RocketPy/pull/889)
 - DOC: Fix documentation build [#908](https://github.com/RocketPy-Team/RocketPy/pull/908)
 - BUG: energy_data plot not working for 3 dof sims [[#906](https://github.com/RocketPy-Team/RocketPy/issues/906)]
 - BUG: Fix CSV column header spacing in FlightDataExporter [#864](https://github.com/RocketPy-Team/RocketPy/issues/864)

rocketpy/rocket/parachute.py

@@ -92,17 +92,25 @@ class Parachute:
         Function of noisy_pressure_signal.
     Parachute.clean_pressure_signal_function : Function
         Function of clean_pressure_signal.
+    Parachute.drag_coefficient : float
+        Drag coefficient of the inflated canopy shape, used only when
+        ``radius`` is not provided to estimate the parachute radius from
+        ``cd_s``: ``R = sqrt(cd_s / (drag_coefficient * pi))``. Typical
+        values: 1.4 for hemispherical canopies (default), 0.75 for flat
+        circular canopies, 1.5 for extended-skirt canopies.
     Parachute.radius : float
         Length of the non-unique semi-axis (radius) of the inflated hemispheroid
-        parachute in meters.
-    Parachute.height : float, None
+        parachute in meters. If not provided at construction time, it is
+        estimated from ``cd_s`` and ``drag_coefficient``.
+    Parachute.height : float
         Length of the unique semi-axis (height) of the inflated hemispheroid
         parachute in meters.
     Parachute.porosity : float
-        Geometric porosity of the canopy (ratio of open area to total canopy area),
-        in [0, 1]. Affects only the added-mass scaling during descent; it does
-        not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass
-        of 1.0 (“neutral” behavior).
+        Geometric porosity of the canopy (ratio of open area to total canopy
+        area), in [0, 1]. Affects only the added-mass scaling during descent;
+        it does not change ``cd_s`` (drag). The default value of 0.0432 is
+        chosen so that the resulting ``added_mass_coefficient`` equals
+        approximately 1.0 ("neutral" added-mass behavior).
     Parachute.added_mass_coefficient : float
         Coefficient used to calculate the added-mass due to dragged air. It is
         calculated from the porosity of the parachute.
@@ -116,9 +124,10 @@ def __init__(
         sampling_rate,
         lag=0,
         noise=(0, 0, 0),
-        radius=1.5,
+        radius=None,
         height=None,
         porosity=0.0432,
+        drag_coefficient=1.4,
     ):
         """Initializes Parachute class.
 
@@ -172,25 +181,83 @@ def __init__(
             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.
+            Length of the non-unique semi-axis (radius) of the inflated
+            hemispheroid parachute. If not provided, it is estimated from
+            ``cd_s`` and ``drag_coefficient`` using:
+            ``radius = sqrt(cd_s / (drag_coefficient * pi))``.
             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
-            Geometric porosity of the canopy (ratio of open area to total canopy area),
-            in [0, 1]. Affects only the added-mass scaling during descent; it does
-            not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass
-            of 1.0 (“neutral” behavior).
+            Geometric porosity of the canopy (ratio of open area to total
+            canopy area), in [0, 1]. Affects only the added-mass scaling
+            during descent; it does not change ``cd_s`` (drag). The default
+            value of 0.0432 is chosen so that the resulting
+            ``added_mass_coefficient`` equals approximately 1.0 ("neutral"
+            added-mass behavior).
+        drag_coefficient : float, optional
+            Drag coefficient of the inflated canopy shape, used only when
+            ``radius`` is not provided. It relates the aerodynamic ``cd_s``
+            to the physical canopy area via
+            ``cd_s = drag_coefficient * pi * radius**2``. Typical values:
+
+            - **1.4** — hemispherical canopy (default, NASA SP-8066)
+            - **0.75** — flat circular canopy
+            - **1.5** — extended-skirt canopy
+
+            Has no effect when ``radius`` is explicitly provided.
         """
+
+        # Save arguments as attributes
         self.name = name
         self.cd_s = cd_s
         self.trigger = trigger
         self.sampling_rate = sampling_rate
         self.lag = lag
         self.noise = noise
+        self.drag_coefficient = drag_coefficient
+        self.porosity = porosity
+
+        # Initialize derived attributes
+        self.radius = self.__resolve_radius(radius, cd_s, drag_coefficient)
+        self.height = self.__resolve_height(height, self.radius)
+        self.added_mass_coefficient = self.__compute_added_mass_coefficient(
+            self.porosity
+        )
+        self.__init_noise(noise)
+        self.__evaluate_trigger_function(trigger)
+
+        # Prints and plots
+        self.prints = _ParachutePrints(self)
+
+    def __resolve_radius(self, radius, cd_s, drag_coefficient):
+        """Resolves parachute radius from input or aerodynamic relation."""
+        if radius is not None:
+            return radius
+
+        # cd_s = Cd * S = Cd * pi * R^2  =>  R = sqrt(cd_s / (Cd * pi))
+        return np.sqrt(cd_s / (drag_coefficient * np.pi))
+
+    def __resolve_height(self, height, radius):
+        """Resolves parachute height defaulting to radius when not provided."""
+        return height or radius
+
+    def __compute_added_mass_coefficient(self, porosity):
+        """Computes the added-mass coefficient from canopy porosity."""
+        return 1.068 * (
+            1 - 1.465 * porosity - 0.25975 * porosity**2 + 1.2626 * porosity**3
+        )
+
+    def __init_noise(self, noise):
+        """Initializes all noise-related attributes.
+
+        Parameters
+        ----------
+        noise : tuple, list
+            List in the format (mean, standard deviation, time-correlation).
+        """
         self.noise_signal = [[-1e-6, np.random.normal(noise[0], noise[1])]]
         self.noisy_pressure_signal = []
         self.clean_pressure_signal = []
@@ -200,32 +267,19 @@ 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.added_mass_coefficient = 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][1]
             + beta * np.random.normal(noise[0], noise[1])
         )
 
-        self.prints = _ParachutePrints(self)
-
-        self.__evaluate_trigger_function(trigger)
-
     def __evaluate_trigger_function(self, trigger):
         """This is used to set the triggerfunc attribute that will be used to
         interact with the Flight class.
         """
         # pylint: disable=unused-argument, function-redefined
-        # The parachute is deployed by a custom function
+
+        # Case 1: The parachute is deployed by a custom function
         if callable(trigger):
             # work around for having added sensors to parachute triggers
             # to avoid breaking changes
@@ -238,26 +292,29 @@ def triggerfunc(p, h, y, sensors):
 
             self.triggerfunc = triggerfunc
 
+        # Case 2: The parachute is deployed at a given height
         elif isinstance(trigger, (int, float)):
             # The parachute is deployed at a given height
-            def triggerfunc(p, h, y, sensors):  # pylint: disable=unused-argument
+            def triggerfunc(p, h, y, sensors):
                 # p = pressure considering parachute noise signal
                 # h = height above ground level considering parachute noise signal
                 # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]
                 return y[5] < 0 and h < trigger
 
             self.triggerfunc = triggerfunc
 
+        # Case 3: The parachute is deployed at apogee
         elif trigger.lower() == "apogee":
             # The parachute is deployed at apogee
-            def triggerfunc(p, h, y, sensors):  # pylint: disable=unused-argument
+            def triggerfunc(p, h, y, sensors):
                 # p = pressure considering parachute noise signal
                 # h = height above ground level considering parachute noise signal
                 # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3]
                 return y[5] < 0
 
             self.triggerfunc = triggerfunc
 
+        # Case 4: Invalid trigger input
         else:
             raise ValueError(
                 f"Unable to set the trigger function for parachute '{self.name}'. "
@@ -289,7 +346,7 @@ def info(self):
     def all_info(self):
         """Prints all information about the Parachute class."""
         self.info()
-        # self.plots.all() # Parachutes still doesn't have plots
+        # self.plots.all() # TODO: Parachutes still doesn't have plots
 
     def to_dict(self, **kwargs):
         allow_pickle = kwargs.get("allow_pickle", True)
@@ -309,6 +366,7 @@ def to_dict(self, **kwargs):
             "lag": self.lag,
             "noise": self.noise,
             "radius": self.radius,
+            "drag_coefficient": self.drag_coefficient,
             "height": self.height,
             "porosity": self.porosity,
         }
@@ -341,7 +399,8 @@ def from_dict(cls, data):
             sampling_rate=data["sampling_rate"],
             lag=data["lag"],
             noise=data["noise"],
-            radius=data.get("radius", 1.5),
+            radius=data.get("radius", None),
+            drag_coefficient=data.get("drag_coefficient", 1.4),
             height=data.get("height", None),
             porosity=data.get("porosity", 0.0432),
         )

rocketpy/rocket/rocket.py

@@ -1502,9 +1502,10 @@ def add_parachute(
         sampling_rate=100,
         lag=0,
         noise=(0, 0, 0),
-        radius=1.5,
+        radius=None,
         height=None,
         porosity=0.0432,
+        drag_coefficient=1.4,
     ):
         """Creates a new parachute, storing its parameters such as
         opening delay, drag coefficients and trigger function.
@@ -1564,26 +1565,34 @@ def add_parachute(
             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.
+            Length of the non-unique semi-axis (radius) of the inflated
+            hemispheroid parachute. If not provided, it is estimated from
+            `cd_s` and `drag_coefficient` using:
+            `radius = sqrt(cd_s / (drag_coefficient * pi))`.
             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
-            Geometric porosity of the canopy (ratio of open area to total canopy area),
-            in [0, 1]. Affects only the added-mass scaling during descent; it does
-            not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass
-            of 1.0 (“neutral” behavior).
+            Geometric porosity of the canopy (ratio of open area to total
+            canopy area), in [0, 1]. Affects only the added-mass scaling
+            during descent; it does not change `cd_s` (drag). The default
+            value of 0.0432 yields an `added_mass_coefficient` of
+            approximately 1.0 ("neutral" added-mass behavior).
+        drag_coefficient : float, optional
+            Drag coefficient of the inflated canopy shape, used only when
+            `radius` is not provided. Typical values: 1.4 for hemispherical
+            canopies (default), 0.75 for flat circular canopies, 1.5 for
+            extended-skirt canopies. Has no effect when `radius` is given.
 
         Returns
         -------
         parachute : Parachute
-            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 containing trigger, sampling_rate, lag, cd_s, noise,
+            radius, drag_coefficient, height, porosity and name. Furthermore,
+            it stores clean_pressure_signal, noise_signal and
+            noisyPressureSignal which are filled in during Flight simulation.
         """
         parachute = Parachute(
             name,
@@ -1595,6 +1604,7 @@ def add_parachute(
             radius,
             height,
             porosity,
+            drag_coefficient,
         )
         self.parachutes.append(parachute)
         return self.parachutes[-1]

rocketpy/stochastic/stochastic_parachute.py

@@ -31,6 +31,9 @@ class StochasticParachute(StochasticModel):
         List with the name of the parachute object. This cannot be randomized.
     radius : tuple, list, int, float
         Radius of the parachute in meters.
+    drag_coefficient : tuple, list, int, float
+        Drag coefficient of the inflated canopy shape, used only when
+        ``radius`` is not provided.
     height : tuple, list, int, float
         Height of the parachute in meters.
     porosity : tuple, list, int, float
@@ -46,6 +49,7 @@ def __init__(
         lag=None,
         noise=None,
         radius=None,
+        drag_coefficient=None,
         height=None,
         porosity=None,
     ):
@@ -74,6 +78,9 @@ def __init__(
             time-correlation).
         radius : tuple, list, int, float
             Radius of the parachute in meters.
+        drag_coefficient : tuple, list, int, float
+            Drag coefficient of the inflated canopy shape, used only when
+            ``radius`` is not provided.
         height : tuple, list, int, float
             Height of the parachute in meters.
         porosity : tuple, list, int, float
@@ -86,6 +93,7 @@ def __init__(
         self.lag = lag
         self.noise = noise
         self.radius = radius
+        self.drag_coefficient = drag_coefficient
         self.height = height
         self.porosity = porosity
 
@@ -100,6 +108,7 @@ def __init__(
             noise=noise,
             name=None,
             radius=radius,
+            drag_coefficient=drag_coefficient,
             height=height,
             porosity=porosity,
         )