ENH: Improve parachute geometric parametrization

rocketpy/rocket/parachute.py

@@ -102,6 +102,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,13 +157,27 @@ 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
+            Radius of the inflated parachute in meters. Default value is 1.5.
+            Units are in meters.
+        parachute_height : float, optional
+            Height of the inflated parachute in meters.
+            Default value is the radius parachute. Units are in meters.
+        porosity : float, optional
+            Porosity of the parachute material, which is a measure of how much air can
+            pass through the parachute material. Default value is 0.0432.
         """
         self.name = name
         self.cd_s = cd_s
         self.trigger = trigger
         self.sampling_rate = sampling_rate
         self.lag = lag
         self.noise = noise
+        self.parachute_radius = parachute_radius
+        if parachute_height is None:
+            parachute_height = parachute_radius
+        self.parachute_height = parachute_height
+        self.porosity = porosity
         self.noise_signal = [[-1e-6, np.random.normal(noise[0], noise[1])]]
         self.noisy_pressure_signal = []
         self.clean_pressure_signal = []
@@ -264,6 +281,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 +310,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/simulation/flight.py

@@ -767,7 +767,16 @@ 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,
@@ -1013,7 +1022,16 @@ 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,
@@ -1961,22 +1979,27 @@ 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
@@ -1987,12 +2010,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)
 
         if post_processing:
             self.__post_processed_variables.append(