ENH: Adding Stability Margin with Mach dependency

rocketpy/Flight.py

@@ -452,8 +452,8 @@ class Flight:
             Expressed as the absolute vale of the magnitude as a function
             of frequency in Hz. Can be called or accessed as array.
 
-        Flight.staticMargin : Function
-            Rocket's static margin during flight in calibers.
+        Flight.stabilityMargin : rocketpy.Function
+            Rocket's stability margin during flight, in calibers.
 
         Fluid Mechanics:
         Flight.streamVelocityX : Function
@@ -2232,10 +2232,26 @@ def attitudeFrequencyResponse(self):
             samplingFrequency=100,
         )
 
-    @cached_property
-    def staticMargin(self):
-        """Static margin of the rocket."""
-        return self.rocket.staticMargin
+    @funcify_method("Time (s)", "Stability Margin (c)", "spline", "constant")
+    def stabilityMargin(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.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        stability : rocketpy.Function
+            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.
+        """
+        s = [(t, self.rocket.stabilityMargin(t, m)) for t, m in self.MachNumber]
+        return Function(s, "Time (s)", "Stability Margin (c)")
 
     # Rail Button Forces
     @funcify_method("Time (s)", "Upper Rail Button Normal Force (N)", "spline", "zero")

rocketpy/Rocket.py

@@ -102,15 +102,14 @@ class Rocket:
         Rocket.aerodynamicSurfaces : Components
             Collection of aerodynamic surfaces of the rocket. Holds Nose cones,
             Fin sets, and Tails.
-        Rocket.cpPosition : float
-            Rocket's center of pressure position relative to the user defined rocket
-            reference system. See `Rocket.centerOfDryMassPosition` for more information
-            regarding the reference system.
-            Expressed in meters.
-        Rocket.staticMargin : float
-            Float value corresponding to rocket static margin when
-            loaded with propellant in units of rocket diameter or
-            calibers.
+        Rocket.centerOfPressure : rocketpy.Function
+            Rocket's center of pressure position relative to the user defined
+            rocket reference system. See `Rocket.coordinateSystemOrientation`
+            for more information regarding the rocket's coordinate system.
+            Expressed in meters as a function of Mach number.
+        Rocket.totalLiftCoeffDer : rocketpy.Function
+            Rocket's total lift coefficient derivative as a function of Mach
+            number.
         Rocket.powerOffDrag : Function
             Rocket's drag coefficient as a function of Mach number when the
             motor is off.
@@ -225,11 +224,6 @@ def __init__(
         # Rail buttons data initialization
         self.rail_buttons = Components()
 
-        self.cpPosition = 0
-        self.staticMargin = Function(
-            lambda x: 0, inputs="Time (s)", outputs="Static Margin (c)"
-        )
-
         # Define aerodynamic drag coefficients
         self.powerOffDrag = Function(
             powerOffDrag,
@@ -245,7 +239,9 @@ def __init__(
             "linear",
             "constant",
         )
-        self.cpPosition = 0  # Set by self.evaluateStaticMargin()
+        self.centerOfPressure = Function(
+            lambda mach: 0
+        )  # This will be set by the self.evaluateCenterOfPressure()
 
         # Create a, possibly, temporary empty motor
         # self.motors = Components()  # currently unused since only one motor is supported
@@ -262,8 +258,8 @@ def __init__(
         self.evaluateReducedMass()
         self.evaluateThrustToWeight()
 
-        # Evaluate static margin (even though no aerodynamic surfaces are present yet)
-        self.evaluateStaticMargin()
+        # Evaluate center of pressure (even if no aerodynamic surfaces is present yet)
+        self.evaluateCenterOfPressure()
 
         # Initialize plots and prints object
         self.prints = _RocketPrints(self)
@@ -388,48 +384,83 @@ def evaluateThrustToWeight(self):
         self.thrustToWeight.setInputs("Time (s)")
         self.thrustToWeight.setOutputs("Thrust/Weight")
 
-    def evaluateStaticMargin(self):
-        """Calculates and returns the rocket's static margin when
-        loaded with propellant. The static margin is saved and returned
-        in units of rocket diameter or calibers. This function also calculates
-        the rocket center of pressure and total lift coefficients.
+    def evaluateCenterOfPressure(self):
+        """Evaluates rocket center of pressure position relative to user defined
+        rocket reference system. It can be called as many times as needed, as it
+        will update the center of pressure function every time it is called. The
+        code will iterate through all aerodynamic surfaces and consider each of
+        their center of pressure position and derivative of the coefficient of
+        lift as a function of Mach number.
 
         Parameters
         ----------
         None
 
         Returns
         -------
-        self.staticMargin : float
-            Float value corresponding to rocket static margin when
-            loaded with propellant in units of rocket diameter or
-            calibers.
+        None
         """
         # Initialize total lift coefficient derivative and center of pressure position
-        self.totalLiftCoeffDer = 0
-        self.cpPosition = 0
+        self.totalLiftCoeffDer = Function(lambda mach: 0)
+        self.centerOfPressure = Function(lambda mach: 0)
 
         # Calculate total lift coefficient derivative and center of pressure
         if len(self.aerodynamicSurfaces) > 0:
             for surface, position in self.aerodynamicSurfaces:
-                self.totalLiftCoeffDer += surface.clalpha(0)
-                self.cpPosition += surface.clalpha(0) * (
+                self.totalLiftCoeffDer += surface.clalpha
+                self.centerOfPressure += surface.clalpha * (
                     position - self._csys * surface.cpz
                 )
-            self.cpPosition /= self.totalLiftCoeffDer
-
-        # Calculate static margin
-        self.staticMargin = (self.centerOfMass - self.cpPosition) / (2 * self.radius)
-        self.staticMargin *= (
-            self._csys
-        )  # Change sign if coordinate system is upside down
-        self.staticMargin.setInputs("Time (s)")
-        self.staticMargin.setOutputs("Static Margin (c)")
-        self.staticMargin.setDiscrete(
-            lower=0, upper=self.motor.burnOutTime, samples=200
-        )
+            self.centerOfPressure /= self.totalLiftCoeffDer
+
+        self.centerOfPressure.setInputs("Mach Number")
+        self.centerOfPressure.setOutputs("Center of Pressure Position (m)")
+        self.totalLiftCoeffDer.setInputs("Mach Number")
+        self.totalLiftCoeffDer.setOutputs("Total Lift Coefficient Derivative")
+
         return None
 
+    def stabilityMargin(self, time, mach):
+        """Calculates the stability margin of the rocket, defined as the
+        distance between the center of pressure and the center of mass, divided
+        by the rocket's diameter. The mach number needs to be specified to
+        calculate the center of pressure position. The time is used to calculate
+        the center of mass position.
+
+        Parameters
+        ----------
+        time : int, float
+            Time at which the center of mass will be calculated.
+        mach : int, float
+            Mach number at which the center of pressure will be calculated.
+
+        Returns
+        -------
+        stabilityMargin : float
+            Stability margin of the rocket, in calibers.
+        """
+        cp = self.centerOfPressure(mach)
+        cm = self.centerOfMass(time)
+        return ((cm - cp) / (2 * self.radius)) * self._csys
+
+    def initialStabilityMargin(self, mach=0):
+        """Calculates the initial stability margin of the rocket, when the
+        rocket is fully loaded with propellant and the motor is off. The mach
+        number is assumed to be zero by default, but can be changed by the user.
+
+        Parameters
+        ----------
+        mach : int, optional
+            Mach number at which the center of pressure will be calculated,
+            by default 0, representing the static margin.
+
+        Returns
+        -------
+        stabilityMargin : float
+            Initial stability margin of the rocket, in calibers.
+        """
+        return self.stabilityMargin(time=0, mach=mach)
+
     def addMotor(self, motor, position):
         """Adds a motor to the rocket.
 
@@ -461,7 +492,7 @@ def addMotor(self, motor, position):
         self.evaluateCenterOfMass()
         self.evaluateReducedMass()
         self.evaluateThrustToWeight()
-        self.evaluateStaticMargin()
+        self.evaluateCenterOfPressure()
         return None
 
     def addSurfaces(self, surfaces, positions):
@@ -498,7 +529,7 @@ def addSurfaces(self, surfaces, positions):
         except TypeError:
             self.aerodynamicSurfaces.add(surfaces, positions)
 
-        self.evaluateStaticMargin()
+        self.evaluateCenterOfPressure()
         return None
 
     def addTail(

rocketpy/plots/flight_plots.py

@@ -721,8 +721,8 @@ def stability_and_control_data(self):
         fig9 = plt.figure(figsize=(9, 6))
 
         ax1 = plt.subplot(211)
-        ax1.plot(self.flight.staticMargin[:, 0], self.flight.staticMargin[:, 1])
-        ax1.set_xlim(0, self.flight.staticMargin[:, 0][-1])
+        ax1.plot(self.flight.stabilityMargin[:, 0], self.flight.stabilityMargin[:, 1])
+        ax1.set_xlim(0, self.flight.stabilityMargin[:, 0][-1])
         ax1.set_title("Static Margin")
         ax1.set_xlabel("Time (s)")
         ax1.set_ylabel("Static Margin (c)")

rocketpy/plots/rocket_plots.py

@@ -66,7 +66,7 @@ def reducedMass(self):
 
         return None
 
-    def staticMargin(self):
+    def stabilityMargin(self):
         """Plots static margin of the rocket as a function of time.
 
         Parameters
@@ -78,7 +78,18 @@ def staticMargin(self):
         None
         """
 
-        self.rocket.staticMargin()
+        # TODO: it would be interesting to make a 3D plot of stability margin
+        # (https://matplotlib.org/stable/gallery/mplot3d/surface3d.html)
+
+        x = np.linspace(0, self.rocket.motor.burnOutTime, 20)
+        y = np.array([self.rocket.stabilityMargin(t, 0) for t in x])
+
+        plt.plot(x, y)
+        plt.xlabel("Time (s)")
+        plt.ylabel("Stability Margin (calibers)")
+        plt.title("Stability Margin (mach = 0)")
+        plt.grid()
+        plt.show()
 
         return None
 
@@ -148,7 +159,7 @@ def all(self):
         self.totalMass()
         self.reducedMass()
         print("\nAerodynamics Plots")
-        self.staticMargin()
+        self.stabilityMargin()
         self.powerOnDrag()
         self.powerOffDrag()
         self.thrustToWeight()

rocketpy/prints/flight_prints.py

@@ -170,7 +170,7 @@ def out_of_rail_conditions(self):
         )
         print(
             "Rail Departure Static Margin: {:.3f} c".format(
-                self.flight.staticMargin(self.flight.outOfRailTime)
+                self.flight.stabilityMargin(self.flight.outOfRailTime)
             )
         )
         print(

rocketpy/prints/rocket_prints.py

@@ -118,18 +118,20 @@ def rocket_aerodynamics_quantities(self):
             cpz = surface.cp[2]
             print(name + " Center of Pressure to CM: {:.3f}".format(cpz) + " m")
         print(
-            "Distance - Center of Pressure to CM: "
-            + "{:.3f}".format(self.rocket.cpPosition)
+            "Distance - Center of Pressure to CM (Mach=0): "
+            + "{:.3f}".format(self.rocket.centerOfPressure(0))
             + " m"
         )
         print(
-            "Initial Static Margin: "
-            + "{:.3f}".format(self.rocket.staticMargin(0))
+            "Initial Stability Margin (Mach=0): "
+            + "{:.3f}".format(self.rocket.initialStabilityMargin(mach=0))
             + " c"
         )
         print(
-            "Final Static Margin: "
-            + "{:.3f}".format(self.rocket.staticMargin(self.rocket.motor.burnOutTime))
+            "Final Stability Margin (Mach=0): "
+            + "{:.3f}".format(
+                self.rocket.stabilityMargin(time=self.rocket.motor.burnOutTime, mach=0)
+            )
             + " c"
         )
 

tests/test_flight.py

@@ -48,7 +48,7 @@ def compute_static_margin_error_given_distance(position, static_margin, rocket):
             length=0.060,
             position=-1.194656,
         )
-        return rocket.staticMargin(0) - static_margin
+        return rocket.stabilityMargin(0, 0) - static_margin
 
     sol = optimize.root_scalar(
         compute_static_margin_error_given_distance,