From 673f4d9d3b7863d8af5e81080402ab8d83749688 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 13:45:05 +0000 Subject: [PATCH 1/2] Implement angle-of-attack dependent stability margin MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add new methods to compute center of pressure and stability margin as a function of angle of attack (α) instead of just using the lift coefficient derivative. This provides a more accurate stability margin that varies with angle of attack, similar to OpenRocket. - Add Rocket.get_cp_position_from_alpha(alpha, mach) method - Add Rocket.get_stability_margin_from_alpha(alpha, mach, time) method - Update Flight.stability_margin to use angle of attack from simulation - Add tests for new functionality Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- rocketpy/rocket/rocket.py | 82 ++++++++++++++++++++++++++++ rocketpy/simulation/flight.py | 24 +++++++- tests/unit/simulation/test_flight.py | 61 +++++++++++++++++++++ 3 files changed, 164 insertions(+), 3 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 9ac1bb86e..6764a4689 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -671,6 +671,88 @@ 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. + """ + # Handle edge case where alpha is effectively zero + if abs(alpha) < 1e-10: + 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) < 1e-10: + 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. diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index a7e0374b2..78af33f03 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -3658,8 +3658,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 ---------- @@ -3671,8 +3671,26 @@ 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): + 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 diff --git a/tests/unit/simulation/test_flight.py b/tests/unit/simulation/test_flight.py index a54c5afdb..35329ab53 100644 --- a/tests/unit/simulation/test_flight.py +++ b/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. From 581ad7a775cbbacf758472cd09af81b020dd3b23 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 13:50:41 +0000 Subject: [PATCH 2/2] Address code review feedback - Add comment explaining degrees to radians conversion in Flight.stability_margin - Replace magic number 1e-10 with named constant epsilon in get_cp_position_from_alpha Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- rocketpy/rocket/rocket.py | 7 +++++-- rocketpy/simulation/flight.py | 1 + 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 6764a4689..574b42e9e 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -697,8 +697,11 @@ def get_cp_position_from_alpha(self, alpha, mach): 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) < 1e-10: + if abs(alpha) < epsilon: return self.cp_position.get_value_opt(mach) total_cn = 0.0 @@ -715,7 +718,7 @@ def get_cp_position_from_alpha(self, alpha, mach): total_cn += cn weighted_cp += cn * surface_cp - if abs(total_cn) < 1e-10: + if abs(total_cn) < epsilon: return self.cp_position.get_value_opt(mach) return weighted_cp / total_cn diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 78af33f03..e1370e445 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -3685,6 +3685,7 @@ def stability_margin(self): 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)