ENH: Improve Stability Margin with Angle of Attack Dependency

rocketpy/rocket/rocket.py

@@ -671,6 +671,91 @@ def evaluate_stability_margin(self):
         )
         return self.stability_margin
 
+    def get_cp_position_from_alpha(self, alpha, mach):
+        """Computes the center of pressure position as a function of the angle
+        of attack and Mach number. This method uses the actual lift coefficient
+        CN(α) instead of the lift coefficient derivative CNα, providing a more
+        accurate CP position that varies with angle of attack.
+
+        Parameters
+        ----------
+        alpha : float
+            Angle of attack in radians.
+        mach : float
+            Mach number.
+
+        Returns
+        -------
+        float
+            Center of pressure position relative to the user-defined rocket
+            reference system, in meters.
+
+        Notes
+        -----
+        The center of pressure is calculated as:
+            CP = Σ(CN(α) × X_cp) / Σ(CN(α))
+        where CN(α) is the lift coefficient as a function of angle of attack
+        and X_cp is the center of pressure position of each aerodynamic surface.
+        """
+        # Threshold for numerical zero checks
+        epsilon = 1e-10
+
+        # Handle edge case where alpha is effectively zero
+        if abs(alpha) < epsilon:
+            return self.cp_position.get_value_opt(mach)
+
+        total_cn = 0.0
+        weighted_cp = 0.0
+
+        for aero_surface, position in self.aerodynamic_surfaces:
+            if isinstance(aero_surface, GenericSurface):
+                continue
+
+            ref_factor = (aero_surface.rocket_radius / self.radius) ** 2
+            cn = ref_factor * aero_surface.cl.get_value_opt(alpha, mach)
+            surface_cp = position.z - self._csys * aero_surface.cpz
+
+            total_cn += cn
+            weighted_cp += cn * surface_cp
+
+        if abs(total_cn) < epsilon:
+            return self.cp_position.get_value_opt(mach)
+
+        return weighted_cp / total_cn
+
+    def get_stability_margin_from_alpha(self, alpha, mach, time):
+        """Computes the stability margin using the angle of attack-dependent
+        center of pressure position.
+
+        Parameters
+        ----------
+        alpha : float
+            Angle of attack in radians.
+        mach : float
+            Mach number.
+        time : float
+            Time in seconds.
+
+        Returns
+        -------
+        float
+            Stability margin in calibers. A positive value indicates the rocket
+            is stable (center of pressure is behind the center of mass).
+
+        Notes
+        -----
+        The stability margin is calculated as:
+            SM = (CG - CP(α)) / d
+        where CG is the center of gravity, CP(α) is the angle of attack-dependent
+        center of pressure, and d is the rocket diameter.
+        """
+        cp_position = self.get_cp_position_from_alpha(alpha, mach)
+        return (
+            (self.center_of_mass.get_value_opt(time) - cp_position)
+            / (2 * self.radius)
+            * self._csys
+        )
+
     def evaluate_static_margin(self):
         """Calculates the static margin of the rocket as a function of time.
 

rocketpy/simulation/flight.py

@@ -3353,8 +3353,8 @@ def static_margin(self):
     def stability_margin(self):
         """Stability margin of the rocket along the flight, it considers the
         variation of the center of pressure position according to the mach
-        number, as well as the variation of the center of gravity position
-        according to the propellant mass evolution.
+        number and angle of attack, as well as the variation of the center of
+        gravity position according to the propellant mass evolution.
 
         Parameters
         ----------
@@ -3366,8 +3366,27 @@ def stability_margin(self):
             Stability margin as a rocketpy.Function of time. The stability margin
             is defined as the distance between the center of pressure and the
             center of gravity, divided by the rocket diameter.
+
+        Notes
+        -----
+        The center of pressure position is computed using the actual lift
+        coefficient CN(α) instead of the lift coefficient derivative CNα. This
+        provides a more accurate stability margin that varies with angle of
+        attack, similar to OpenRocket's implementation.
         """
-        return [(t, self.rocket.stability_margin(m, t)) for t, m in self.mach_number]
+        aoa_values = self.angle_of_attack.y_array
+        mach_values = self.mach_number.y_array
+        time_values = self.time
+
+        results = []
+        for i, t in enumerate(time_values):
+            # Convert angle of attack from degrees to radians
+            alpha_rad = np.deg2rad(aoa_values[i])
+            mach = mach_values[i]
+            sm = self.rocket.get_stability_margin_from_alpha(alpha_rad, mach, t)
+            results.append((t, sm))
+
+        return results
 
     # Rail Button Forces
 

tests/unit/simulation/test_flight.py

@@ -166,6 +166,67 @@ def test_out_of_rail_stability_margin(flight_calisto_custom_wind):
     assert np.isclose(res, 2.14, atol=0.1)
 
 
+def test_stability_margin_uses_angle_of_attack(flight_calisto_custom_wind):
+    """Test that the stability margin calculation accounts for angle of attack.
+
+    The stability margin should use the actual lift coefficient CN(α) instead of
+    just the lift coefficient derivative CNα. This provides a more accurate
+    stability margin that varies with angle of attack.
+
+    Parameters
+    ----------
+    flight_calisto_custom_wind : rocketpy.Flight
+    """
+    # Get stability margin at various times
+    stability_array = flight_calisto_custom_wind.stability_margin[:, 1]
+    time_array = flight_calisto_custom_wind.stability_margin[:, 0]
+
+    # Verify stability margin is a reasonable Function
+    assert len(stability_array) == len(time_array)
+    assert len(stability_array) > 0
+
+    # Verify the rocket's get_stability_margin_from_alpha method works
+    rocket = flight_calisto_custom_wind.rocket
+    alpha = np.deg2rad(5)  # 5 degrees angle of attack
+    mach = 0.5
+    time = 0.0
+
+    sm_from_alpha = rocket.get_stability_margin_from_alpha(alpha, mach, time)
+    sm_from_mach = rocket.stability_margin(mach, time)
+
+    # Both should return reasonable stability margins (positive for stable rocket)
+    assert isinstance(sm_from_alpha, float)
+    assert isinstance(sm_from_mach, float)
+    assert sm_from_alpha > 0  # Should be stable
+    assert sm_from_mach > 0  # Should be stable
+
+
+def test_cp_position_from_alpha_edge_cases(flight_calisto_custom_wind):
+    """Test edge cases for the get_cp_position_from_alpha method.
+
+    Parameters
+    ----------
+    flight_calisto_custom_wind : rocketpy.Flight
+    """
+    rocket = flight_calisto_custom_wind.rocket
+    mach = 0.5
+
+    # Test with zero angle of attack - should fall back to standard cp_position
+    cp_zero_alpha = rocket.get_cp_position_from_alpha(0.0, mach)
+    cp_standard = rocket.cp_position.get_value_opt(mach)
+    assert np.isclose(cp_zero_alpha, cp_standard)
+
+    # Test with small positive angle of attack
+    alpha_small = np.deg2rad(1)
+    cp_small = rocket.get_cp_position_from_alpha(alpha_small, mach)
+    assert isinstance(cp_small, float)
+
+    # Test with larger angle of attack
+    alpha_large = np.deg2rad(10)
+    cp_large = rocket.get_cp_position_from_alpha(alpha_large, mach)
+    assert isinstance(cp_large, float)
+
+
 def test_export_sensor_data(flight_calisto_with_sensors):
     """Test the export of sensor data.