diff --git a/rocketpy/motors/cluster_motor.py b/rocketpy/motors/cluster_motor.py
new file mode 100644
index 000000000..6a71dd40c
--- /dev/null
+++ b/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",
+            )
diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py
index e208c775f..a01a62dfa 100644
--- a/rocketpy/plots/rocket_plots.py
+++ b/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 = (
diff --git a/tests/integration/motors/test_cluster_motor.py b/tests/integration/motors/test_cluster_motor.py
new file mode 100644
index 000000000..4c7ce20d0
--- /dev/null
+++ b/tests/integration/motors/test_cluster_motor.py
@@ -0,0 +1,149 @@
+import pytest
+import numpy as np
+from rocketpy import SolidMotor, Function
+from rocketpy.motors.cluster_motor import ClusterMotor
+
+
+@pytest.fixture
+def base_motor():
+    """
+    Creates a simplified SolidMotor for testing purposes.
+    Properties:
+    - Constant Thrust: 1000 N
+    - Burn time: 5 s
+    - Dry mass: 10 kg
+    - Dry Inertia: (1.0, 1.0, 0.1)
+    """
+    thrust_curve = Function(lambda t: 1000 if t < 5 else 0, "Time (s)", "Thrust (N)")
+
+    return SolidMotor(
+        thrust_source=thrust_curve,
+        burn_time=5,
+        dry_mass=10.0,
+        dry_inertia=(1.0, 1.0, 0.1),  # Ixx, Iyy, Izz
+        grain_number=1,
+        grain_density=1000,
+        grain_outer_radius=0.05,
+        grain_initial_inner_radius=0.02,
+        grain_initial_height=0.5,
+        coordinate_system_orientation="nozzle_to_combustion_chamber",
+        nozzle_radius=0.02,
+        grain_separation=0.001,
+        grains_center_of_mass_position=0.25,
+        center_of_dry_mass_position=0.25,
+    )
+
+
+def test_cluster_initialization(base_motor):
+    """
+    Tests if the ClusterMotor initializes basic attributes correctly.
+    """
+    N = 3
+    R = 0.5
+    cluster = ClusterMotor(motor=base_motor, number=N, radius=R)
+
+    assert cluster.number == N
+    assert cluster.radius == R
+    assert cluster.grain_outer_radius == base_motor.grain_outer_radius
+
+
+def test_cluster_mass_and_thrust_scaling(base_motor):
+    """
+    Tests if scalar and derived properties are correctly multiplied by N and that functional properties preserve their Function behavior
+    """
+    N = 4
+    R = 0.2
+    cluster = ClusterMotor(motor=base_motor, number=N, radius=R)
+
+    # 1. Check Thrust Scaling
+    # Thrust at t=1 should be N * single_motor_thrust
+    assert np.isclose(cluster.thrust(1), base_motor.thrust(1) * N)
+
+    # 2. Check Dry Mass Scaling
+    assert np.isclose(cluster.dry_mass, base_motor.dry_mass * N)
+
+    # 3. Check Propellant Mass Scaling
+    assert np.isclose(cluster.propellant_mass(0), base_motor.propellant_mass(0) * N)
+    assert np.isclose(cluster.total_impulse, base_motor.total_impulse * N)
+    assert np.isclose(cluster.average_thrust, base_motor.average_thrust * N)
+
+
+def test_cluster_dry_inertia_steiner_theorem(base_motor):
+    """
+    Tests the implementation of the Parallel Axis Theorem (Huygens-Steiner)
+    for the static (dry) mass of the cluster.
+
+    Theoretical Formulas:
+    I_zz_cluster = N * I_zz_local + N * m * R^2
+    I_xx_cluster = N * I_xx_local + (N/2) * m * R^2  (Radial symmetry approximation)
+    """
+    N = 3
+    R = 1.0  # 1 meter radius for simpler checking
+    cluster = ClusterMotor(motor=base_motor, number=N, radius=R)
+
+    m_dry = base_motor.dry_mass
+    Ixx_loc = base_motor.dry_I_11
+    Izz_loc = base_motor.dry_I_33
+
+    # Expected Izz (Longitudinal / Roll)
+    expected_Izz = N * Izz_loc + N * m_dry * (R**2)
+
+    # Expected Ixx/Iyy (Transverse / Pitch / Yaw)
+    expected_I_trans = N * Ixx_loc + (N / 2) * m_dry * (R**2)
+
+    assert np.isclose(cluster.dry_I_33, expected_Izz)
+    assert np.isclose(cluster.dry_I_11, expected_I_trans)
+    assert np.isclose(cluster.dry_I_22, expected_I_trans)
+
+
+def test_cluster_propellant_inertia_dynamic(base_motor):
+    """
+    Tests if the Steiner theorem is correctly applied dynamically
+    to the changing propellant mass over time.
+    """
+    N = 2
+    R = 0.5
+    cluster = ClusterMotor(motor=base_motor, number=N, radius=R)
+
+    t = 0  # Check at t=0
+
+    m_prop = base_motor.propellant_mass(t)
+    Ixx_prop_loc = base_motor.propellant_I_11(t)
+    Izz_prop_loc = base_motor.propellant_I_33(t)
+
+    # Expected Dynamic Ixx
+    # Ixx_term1 (Local rotation) + Ixx_term2 (Parallel axis offset)
+    expected_Ixx = (Ixx_prop_loc * N) + (m_prop * 0.5 * N * R**2)
+
+    # Expected Dynamic Izz
+    expected_Izz = (Izz_prop_loc * N) + (m_prop * N * R**2)
+
+    assert np.isclose(cluster.propellant_I_11(t), expected_Ixx)
+    assert np.isclose(cluster.propellant_I_33(t), expected_Izz)
+
+def test_cluster_invalid_inputs(base_motor):
+    """Tests if the validation raises errors for bad inputs."""
+    with pytest.raises(ValueError):
+        ClusterMotor(motor=base_motor, number=1, radius=0.5) # N < 2
+    with pytest.raises(ValueError):
+        ClusterMotor(motor=base_motor, number=2, radius=-1.0) # Radius < 0
+    with pytest.raises(TypeError):
+        ClusterMotor(motor=base_motor, number="two", radius=0.5) # N is string
+
+def test_cluster_methods_and_setters(base_motor):
+    """Touches the display methods and setters to ensure coverage."""
+    cluster = ClusterMotor(motor=base_motor, number=2, radius=0.5)
+    
+    # 1. Touch the info method
+    cluster.info()
+    
+    # 2. Touch the draw method (without showing the plot to avoid blocking tests)
+    cluster.draw_cluster_layout(show=False)
+    cluster.draw_cluster_layout(rocket_radius=0.1, show=False)
+    
+    # 3. Touch a few setters
+    cluster.propellant_mass = 50.0
+    assert cluster.propellant_mass == 50.0
+    
+    cluster.propellant_I_11 = 2.0
+    assert cluster.propellant_I_11 == 2.0
\ No newline at end of file