Feature/cluster motor structure

rocketpy/motors/cluster_motor.py

@@ -0,0 +1,263 @@
+# pylint: disable=invalid-name
+import matplotlib.pyplot as plt
+import numpy as np
+from rocketpy import Function
+from rocketpy.motors import Motor
+
+
+class ClusterMotor(Motor):
+    """
+    A class representing a cluster of N identical motors arranged symmetrically.
+
+    This class aggregates the physical properties (thrust, mass, inertia) of
+    multiple motors using the Parallel Axis Theorem (Huygens-Steiner theorem).
+
+    Attributes
+    ----------
+    motor : SolidMotor
+        The single motor instance used in the cluster.
+    number : int
+        The number of motors in the cluster.
+    radius : float
+        The radial distance from the rocket's central axis to the center of each motor.
+    """
+
+    def __init__(self, motor, number, radius):
+        """
+        Initialize the ClusterMotor.
+
+        Parameters
+        ----------
+        motor : SolidMotor
+            The base motor to be clustered.
+        number : int
+            Number of motors. Must be >= 2.
+        radius : float
+            Distance from center of rocket to center of motor (m).
+        """
+        if not isinstance(number, int):
+            raise TypeError(f"number must be an int, got {type(number).__name__}")
+        if number < 2:
+            raise ValueError("number must be >= 2 for a ClusterMotor")
+        if not isinstance(radius, (int, float)):
+            raise TypeError(
+                f"radius must be a real number, got {type(radius).__name__}"
+            )
+        if radius < 0:
+            raise ValueError("radius must be non-negative")
+
+        self.motor = motor
+        self.number = number
+        self.radius = float(radius)
+        dry_inertia_cluster = self._calculate_dry_inertia()
+
+        # Use a thrust source scaled by the number of motors so that
+        # all thrust-derived quantities computed by the base Motor class
+        # correspond to the full cluster rather than a single motor.
+        scaled_thrust_source = motor.thrust * number
+
+        super().__init__(
+            thrust_source=scaled_thrust_source,
+            nozzle_radius=motor.nozzle_radius,
+            burn_time=motor.burn_time,
+            dry_mass=motor.dry_mass * number,
+            dry_inertia=dry_inertia_cluster,
+            center_of_dry_mass_position=motor.center_of_dry_mass_position,
+            coordinate_system_orientation=motor.coordinate_system_orientation,
+            interpolation_method="linear",
+        )
+
+        self._setup_grain_properties()
+        self._propellant_mass = self.motor.propellant_mass * self.number
+        self._propellant_initial_mass = self.number * self.motor.propellant_initial_mass
+        self._center_of_propellant_mass = self.motor.center_of_propellant_mass
+        self._evaluate_propellant_inertia()
+
+    def _evaluate_propellant_inertia(self):
+        """Calculates the dynamic inertia of the propellant using Steiner's theorem."""
+        Ixx_term1 = self.motor.propellant_I_11 * self.number
+        Ixx_term2 = self.motor.propellant_mass * (0.5 * self.number * self.radius**2)
+        self._propellant_I_11 = Ixx_term1 + Ixx_term2
+        self._propellant_I_22 = self._propellant_I_11
+
+        Izz_term1 = self.motor.propellant_I_33 * self.number
+        Izz_term2 = self.motor.propellant_mass * (self.number * self.radius**2)
+        self._propellant_I_33 = Izz_term1 + Izz_term2
+
+        zero_func = Function(0)
+        self._propellant_I_12 = zero_func
+        self._propellant_I_13 = zero_func
+        self._propellant_I_23 = zero_func
+
+    def _setup_grain_properties(self):
+        """Copies the grain properties from the base motor."""
+        self.throat_radius = self.motor.throat_radius
+        self.grain_number = self.motor.grain_number
+        self.grain_density = self.motor.grain_density
+        self.grain_outer_radius = self.motor.grain_outer_radius
+        self.grain_initial_inner_radius = self.motor.grain_initial_inner_radius
+        self.grain_initial_height = self.motor.grain_initial_height
+        self.grains_center_of_mass_position = self.motor.grains_center_of_mass_position
+
+    @property
+    def thrust(self):
+        return self._thrust
+
+    @thrust.setter
+    def thrust(self, value):
+        self._thrust = value
+
+    @property
+    def propellant_mass(self):
+        return self._propellant_mass
+
+    @propellant_mass.setter
+    def propellant_mass(self, value):
+        self._propellant_mass = value
+
+    @property
+    def propellant_initial_mass(self):
+        return self._propellant_initial_mass
+
+    @propellant_initial_mass.setter
+    def propellant_initial_mass(self, value):
+        self._propellant_initial_mass = value
+
+    @property
+    def center_of_propellant_mass(self):
+        return self._center_of_propellant_mass
+
+    @center_of_propellant_mass.setter
+    def center_of_propellant_mass(self, value):
+        self._center_of_propellant_mass = value
+
+    @property
+    def propellant_I_11(self):
+        return self._propellant_I_11
+
+    @propellant_I_11.setter
+    def propellant_I_11(self, value):
+        self._propellant_I_11 = value
+
+    @property
+    def propellant_I_22(self):
+        return self._propellant_I_22
+
+    @propellant_I_22.setter
+    def propellant_I_22(self, value):
+        self._propellant_I_22 = value
+
+    @property
+    def propellant_I_33(self):
+        return self._propellant_I_33
+
+    @propellant_I_33.setter
+    def propellant_I_33(self, value):
+        self._propellant_I_33 = value
+
+    @property
+    def propellant_I_12(self):
+        return self._propellant_I_12
+
+    @propellant_I_12.setter
+    def propellant_I_12(self, value):
+        self._propellant_I_12 = value
+
+    @property
+    def propellant_I_13(self):
+        return self._propellant_I_13
+
+    @propellant_I_13.setter
+    def propellant_I_13(self, value):
+        self._propellant_I_13 = value
+
+    @property
+    def propellant_I_23(self):
+        return self._propellant_I_23
+
+    @propellant_I_23.setter
+    def propellant_I_23(self, value):
+        self._propellant_I_23 = value
+
+    @property
+    def exhaust_velocity(self):
+        return self.motor.exhaust_velocity
+
+    def _calculate_dry_inertia(self):
+        Ixx_loc = self.motor.dry_I_11
+        Iyy_loc = self.motor.dry_I_22
+        Izz_loc = self.motor.dry_I_33
+        m_dry = self.motor.dry_mass
+
+        Izz_cluster = self.number * Izz_loc + self.number * m_dry * (self.radius**2)
+        Ixx_cluster = self.number * Ixx_loc + (self.number / 2) * m_dry * (
+            self.radius**2
+        )
+        Iyy_cluster = self.number * Iyy_loc + (self.number / 2) * m_dry * (
+            self.radius**2
+        )
+
+        return (Ixx_cluster, Iyy_cluster, Izz_cluster)
+
+    def info(self, *args, **kwargs):
+        print("Cluster Configuration:")
+        print(f" - Motors: {self.number} x {type(self.motor).__name__}")
+        print(f" - Radial Distance: {self.radius} m")
+        return self.motor.info(*args, **kwargs)
+
+    def draw_cluster_layout(self, rocket_radius=None, show=True):
+        """Draw the geometric layout of the clustered motors."""
+        fig, ax = plt.subplots(figsize=(6, 6))
+        ax.plot(0, 0, "k+", markersize=10, label="Central axis")
+        if rocket_radius:
+            rocket_tube = plt.Circle(
+                (0, 0),
+                rocket_radius,
+                color="black",
+                fill=False,
+                linestyle="--",
+                linewidth=2,
+                label="Rocket",
+            )
+            ax.add_patch(rocket_tube)
+            limit = rocket_radius * 1.2
+        else:
+            limit = self.radius * 2
+        self._draw_engines(ax)
+        ax.set_aspect("equal", "box")
+        ax.set_xlim(-limit, limit)
+        ax.set_ylim(-limit, limit)
+        ax.set_xlabel("Position X (m)")
+        ax.set_ylabel("Position Y (m)")
+        ax.set_title(f"Cluster Configuration : {self.number} engines")
+        ax.grid(True, linestyle=":", alpha=0.6)
+        ax.legend(loc="upper right")
+        if show:
+            plt.show()
+        return fig, ax
+
+    def _draw_engines(self, ax):
+        """Draws the individual engines of the cluster."""
+        motor_outer_radius = self.grain_outer_radius
+        angles = np.linspace(0, 2 * np.pi, self.number, endpoint=False)
+
+        for i, angle in enumerate(angles):
+            x = self.radius * np.cos(angle)
+            y = self.radius * np.sin(angle)
+            motor_circle = plt.Circle(
+                (x, y),
+                motor_outer_radius,
+                color="red",
+                alpha=0.5,
+                label="Engine" if i == 0 else "",
+            )
+            ax.add_patch(motor_circle)
+            ax.text(
+                x,
+                y,
+                str(i + 1),
+                color="white",
+                ha="center",
+                va="center",
+                fontweight="bold",
+            )

rocketpy/plots/rocket_plots.py

@@ -1,7 +1,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from rocketpy.motors import EmptyMotor, HybridMotor, LiquidMotor, SolidMotor
+from rocketpy.motors import HybridMotor, LiquidMotor, SolidMotor
 from rocketpy.rocket.aero_surface import Fins, NoseCone, Tail
 from rocketpy.rocket.aero_surface.generic_surface import GenericSurface
 
@@ -420,57 +420,84 @@ def _draw_tubes(self, ax, drawn_surfaces, vis_args):
     def _draw_motor(self, last_radius, last_x, ax, vis_args):
         """Draws the motor from motor patches"""
         total_csys = self.rocket._csys * self.rocket.motor._csys
+        is_cluster = hasattr(self.rocket.motor, "number")
+        base_motor = self.rocket.motor.motor if is_cluster else self.rocket.motor
+
+        if is_cluster:
+            angles = np.linspace(0, 2 * np.pi, self.rocket.motor.number, endpoint=False)
+            y_offsets = self.rocket.motor.radius * np.cos(angles)
+        else:
+            y_offsets = [0]
         nozzle_position = (
-            self.rocket.motor_position + self.rocket.motor.nozzle_position * total_csys
+            self.rocket.motor_position + base_motor.nozzle_position * total_csys
         )
-
         # Get motor patches translated to the correct position
         motor_patches = self._generate_motor_patches(total_csys, ax)
-
         # Draw patches
-        if not isinstance(self.rocket.motor, EmptyMotor):
-            # Add nozzle last so it is in front of the other patches
-            nozzle = self.rocket.motor.plots._generate_nozzle(
-                translate=(nozzle_position, 0), csys=self.rocket._csys
-            )
-            motor_patches += [nozzle]
+        if type(self.rocket.motor).__name__ != "EmptyMotor":
+            for y_off in y_offsets:
+                nozzle = base_motor.plots._generate_nozzle(
+                    translate=(nozzle_position, y_off), csys=self.rocket._csys
+                )
+                if y_off != y_offsets[0]:
+                    nozzle.set_label("_nolegend_")
+                motor_patches.append(nozzle)
 
-            outline = self.rocket.motor.plots._generate_motor_region(
+            outline = base_motor.plots._generate_motor_region(
                 list_of_patches=motor_patches
             )
-            # add outline first so it is behind the other patches
-            ax.add_patch(outline)
+            if not is_cluster:
+                ax.add_patch(outline)
+
             for patch in motor_patches:
+                if is_cluster:
+                    patch.set_alpha(0.6)
                 ax.add_patch(patch)
-
         self._draw_nozzle_tube(last_radius, last_x, nozzle_position, ax, vis_args)
 
-    def _generate_motor_patches(self, total_csys, ax):  # pylint: disable=unused-argument
+    def _generate_motor_patches(
+        self, total_csys, ax
+    ):  # pylint: disable=unused-argument
         """Generates motor patches for drawing"""
         motor_patches = []
 
-        if isinstance(self.rocket.motor, SolidMotor):
+        is_cluster = hasattr(self.rocket.motor, "number")
+        base_motor = self.rocket.motor.motor if is_cluster else self.rocket.motor
+
+        if isinstance(base_motor, SolidMotor):
+            y_offsets = (
+                self.rocket.motor.radius
+                * np.cos(
+                    np.linspace(0, 2 * np.pi, self.rocket.motor.number, endpoint=False)
+                )
+                if is_cluster
+                else [0]
+            )
             grains_cm_position = (
                 self.rocket.motor_position
-                + self.rocket.motor.grains_center_of_mass_position * total_csys
-            )
-            ax.scatter(
-                grains_cm_position,
-                0,
-                color="brown",
-                label="Grains Center of Mass",
-                s=8,
-                zorder=10,
+                + base_motor.grains_center_of_mass_position * total_csys
             )
+            for y_off in y_offsets:
+                ax.scatter(
+                    grains_cm_position,
+                    y_off,
+                    color="brown",
+                    label="Grains Center of Mass" if y_off == y_offsets[0] else "",
+                    s=8,
+                    zorder=10,
+                )
 
-            chamber = self.rocket.motor.plots._generate_combustion_chamber(
-                translate=(grains_cm_position, 0), label=None
-            )
-            grains = self.rocket.motor.plots._generate_grains(
-                translate=(grains_cm_position, 0)
-            )
+                chamber = base_motor.plots._generate_combustion_chamber(
+                    translate=(grains_cm_position, y_off), label=None
+                )
+                grains = base_motor.plots._generate_grains(
+                    translate=(grains_cm_position, y_off)
+                )
+                if y_off != y_offsets[0]:
+                    for grain in grains:
+                        grain.set_label("_nolegend_")
 
-            motor_patches += [chamber, *grains]
+                motor_patches += [chamber, *grains]
 
         elif isinstance(self.rocket.motor, HybridMotor):
             grains_cm_position = (