From 4c98f90bbe2815ac3830271d31a5b2a2b7a7232b Mon Sep 17 00:00:00 2001
From: Guilherme <guilherme_fernandes@usp.br>
Date: Mon, 12 Jun 2023 00:41:00 +0200
Subject: [PATCH 01/16] ENH: stabilityMargin as a function of mach and time

---
 rocketpy/Rocket.py | 119 ++++++++++++++++++++++++++++-----------------
 1 file changed, 75 insertions(+), 44 deletions(-)

diff --git a/rocketpy/Rocket.py b/rocketpy/Rocket.py
index 76a8a82a1..5fdb2da6d 100644
--- a/rocketpy/Rocket.py
+++ b/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,11 +384,13 @@ 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
         ----------
@@ -400,36 +398,69 @@ def evaluateStaticMargin(self):
 
         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(

From f29f405c6e0d99f72477579da658e25cd1e2fea2 Mon Sep 17 00:00:00 2001
From: Guilherme <guilherme_fernandes@usp.br>
Date: Mon, 12 Jun 2023 00:41:27 +0200
Subject: [PATCH 02/16] ENH: adjust Flight class to new stability

---
 rocketpy/Flight.py | 28 ++++++++++++++++++++++------
 1 file changed, 22 insertions(+), 6 deletions(-)

diff --git a/rocketpy/Flight.py b/rocketpy/Flight.py
index 78ecbbf0d..aa472f656 100644
--- a/rocketpy/Flight.py
+++ b/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")

From 2f32686e5386f0ab628af4e82b5cef0d8bed3e73 Mon Sep 17 00:00:00 2001
From: Guilherme <guilherme_fernandes@usp.br>
Date: Mon, 12 Jun 2023 01:04:34 +0200
Subject: [PATCH 03/16] ENH: update plots and prints

---
 rocketpy/plots/flight_plots.py   |  4 ++--
 rocketpy/plots/rocket_plots.py   | 17 ++++++++++++++---
 rocketpy/prints/flight_prints.py |  2 +-
 rocketpy/prints/rocket_prints.py | 14 ++++++++------
 4 files changed, 25 insertions(+), 12 deletions(-)

diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py
index b2f6f5c0e..3469bd6e6 100644
--- a/rocketpy/plots/flight_plots.py
+++ b/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)")
diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py
index 1826be771..2fa33dc0c 100644
--- a/rocketpy/plots/rocket_plots.py
+++ b/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()
diff --git a/rocketpy/prints/flight_prints.py b/rocketpy/prints/flight_prints.py
index e8aec1cd3..3084f19a5 100644
--- a/rocketpy/prints/flight_prints.py
+++ b/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(
diff --git a/rocketpy/prints/rocket_prints.py b/rocketpy/prints/rocket_prints.py
index 2ba50f87b..cc14cf6aa 100644
--- a/rocketpy/prints/rocket_prints.py
+++ b/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"
         )
 

From 6fecf20ebeed969a366184aaf5d17adfb22a75cc Mon Sep 17 00:00:00 2001
From: Guilherme <guilherme_fernandes@usp.br>
Date: Mon, 12 Jun 2023 01:04:59 +0200
Subject: [PATCH 04/16] TST: adjusting tests to the new stability margin

---
 tests/test_flight.py |   2 +-
 tests/test_rocket.py | 118 +++++++++++++++++++++++--------------------
 2 files changed, 65 insertions(+), 55 deletions(-)

diff --git a/tests/test_flight.py b/tests/test_flight.py
index ac188492d..9d740e40e 100644
--- a/tests/test_flight.py
+++ b/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,
diff --git a/tests/test_rocket.py b/tests/test_rocket.py
index 8abbd5df8..16a7f4784 100644
--- a/tests/test_rocket.py
+++ b/tests/test_rocket.py
@@ -81,7 +81,7 @@ def mainTrigger(p, h, y):
         noise=(0, 8.3, 0.5),
     )
 
-    static_margin = test_rocket.staticMargin(0)
+    static_margin = test_rocket.stabilityMargin(time=0, mach=0)
 
     # Check if allInfo and static_method methods are working properly
     assert test_rocket.allInfo() == None or not abs(static_margin - 2.05) < 0.01
@@ -151,7 +151,7 @@ def test_coordinate_system_orientation(mock_show):
         4, span=0.100, rootChord=0.120, tipChord=0.040, position=-1.04956
     )
 
-    static_margin_tail_to_nose = rocket_tail_to_nose.staticMargin(0)
+    static_margin_tail_to_nose = rocket_tail_to_nose.stabilityMargin(0, 0)
 
     rocket_nose_to_tail = Rocket(
         radius=127 / 2000,
@@ -173,7 +173,7 @@ def test_coordinate_system_orientation(mock_show):
         4, span=0.100, rootChord=0.120, tipChord=0.040, position=1.04956
     )
 
-    static_margin_nose_to_tail = rocket_nose_to_tail.staticMargin(0)
+    static_margin_nose_to_tail = rocket_nose_to_tail.stabilityMargin(0, 0)
 
     assert (
         rocket_tail_to_nose.allInfo() == None
@@ -256,7 +256,7 @@ def mainTrigger(p, h, y):
         noise=(0, 8.3, 0.5),
     )
 
-    static_margin = test_rocket.staticMargin(0)
+    static_margin = test_rocket.stabilityMargin(time=0, mach=0)
 
     assert test_rocket.allInfo() == None or not abs(static_margin - 2.30) < 0.01
     assert FinSet.draw() == None
@@ -349,23 +349,23 @@ def mainTrigger(p, h, y):
         noise=(0, 8.3, 0.5),
     )
 
-    static_margin = test_rocket.staticMargin(0)
+    static_margin = test_rocket.stabilityMargin(time=0, mach=0)
 
     assert test_rocket.allInfo() == None or not abs(static_margin - 2.03) < 0.01
 
 
 def test_evaluate_static_margin_assert_cp_equals_cm(kg, m, dimensionless_rocket):
     rocket = dimensionless_rocket
-    rocket.evaluateStaticMargin()
+    rocket.evaluateCenterOfPressure()
 
     burnOutTime = rocket.motor.burnOutTime
 
-    assert rocket.centerOfMass(0) / (2 * rocket.radius) == rocket.staticMargin(0)
+    assert rocket.centerOfMass(0) / (2 * rocket.radius) == rocket.stabilityMargin(0, 0)
     assert pytest.approx(
         rocket.centerOfMass(burnOutTime) / (2 * rocket.radius), 1e-12
-    ) == pytest.approx(rocket.staticMargin(burnOutTime), 1e-12)
-    assert rocket.totalLiftCoeffDer == 0
-    assert rocket.cpPosition == 0
+    ) == pytest.approx(rocket.stabilityMargin(time=burnOutTime, mach=0), 1e-12)
+    assert rocket.totalLiftCoeffDer(0) == 0
+    assert rocket.centerOfPressure(0) == 0
 
 
 @pytest.mark.parametrize(
@@ -384,27 +384,29 @@ def test_add_nose_assert_cp_cm_plus_nose(k, type, rocket, dimensionless_rocket,
     clalpha = 2
 
     static_margin_initial = (rocket.centerOfMass(0) - cpz) / (2 * rocket.radius)
-    assert static_margin_initial == pytest.approx(rocket.staticMargin(0), 1e-12)
+    assert static_margin_initial == pytest.approx(rocket.stabilityMargin(0, 0), 1e-12)
 
     static_margin_final = (rocket.centerOfMass(np.inf) - cpz) / (2 * rocket.radius)
-    assert static_margin_final == pytest.approx(rocket.staticMargin(np.inf), 1e-12)
+    assert static_margin_final == pytest.approx(
+        rocket.stabilityMargin(time=np.inf, mach=0), 1e-12
+    )
 
-    assert clalpha == pytest.approx(rocket.totalLiftCoeffDer, 1e-12)
-    assert rocket.cpPosition == pytest.approx(cpz, 1e-12)
+    assert clalpha == pytest.approx(rocket.totalLiftCoeffDer(0), 1e-12)
+    assert rocket.centerOfPressure(0) == pytest.approx(cpz, 1e-12)
 
     dimensionless_rocket.addNose(length=0.55829 * m, kind=type, position=1.278 * m)
-    assert pytest.approx(dimensionless_rocket.staticMargin(0), 1e-12) == pytest.approx(
-        rocket.staticMargin(0), 1e-12
-    )
     assert pytest.approx(
-        dimensionless_rocket.staticMargin(np.inf), 1e-12
-    ) == pytest.approx(rocket.staticMargin(np.inf), 1e-12)
+        dimensionless_rocket.stabilityMargin(0, 0), 1e-12
+    ) == pytest.approx(rocket.stabilityMargin(0, 0), 1e-12)
     assert pytest.approx(
-        dimensionless_rocket.totalLiftCoeffDer, 1e-12
-    ) == pytest.approx(rocket.totalLiftCoeffDer, 1e-12)
-    assert pytest.approx(dimensionless_rocket.cpPosition / m, 1e-12) == pytest.approx(
-        rocket.cpPosition, 1e-12
-    )
+        dimensionless_rocket.stabilityMargin(time=np.inf, mach=0), 1e-12
+    ) == pytest.approx(rocket.stabilityMargin(time=np.inf, mach=0), 1e-12)
+    assert pytest.approx(
+        dimensionless_rocket.totalLiftCoeffDer(0), 1e-12
+    ) == pytest.approx(rocket.totalLiftCoeffDer(0), 1e-12)
+    assert pytest.approx(
+        dimensionless_rocket.centerOfPressure(0) / m, 1e-12
+    ) == pytest.approx(rocket.centerOfPressure(0), 1e-12)
 
 
 def test_add_tail_assert_cp_cm_plus_tail(rocket, dimensionless_rocket, m):
@@ -421,12 +423,14 @@ def test_add_tail_assert_cp_cm_plus_tail(rocket, dimensionless_rocket, m):
     )
 
     static_margin_initial = (rocket.centerOfMass(0) - cpz) / (2 * rocket.radius)
-    assert static_margin_initial == pytest.approx(rocket.staticMargin(0), 1e-12)
+    assert static_margin_initial == pytest.approx(rocket.stabilityMargin(0, 0), 1e-12)
 
     static_margin_final = (rocket.centerOfMass(np.inf) - cpz) / (2 * rocket.radius)
-    assert static_margin_final == pytest.approx(rocket.staticMargin(np.inf), 1e-12)
-    assert np.abs(clalpha) == pytest.approx(np.abs(rocket.totalLiftCoeffDer), 1e-8)
-    assert rocket.cpPosition == cpz
+    assert static_margin_final == pytest.approx(
+        rocket.stabilityMargin(time=np.inf, mach=0), 1e-12
+    )
+    assert np.abs(clalpha) == pytest.approx(np.abs(rocket.totalLiftCoeffDer(0)), 1e-8)
+    assert rocket.centerOfPressure(0) == cpz
 
     dimensionless_rocket.addTail(
         topRadius=0.0635 * m,
@@ -434,18 +438,18 @@ def test_add_tail_assert_cp_cm_plus_tail(rocket, dimensionless_rocket, m):
         length=0.060 * m,
         position=-1.194656 * m,
     )
-    assert pytest.approx(dimensionless_rocket.staticMargin(0), 1e-12) == pytest.approx(
-        rocket.staticMargin(0), 1e-12
-    )
     assert pytest.approx(
-        dimensionless_rocket.staticMargin(np.inf), 1e-12
-    ) == pytest.approx(rocket.staticMargin(np.inf), 1e-12)
+        dimensionless_rocket.stabilityMargin(0, 0), 1e-12
+    ) == pytest.approx(rocket.stabilityMargin(0, 0), 1e-12)
     assert pytest.approx(
-        dimensionless_rocket.totalLiftCoeffDer, 1e-12
-    ) == pytest.approx(rocket.totalLiftCoeffDer, 1e-12)
-    assert pytest.approx(dimensionless_rocket.cpPosition / m, 1e-12) == pytest.approx(
-        rocket.cpPosition, 1e-12
-    )
+        dimensionless_rocket.stabilityMargin(time=np.inf, mach=0), 1e-12
+    ) == pytest.approx(rocket.stabilityMargin(time=np.inf, mach=0), 1e-12)
+    assert pytest.approx(
+        dimensionless_rocket.totalLiftCoeffDer(0), 1e-12
+    ) == pytest.approx(rocket.totalLiftCoeffDer(0), 1e-12)
+    assert pytest.approx(
+        dimensionless_rocket.centerOfPressure(0) / m, 1e-12
+    ) == pytest.approx(rocket.centerOfPressure(0), 1e-12)
 
 
 @pytest.mark.parametrize(
@@ -477,7 +481,9 @@ def test_add_trapezoidal_fins_sweep_angle(
     assert cl_alpha == pytest.approx(expected_clalpha, 0.01)
 
     # Check rocket's center of pressure (just double checking)
-    assert translate - rocket.cpPosition == pytest.approx(expected_cpz_cm, 0.01)
+    assert translate - rocket.centerOfPressure(0) == pytest.approx(
+        expected_cpz_cm, 0.01
+    )
 
 
 @pytest.mark.parametrize(
@@ -509,7 +515,9 @@ def test_add_trapezoidal_fins_sweep_length(
     assert cl_alpha == pytest.approx(expected_clalpha, 0.01)
 
     # Check rocket's center of pressure (just double checking)
-    assert translate - rocket.cpPosition == pytest.approx(expected_cpz_cm, 0.01)
+    assert translate - rocket.centerOfPressure(0) == pytest.approx(
+        expected_cpz_cm, 0.01
+    )
 
     assert isinstance(rocket.aerodynamicSurfaces[0].component, NoseCone)
 
@@ -539,13 +547,15 @@ def test_add_fins_assert_cp_cm_plus_fins(rocket, dimensionless_rocket, m):
     clalpha *= 1 + rocket.radius / (0.1 + rocket.radius)
 
     static_margin_initial = (rocket.centerOfMass(0) - cpz) / (2 * rocket.radius)
-    assert static_margin_initial == pytest.approx(rocket.staticMargin(0), 1e-12)
+    assert static_margin_initial == pytest.approx(rocket.stabilityMargin(0, 0), 1e-12)
 
     static_margin_final = (rocket.centerOfMass(np.inf) - cpz) / (2 * rocket.radius)
-    assert static_margin_final == pytest.approx(rocket.staticMargin(np.inf), 1e-12)
+    assert static_margin_final == pytest.approx(
+        rocket.stabilityMargin(time=np.inf, mach=0), 1e-12
+    )
 
-    assert np.abs(clalpha) == pytest.approx(np.abs(rocket.totalLiftCoeffDer), 1e-12)
-    assert rocket.cpPosition == pytest.approx(cpz, 1e-12)
+    assert np.abs(clalpha) == pytest.approx(np.abs(rocket.totalLiftCoeffDer(0)), 1e-12)
+    assert rocket.centerOfPressure(0) == pytest.approx(cpz, 1e-12)
 
     dimensionless_rocket.addTrapezoidalFins(
         4,
@@ -554,18 +564,18 @@ def test_add_fins_assert_cp_cm_plus_fins(rocket, dimensionless_rocket, m):
         tipChord=0.040 * m,
         position=-1.04956 * m,
     )
-    assert pytest.approx(dimensionless_rocket.staticMargin(0), 1e-12) == pytest.approx(
-        rocket.staticMargin(0), 1e-12
-    )
     assert pytest.approx(
-        dimensionless_rocket.staticMargin(np.inf), 1e-12
-    ) == pytest.approx(rocket.staticMargin(np.inf), 1e-12)
+        dimensionless_rocket.stabilityMargin(0, 0), 1e-12
+    ) == pytest.approx(rocket.stabilityMargin(0, 0), 1e-12)
     assert pytest.approx(
-        dimensionless_rocket.totalLiftCoeffDer, 1e-12
-    ) == pytest.approx(rocket.totalLiftCoeffDer, 1e-12)
-    assert pytest.approx(dimensionless_rocket.cpPosition / m, 1e-12) == pytest.approx(
-        rocket.cpPosition, 1e-12
-    )
+        dimensionless_rocket.stabilityMargin(time=np.inf, mach=0), 1e-12
+    ) == pytest.approx(rocket.stabilityMargin(time=np.inf, mach=0), 1e-12)
+    assert pytest.approx(
+        dimensionless_rocket.totalLiftCoeffDer(0), 1e-12
+    ) == pytest.approx(rocket.totalLiftCoeffDer(0), 1e-12)
+    assert pytest.approx(
+        dimensionless_rocket.centerOfPressure(0) / m, 1e-12
+    ) == pytest.approx(rocket.centerOfPressure(0), 1e-12)
 
 
 def test_add_cm_eccentricity_assert_properties_set(rocket):

From c2c5b615b2dcef64885d215853c612bd0696985d Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 29 Sep 2023 18:55:34 +0200
Subject: [PATCH 05/16] ENH: improve stability margin implementation and docs

---
 rocketpy/rocket/rocket.py     | 126 +++++++++++++++++++++++++---------
 rocketpy/simulation/flight.py |   7 +-
 2 files changed, 97 insertions(+), 36 deletions(-)

diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index 0ebee9ae1..3ffde63a2 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -81,11 +81,20 @@ class Rocket:
     Rocket.aerodynamic_surfaces : list
         Collection of aerodynamic surfaces of the rocket. Holds Nose cones,
         Fin sets, and Tails.
-    Rocket.cp_position : float
-        Rocket's center of pressure position relative to the user defined rocket
-        reference system. See
-        :doc:`Positions and Coordinate Systems </user/positions>` for more
-        information regarding the reference system. Expressed in meters.
+    Rocket.cp_position : Function
+        Function of Mach number expressing the rocket's center of pressure
+        position relative to user defined rocket reference system.
+        See :doc:`Positions and Coordinate Systems </user/positions>`
+        for more information.
+    Rocket.stability_margin : Function
+        Stability margin of the rocket, in calibers, as a function of mach
+        number and time. Stability margin is defined as the distance between
+        the center of pressure and the center of mass, divided by the
+        rocket's diameter.
+    Rocket.static_margin : Function
+        Static margin of the rocket, in calibers, as a function of time. Static
+        margin is defined as the distance between the center of pressure and the
+        center of mass, divided by the rocket's diameter.
     Rocket.static_margin : float
         Float value corresponding to rocket static margin when
         loaded with propellant in units of rocket diameter or calibers.
@@ -221,9 +230,23 @@ def __init__(
         # Rail buttons data initialization
         self.rail_buttons = Components()
 
-        self.cp_position = 0
+        self.cp_position = Function(
+            lambda mach: 0,
+            inputs="Mach Number",
+            outputs="Center of Pressure Position (m)",
+        )
+        self.total_lift_coeff_der = Function(
+            lambda mach: 0,
+            inputs="Mach Number",
+            outputs="Total Lift Coefficient Derivative",
+        )
         self.static_margin = Function(
-            lambda x: 0, inputs="Time (s)", outputs="Static Margin (c)"
+            lambda time: 0, inputs="Time (s)", outputs="Static Margin (c)"
+        )
+        self.stability_margin = Function(
+            lambda mach, time: 0,
+            inputs=["Mach", "Time (s)"],
+            outputs="Stability Margin (c)",
         )
 
         # Define aerodynamic drag coefficients
@@ -241,7 +264,6 @@ def __init__(
             "linear",
             "constant",
         )
-        self.cp_position = 0  # Set by self.evaluate_static_margin()
 
         # Create a, possibly, temporary empty motor
         # self.motors = Components()  # currently unused since only one motor is supported
@@ -261,7 +283,12 @@ def __init__(
         self.evaluate_reduced_mass()
         self.evaluate_thrust_to_weight()
 
-        # Evaluate static margin (even though no aerodynamic surfaces are present yet)
+        # Evaluate stability (even though no aerodynamic surfaces are present yet)
+        # TODO: should we have a funtion call evaluate_stability that evaluates all stability related functions?
+        # or maybe have these three functions' code inside a function called evaluate_stability?
+        # because there isnt much reason to call them independently
+        self.evaluate_center_of_pressure()
+        self.evaluate_stability_margin()
         self.evaluate_static_margin()
 
         # Initialize plots and prints object
@@ -420,46 +447,77 @@ def evaluate_thrust_to_weight(self):
         self.thrust_to_weight.set_inputs("Time (s)")
         self.thrust_to_weight.set_outputs("Thrust/Weight")
 
-    def evaluate_static_margin(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 evaluate_center_of_pressure(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.
 
         Returns
         -------
-        self.static_margin : float
-            Float value corresponding to rocket static margin when
-            loaded with propellant in units of rocket diameter or
-            calibers.
+        self.cp_position : Function
+            Function of Mach number expressing the rocket's center of pressure
+            position relative to user defined rocket reference system.
+            See :doc:`Positions and Coordinate Systems </user/positions>`
+            for more information.
         """
-        # Initialize total lift coefficient derivative and center of pressure position
-        self.total_lift_coeff_der = 0
-        self.cp_position = 0
+        # Re-Initialize total lift coefficient derivative and center of pressure position
+        self.total_lift_coeff_der.set_source(lambda mach: 0)
+        self.cp_position.set_source(lambda mach: 0)
 
         # Calculate total lift coefficient derivative and center of pressure
         if len(self.aerodynamic_surfaces) > 0:
             for aero_surface, position in self.aerodynamic_surfaces:
-                self.total_lift_coeff_der += aero_surface.clalpha(0)
-                self.cp_position += aero_surface.clalpha(0) * (
+                self.total_lift_coeff_der += aero_surface.clalpha
+                self.cp_position += aero_surface.clalpha * (
                     position - self._csys * aero_surface.cpz
                 )
             self.cp_position /= self.total_lift_coeff_der
 
+        return self.cp_position
+
+    def evaluate_stability_margin(self):
+        """Calculates the stability margin of the rocket as a function of mach
+        number and time.
+
+        Returns
+        -------
+        stability_margin : Function
+            Stability margin of the rocket, in calibers, as a function of mach
+            number and time. Stability margin is defined as the distance between
+            the center of pressure and the center of mass, divided by the
+            rocket's diameter.
+        """
+        self.stability_margin.set_source(
+            lambda mach, time: (
+                (self.center_of_mass(time) - self.cp_position(mach)) / (2 * self.radius)
+            )
+            * self._csys
+        )
+        return self.stability_margin
+
+    def evaluate_static_margin(self):
+        """Calculates the static margin of the rocket as a function of time.
+
+        Returns
+        -------
+        static_margin : Function
+            Static margin of the rocket, in calibers, as a function of time.
+            Static margin is defined as the distance between the center of
+            pressure and the center of mass, divided by the rocket's diameter.
+        """
         # Calculate static margin
-        self.static_margin = (self.center_of_mass - self.cp_position) / (
-            2 * self.radius
+        self.static_margin.set_source(
+            lambda time: (self.center_of_mass(time) - self.cp_position(0))
+            / (2 * self.radius)
+            * self._csys
         )
-        self.static_margin *= (
-            self._csys
-        )  # Change sign if coordinate system is upside down
-        self.static_margin.set_inputs("Time (s)")
-        self.static_margin.set_outputs("Static Margin (c)")
-        self.static_margin.set_title("Static Margin")
         self.static_margin.set_discrete(
             lower=0, upper=self.motor.burn_out_time, samples=200
         )
-        return None
+        return self.static_margin
 
     def evaluate_dry_inertias(self):
         """Calculates and returns the rocket's dry inertias relative to
@@ -658,6 +716,8 @@ def add_motor(self, motor, position):
         self.evaluate_inertias()
         self.evaluate_reduced_mass()
         self.evaluate_thrust_to_weight()
+        self.evaluate_center_of_pressure()
+        self.evaluate_stability_margin()
         self.evaluate_static_margin()
         return None
 
@@ -696,6 +756,8 @@ def add_surfaces(self, surfaces, positions):
         except TypeError:
             self.aerodynamic_surfaces.add(surfaces, positions)
 
+        self.evaluate_center_of_pressure()
+        self.evaluate_stability_margin()
         self.evaluate_static_margin()
         return None
 
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index ef127074e..dd103fb71 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -2566,8 +2566,8 @@ def static_margin(self):
         """Static margin of the rocket."""
         return self.rocket.static_margin
 
-    @funcify_method("Time (s)", "Stability Margin (c)", "spline", "constant")
-    def stabilityMargin(self):
+    @funcify_method("Time (s)", "Stability Margin (c)", "linear", "zero")
+    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
@@ -2584,8 +2584,7 @@ def stabilityMargin(self):
             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)")
+        return [(t, self.rocket.stability_margin(m, t)) for t, m in self.mach_number]
 
     # Rail Button Forces
     @cached_property

From 78e15ab9e5e607adc153827ef6d5d3b2d5055922 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 29 Sep 2023 18:56:02 +0200
Subject: [PATCH 06/16] ENH: adapt prints and plots to new changes

---
 rocketpy/plots/flight_plots.py   |  6 +++---
 rocketpy/plots/rocket_plots.py   | 15 ++++++++++++++-
 rocketpy/prints/rocket_prints.py |  2 +-
 3 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py
index 5b012300f..348b6229f 100644
--- a/rocketpy/plots/flight_plots.py
+++ b/rocketpy/plots/flight_plots.py
@@ -714,9 +714,9 @@ def stability_and_control_data(self):
         fig9 = plt.figure(figsize=(9, 6))
 
         ax1 = plt.subplot(211)
-        ax1.plot(self.flight.static_margin[:, 0], self.flight.static_margin[:, 1])
-        ax1.set_xlim(0, self.flight.static_margin[:, 0][-1])
-        ax1.set_title("Static Margin")
+        ax1.plot(self.flight.stability_margin[:, 0], self.flight.stability_margin[:, 1])
+        ax1.set_xlim(0, self.flight.stability_margin[:, 0][-1])
+        ax1.set_title("Stability Margin")
         ax1.set_xlabel("Time (s)")
         ax1.set_ylabel("Static Margin (c)")
         ax1.grid()
diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py
index 9764eb467..4b84fd60d 100644
--- a/rocketpy/plots/rocket_plots.py
+++ b/rocketpy/plots/rocket_plots.py
@@ -12,7 +12,7 @@ class _RocketPlots:
 
     """
 
-    def __init__(self, rocket) -> None:
+    def __init__(self, rocket):
         """Initializes _RocketPlots class.
 
         Parameters
@@ -65,6 +65,18 @@ def static_margin(self):
 
         return None
 
+    def stability_margin(self):
+        """Plots static margin of the rocket as a function of time.
+
+        Returns
+        -------
+        None
+        """
+
+        self.rocket.stability_margin()
+
+        return None
+
     def power_on_drag(self):
         """Plots power on drag of the rocket as a function of time.
 
@@ -119,6 +131,7 @@ def all(self):
         self.reduced_mass()
         print("\nAerodynamics Plots")
         self.static_margin()
+        self.stability_margin()
         self.power_on_drag()
         self.power_off_drag()
         self.thrust_to_weight()
diff --git a/rocketpy/prints/rocket_prints.py b/rocketpy/prints/rocket_prints.py
index f8d098c3d..4cfe32a1a 100644
--- a/rocketpy/prints/rocket_prints.py
+++ b/rocketpy/prints/rocket_prints.py
@@ -148,7 +148,7 @@ def rocket_aerodynamics_quantities(self):
             )
         print(
             "Distance - Center of Pressure to Center of Dry Mass: "
-            + "{:.3f}".format(self.rocket.center_of_mass(0) - self.rocket.cp_position)
+            + "{:.3f}".format(self.rocket.center_of_mass(0) - self.rocket.cp_position(0))
             + " m"
         )
         print(

From 720352836fbb7e45738dd6f168e011d9d310474e Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 29 Sep 2023 18:56:18 +0200
Subject: [PATCH 07/16] ENH: adapt utilities for new changes

---
 rocketpy/utilities.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py
index 48532d5b6..cb7d75960 100644
--- a/rocketpy/utilities.py
+++ b/rocketpy/utilities.py
@@ -555,6 +555,7 @@ def apogee(mass):
         rocket.evaluate_center_of_mass()
         rocket.evaluate_reduced_mass()
         rocket.evaluate_thrust_to_weight()
+        rocket.evaluate_center_of_pressure()
         rocket.evaluate_static_margin()
         # Then we can run the flight simulation
         test_flight = Flight(
@@ -619,6 +620,7 @@ def liftoff_speed(mass):
         rocket.evaluate_center_of_mass()
         rocket.evaluate_reduced_mass()
         rocket.evaluate_thrust_to_weight()
+        rocket.evaluate_center_of_pressure()
         rocket.evaluate_static_margin()
         # Then we can run the flight simulation
         test_flight = Flight(

From 80177a11c12239a21d6881ad7227fffdbaef347f Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 29 Sep 2023 18:56:39 +0200
Subject: [PATCH 08/16] TST: adapt tests for new changes

---
 tests/test_rocket.py | 31 ++++++++++++++++---------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/tests/test_rocket.py b/tests/test_rocket.py
index 0eebdf1f7..bebfec961 100644
--- a/tests/test_rocket.py
+++ b/tests/test_rocket.py
@@ -159,6 +159,7 @@ def test_airfoil(
 
 def test_evaluate_static_margin_assert_cp_equals_cm(dimensionless_calisto):
     rocket = dimensionless_calisto
+    rocket.evaluate_center_of_pressure()
     rocket.evaluate_static_margin()
 
     burn_time = rocket.motor.burn_time
@@ -170,7 +171,7 @@ def test_evaluate_static_margin_assert_cp_equals_cm(dimensionless_calisto):
         rocket.center_of_mass(burn_time[1]) / (2 * rocket.radius), 1e-8
     ) == pytest.approx(rocket.static_margin(burn_time[1]), 1e-8)
     assert pytest.approx(rocket.total_lift_coeff_der, 1e-8) == pytest.approx(0, 1e-8)
-    assert pytest.approx(rocket.cp_position, 1e-8) == pytest.approx(0, 1e-8)
+    assert pytest.approx(rocket.cp_position(0), 1e-8) == pytest.approx(0, 1e-8)
 
 
 @pytest.mark.parametrize(
@@ -189,7 +190,7 @@ def test_add_nose_assert_cp_cm_plus_nose(k, type, calisto, dimensionless_calisto
     assert static_margin_final == pytest.approx(calisto.static_margin(np.inf), 1e-8)
 
     assert clalpha == pytest.approx(calisto.total_lift_coeff_der, 1e-8)
-    assert calisto.cp_position == pytest.approx(cpz, 1e-8)
+    assert calisto.cp_position(0) == pytest.approx(cpz, 1e-8)
 
     dimensionless_calisto.add_nose(length=0.55829 * m, kind=type, position=(1.160) * m)
     assert pytest.approx(dimensionless_calisto.static_margin(0), 1e-8) == pytest.approx(
@@ -201,9 +202,9 @@ def test_add_nose_assert_cp_cm_plus_nose(k, type, calisto, dimensionless_calisto
     assert pytest.approx(
         dimensionless_calisto.total_lift_coeff_der, 1e-8
     ) == pytest.approx(calisto.total_lift_coeff_der, 1e-8)
-    assert pytest.approx(dimensionless_calisto.cp_position / m, 1e-8) == pytest.approx(
-        calisto.cp_position, 1e-8
-    )
+    assert pytest.approx(
+        dimensionless_calisto.cp_position(0) / m, 1e-8
+    ) == pytest.approx(calisto.cp_position(0), 1e-8)
 
 
 def test_add_tail_assert_cp_cm_plus_tail(calisto, dimensionless_calisto, m):
@@ -225,7 +226,7 @@ def test_add_tail_assert_cp_cm_plus_tail(calisto, dimensionless_calisto, m):
     static_margin_final = (calisto.center_of_mass(np.inf) - cpz) / (2 * calisto.radius)
     assert static_margin_final == pytest.approx(calisto.static_margin(np.inf), 1e-8)
     assert np.abs(clalpha) == pytest.approx(np.abs(calisto.total_lift_coeff_der), 1e-8)
-    assert calisto.cp_position == cpz
+    assert calisto.cp_position(0) == cpz
 
     dimensionless_calisto.add_tail(
         top_radius=0.0635 * m,
@@ -242,9 +243,9 @@ def test_add_tail_assert_cp_cm_plus_tail(calisto, dimensionless_calisto, m):
     assert pytest.approx(
         dimensionless_calisto.total_lift_coeff_der, 1e-8
     ) == pytest.approx(calisto.total_lift_coeff_der, 1e-8)
-    assert pytest.approx(dimensionless_calisto.cp_position / m, 1e-8) == pytest.approx(
-        calisto.cp_position, 1e-8
-    )
+    assert pytest.approx(
+        dimensionless_calisto.cp_position(0) / m, 1e-8
+    ) == pytest.approx(calisto.cp_position(0), 1e-8)
 
 
 @pytest.mark.parametrize(
@@ -280,7 +281,7 @@ def test_add_trapezoidal_fins_sweep_angle(
     assert cl_alpha == pytest.approx(expected_clalpha, 0.01)
 
     # Check rocket's center of pressure (just double checking)
-    assert translate - calisto.cp_position == pytest.approx(expected_cpz_cm, 0.01)
+    assert translate - calisto.cp_position(0) == pytest.approx(expected_cpz_cm, 0.01)
 
 
 @pytest.mark.parametrize(
@@ -320,7 +321,7 @@ def test_add_trapezoidal_fins_sweep_length(
     assert cl_alpha == pytest.approx(expected_clalpha, 0.01)
 
     # Check rocket's center of pressure (just double checking)
-    assert translate - calisto.cp_position == pytest.approx(expected_cpz_cm, 0.01)
+    assert translate - calisto.cp_position(0) == pytest.approx(expected_cpz_cm, 0.01)
 
     assert isinstance(calisto.aerodynamic_surfaces[0].component, NoseCone)
 
@@ -356,7 +357,7 @@ def test_add_fins_assert_cp_cm_plus_fins(calisto, dimensionless_calisto, m):
     assert static_margin_final == pytest.approx(calisto.static_margin(np.inf), 1e-8)
 
     assert np.abs(clalpha) == pytest.approx(np.abs(calisto.total_lift_coeff_der), 1e-8)
-    assert calisto.cp_position == pytest.approx(cpz, 1e-8)
+    assert calisto.cp_position(0) == pytest.approx(cpz, 1e-8)
 
     dimensionless_calisto.add_trapezoidal_fins(
         4,
@@ -374,9 +375,9 @@ def test_add_fins_assert_cp_cm_plus_fins(calisto, dimensionless_calisto, m):
     assert pytest.approx(
         dimensionless_calisto.total_lift_coeff_der, 1e-8
     ) == pytest.approx(calisto.total_lift_coeff_der, 1e-8)
-    assert pytest.approx(dimensionless_calisto.cp_position / m, 1e-8) == pytest.approx(
-        calisto.cp_position, 1e-8
-    )
+    assert pytest.approx(
+        dimensionless_calisto.cp_position(0) / m, 1e-8
+    ) == pytest.approx(calisto.cp_position(0), 1e-8)
 
 
 def test_add_cm_eccentricity_assert_properties_set(calisto):

From c31c4e70753036c087c63b9d380efcf63a68c857 Mon Sep 17 00:00:00 2001
From: Lint Action <lint-action@samuelmeuli.com>
Date: Fri, 29 Sep 2023 22:56:20 +0000
Subject: [PATCH 09/16] Fix code style issues with Black

---
 rocketpy/prints/rocket_prints.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/rocketpy/prints/rocket_prints.py b/rocketpy/prints/rocket_prints.py
index 4cfe32a1a..d8a4532db 100644
--- a/rocketpy/prints/rocket_prints.py
+++ b/rocketpy/prints/rocket_prints.py
@@ -148,7 +148,9 @@ def rocket_aerodynamics_quantities(self):
             )
         print(
             "Distance - Center of Pressure to Center of Dry Mass: "
-            + "{:.3f}".format(self.rocket.center_of_mass(0) - self.rocket.cp_position(0))
+            + "{:.3f}".format(
+                self.rocket.center_of_mass(0) - self.rocket.cp_position(0)
+            )
             + " m"
         )
         print(

From eab0ab0b253b16c75509b5f338ed19aaf694b03e Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Sat, 30 Sep 2023 00:59:23 +0200
Subject: [PATCH 10/16] MAINT: remove unecessary todo

---
 rocketpy/rocket/rocket.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index 3ffde63a2..ea4d3d7c1 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -284,9 +284,6 @@ def __init__(
         self.evaluate_thrust_to_weight()
 
         # Evaluate stability (even though no aerodynamic surfaces are present yet)
-        # TODO: should we have a funtion call evaluate_stability that evaluates all stability related functions?
-        # or maybe have these three functions' code inside a function called evaluate_stability?
-        # because there isnt much reason to call them independently
         self.evaluate_center_of_pressure()
         self.evaluate_stability_margin()
         self.evaluate_static_margin()

From 80c712ce10e93743bed2091da8ad7ce8a8edbbb3 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Sat, 30 Sep 2023 01:11:49 +0200
Subject: [PATCH 11/16] TST: fix total_lift_coeff_der test issues

---
 tests/test_rocket.py | 24 ++++++++++++++----------
 1 file changed, 14 insertions(+), 10 deletions(-)

diff --git a/tests/test_rocket.py b/tests/test_rocket.py
index bebfec961..185ce12e6 100644
--- a/tests/test_rocket.py
+++ b/tests/test_rocket.py
@@ -170,7 +170,7 @@ def test_evaluate_static_margin_assert_cp_equals_cm(dimensionless_calisto):
     assert pytest.approx(
         rocket.center_of_mass(burn_time[1]) / (2 * rocket.radius), 1e-8
     ) == pytest.approx(rocket.static_margin(burn_time[1]), 1e-8)
-    assert pytest.approx(rocket.total_lift_coeff_der, 1e-8) == pytest.approx(0, 1e-8)
+    assert pytest.approx(rocket.total_lift_coeff_der(0), 1e-8) == pytest.approx(0, 1e-8)
     assert pytest.approx(rocket.cp_position(0), 1e-8) == pytest.approx(0, 1e-8)
 
 
@@ -189,7 +189,7 @@ def test_add_nose_assert_cp_cm_plus_nose(k, type, calisto, dimensionless_calisto
     static_margin_final = (calisto.center_of_mass(np.inf) - cpz) / (2 * calisto.radius)
     assert static_margin_final == pytest.approx(calisto.static_margin(np.inf), 1e-8)
 
-    assert clalpha == pytest.approx(calisto.total_lift_coeff_der, 1e-8)
+    assert clalpha == pytest.approx(calisto.total_lift_coeff_der(0), 1e-8)
     assert calisto.cp_position(0) == pytest.approx(cpz, 1e-8)
 
     dimensionless_calisto.add_nose(length=0.55829 * m, kind=type, position=(1.160) * m)
@@ -200,8 +200,8 @@ def test_add_nose_assert_cp_cm_plus_nose(k, type, calisto, dimensionless_calisto
         dimensionless_calisto.static_margin(np.inf), 1e-8
     ) == pytest.approx(calisto.static_margin(np.inf), 1e-8)
     assert pytest.approx(
-        dimensionless_calisto.total_lift_coeff_der, 1e-8
-    ) == pytest.approx(calisto.total_lift_coeff_der, 1e-8)
+        dimensionless_calisto.total_lift_coeff_der(0), 1e-8
+    ) == pytest.approx(calisto.total_lift_coeff_der(0), 1e-8)
     assert pytest.approx(
         dimensionless_calisto.cp_position(0) / m, 1e-8
     ) == pytest.approx(calisto.cp_position(0), 1e-8)
@@ -225,7 +225,9 @@ def test_add_tail_assert_cp_cm_plus_tail(calisto, dimensionless_calisto, m):
 
     static_margin_final = (calisto.center_of_mass(np.inf) - cpz) / (2 * calisto.radius)
     assert static_margin_final == pytest.approx(calisto.static_margin(np.inf), 1e-8)
-    assert np.abs(clalpha) == pytest.approx(np.abs(calisto.total_lift_coeff_der), 1e-8)
+    assert np.abs(clalpha) == pytest.approx(
+        np.abs(calisto.total_lift_coeff_der(0)), 1e-8
+    )
     assert calisto.cp_position(0) == cpz
 
     dimensionless_calisto.add_tail(
@@ -241,8 +243,8 @@ def test_add_tail_assert_cp_cm_plus_tail(calisto, dimensionless_calisto, m):
         dimensionless_calisto.static_margin(np.inf), 1e-8
     ) == pytest.approx(calisto.static_margin(np.inf), 1e-8)
     assert pytest.approx(
-        dimensionless_calisto.total_lift_coeff_der, 1e-8
-    ) == pytest.approx(calisto.total_lift_coeff_der, 1e-8)
+        dimensionless_calisto.total_lift_coeff_der(0), 1e-8
+    ) == pytest.approx(calisto.total_lift_coeff_der(0), 1e-8)
     assert pytest.approx(
         dimensionless_calisto.cp_position(0) / m, 1e-8
     ) == pytest.approx(calisto.cp_position(0), 1e-8)
@@ -356,7 +358,9 @@ def test_add_fins_assert_cp_cm_plus_fins(calisto, dimensionless_calisto, m):
     static_margin_final = (calisto.center_of_mass(np.inf) - cpz) / (2 * calisto.radius)
     assert static_margin_final == pytest.approx(calisto.static_margin(np.inf), 1e-8)
 
-    assert np.abs(clalpha) == pytest.approx(np.abs(calisto.total_lift_coeff_der), 1e-8)
+    assert np.abs(clalpha) == pytest.approx(
+        np.abs(calisto.total_lift_coeff_der(0)), 1e-8
+    )
     assert calisto.cp_position(0) == pytest.approx(cpz, 1e-8)
 
     dimensionless_calisto.add_trapezoidal_fins(
@@ -373,8 +377,8 @@ def test_add_fins_assert_cp_cm_plus_fins(calisto, dimensionless_calisto, m):
         dimensionless_calisto.static_margin(np.inf), 1e-8
     ) == pytest.approx(calisto.static_margin(np.inf), 1e-8)
     assert pytest.approx(
-        dimensionless_calisto.total_lift_coeff_der, 1e-8
-    ) == pytest.approx(calisto.total_lift_coeff_der, 1e-8)
+        dimensionless_calisto.total_lift_coeff_der(0), 1e-8
+    ) == pytest.approx(calisto.total_lift_coeff_der(0), 1e-8)
     assert pytest.approx(
         dimensionless_calisto.cp_position(0) / m, 1e-8
     ) == pytest.approx(calisto.cp_position(0), 1e-8)

From e110fffe77807d8ddb966ca6049e0328d2baaab9 Mon Sep 17 00:00:00 2001
From: Gui-FernandesBR <guilherme_fernandes@usp.br>
Date: Sat, 30 Sep 2023 00:00:48 -0300
Subject: [PATCH 12/16] ENH: improve stability related plots and prints

---
 rocketpy/plots/flight_plots.py   | 23 ++++++++++++++++++++++-
 rocketpy/plots/rocket_plots.py   |  7 ++++++-
 rocketpy/prints/flight_prints.py | 29 +++++++++++++++++++++++++++--
 rocketpy/prints/rocket_prints.py | 24 ++++++++++++++----------
 rocketpy/simulation/flight.py    | 23 +++++++++++++++++++++++
 5 files changed, 92 insertions(+), 14 deletions(-)

diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py
index 348b6229f..4668eb13b 100644
--- a/rocketpy/plots/flight_plots.py
+++ b/rocketpy/plots/flight_plots.py
@@ -718,7 +718,28 @@ def stability_and_control_data(self):
         ax1.set_xlim(0, self.flight.stability_margin[:, 0][-1])
         ax1.set_title("Stability Margin")
         ax1.set_xlabel("Time (s)")
-        ax1.set_ylabel("Static Margin (c)")
+        ax1.set_ylabel("Stability Margin (c)")
+        ax1.set_xlim(0, self.first_event_time)
+        ax1.axvline(
+            x=self.flight.out_of_rail_time,
+            color="r",
+            linestyle="--",
+            label="Out of Rail Time",
+        )
+        ax1.axvline(
+            x=self.flight.rocket.motor.burn_out_time,
+            color="g",
+            linestyle="--",
+            label="Burn Out Time",
+        )
+
+        ax1.axvline(
+            x=self.flight.apogee_time,
+            color="m",
+            linestyle="--",
+            label="Apogee Time",
+        )
+        ax1.legend()
         ax1.grid()
 
         ax2 = plt.subplot(212)
diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py
index 4b84fd60d..9aacdb17e 100644
--- a/rocketpy/plots/rocket_plots.py
+++ b/rocketpy/plots/rocket_plots.py
@@ -73,7 +73,12 @@ def stability_margin(self):
         None
         """
 
-        self.rocket.stability_margin()
+        self.rocket.stability_margin.plot2D(
+            lower=0,
+            upper=[2, self.rocket.motor.burn_out_time],  # Mach 2 and burnout
+            samples=[20, 20],
+            disp_type="surface",
+        )
 
         return None
 
diff --git a/rocketpy/prints/flight_prints.py b/rocketpy/prints/flight_prints.py
index 9614416ad..4be5ca994 100644
--- a/rocketpy/prints/flight_prints.py
+++ b/rocketpy/prints/flight_prints.py
@@ -155,8 +155,8 @@ def out_of_rail_conditions(self):
             )
         )
         print(
-            "Rail Departure Static Margin: {:.3f} c".format(
-                self.flight.static_margin(self.flight.out_of_rail_time)
+            "Rail Departure Stability Margin: {:.3f} c".format(
+                self.flight.stability_margin(self.flight.out_of_rail_time)
             )
         )
         print(
@@ -362,6 +362,11 @@ def maximum_values(self):
                 self.flight.max_acceleration_power_off_time,
             )
         )
+        print(
+            "Maximum Stability Margin: {:.3f} c at {:.2f} s".format(
+                self.flight.max_stability_margin, self.flight.max_stability_margin_time
+            )
+        )
 
         if (
             len(self.flight.rocket.rail_buttons) == 0
@@ -391,6 +396,22 @@ def maximum_values(self):
             )
         return None
 
+    def stability_margin(self):
+        """Prints out the maximum and minimum stability margin available
+        about the flight."""
+        print("\nStability Margin\n")
+        print(
+            "Maximum Stability Margin: {:.3f} c at {:.2f} s".format(
+                self.flight.max_stability_margin, self.flight.max_stability_margin_time
+            )
+        )
+        print(
+            "Minimum Stability Margin: {:.3f} c at {:.2f} s".format(
+                self.flight.min_stability_margin, self.flight.min_stability_margin_time
+            )
+        )
+        return None
+
     def all(self):
         """Prints out all data available about the Flight.
 
@@ -431,6 +452,10 @@ def all(self):
         self.impact_conditions()
         print()
 
+        # Print stability margin
+        self.stability_margin()
+        print()
+
         # Print maximum values
         self.maximum_values()
         print()
diff --git a/rocketpy/prints/rocket_prints.py b/rocketpy/prints/rocket_prints.py
index d8a4532db..4782de867 100644
--- a/rocketpy/prints/rocket_prints.py
+++ b/rocketpy/prints/rocket_prints.py
@@ -135,36 +135,40 @@ def rocket_aerodynamics_quantities(self):
                 + "/rad"
             )
 
-        print("\nAerodynamics Center of Pressure\n")
+        print("\nCenter of Pressure\n")
         for surface, position in self.rocket.aerodynamic_surfaces:
             name = surface.name
-            cpz = surface.cp[2]
+            cpz = surface.cp[2]  # relative to the user defined coordinate system
             print(
                 name
-                + " Center of Pressure to CM: {:.3f}".format(
+                + " Center of Pressure position: {:.3f}".format(
                     position - self.rocket._csys * cpz
                 )
                 + " m"
             )
+        print("\nStability\n")
         print(
-            "Distance - Center of Pressure to Center of Dry Mass: "
-            + "{:.3f}".format(
-                self.rocket.center_of_mass(0) - self.rocket.cp_position(0)
-            )
-            + " m"
+            f"Center of Mass position (time=0): {self.rocket.center_of_mass(0):.3f} m"
         )
         print(
-            "Initial Static Margin: "
+            "Initial Static Margin (mach=0, time=0): "
             + "{:.3f}".format(self.rocket.static_margin(0))
             + " c"
         )
         print(
-            "Final Static Margin: "
+            "Final Static Margin (mach=0, time=burn_out): "
             + "{:.3f}".format(
                 self.rocket.static_margin(self.rocket.motor.burn_out_time)
             )
             + " c"
         )
+        print(
+            "Rocket Center of Mass (time=0) - Center of Pressure (mach=0): "
+            + "{:.3f}".format(
+                abs(self.rocket.center_of_mass(0) - self.rocket.cp_position(0))
+            )
+            + " m\n"
+        )
 
         return None
 
diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index dd103fb71..327b66047 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -2334,6 +2334,29 @@ def max_mach_number(self):
         """Maximum Mach number."""
         return self.mach_number(self.max_mach_number_time)
 
+    # Stability Margin
+    @cached_property
+    def max_stability_margin_time(self):
+        """Time of maximum stability margin."""
+        max_stability_margin_time_index = np.argmax(self.stability_margin[:, 1])
+        return self.stability_margin[max_stability_margin_time_index, 0]
+    
+    @cached_property
+    def max_stability_margin(self):
+        """Maximum stability margin."""
+        return self.stability_margin(self.max_stability_margin_time)
+
+    @cached_property
+    def min_stability_margin_time(self):
+        """Time of minimum stability margin."""
+        min_stability_margin_time_index = np.argmin(self.stability_margin[:, 1])
+        return self.stability_margin[min_stability_margin_time_index, 0]
+    
+    @cached_property
+    def min_stability_margin(self):
+        """Minimum stability margin."""
+        return self.stability_margin(self.min_stability_margin_time)
+
     # Reynolds Number
     @funcify_method("Time (s)", "Reynolds Number", "spline", "zero")
     def reynolds_number(self):

From e70825ca7bd957d3796874b899ec0ae6a0154ab8 Mon Sep 17 00:00:00 2001
From: Lint Action <lint-action@samuelmeuli.com>
Date: Sat, 30 Sep 2023 03:01:52 +0000
Subject: [PATCH 13/16] Fix code style issues with Black

---
 rocketpy/simulation/flight.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index 327b66047..f04f91533 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -2340,7 +2340,7 @@ def max_stability_margin_time(self):
         """Time of maximum stability margin."""
         max_stability_margin_time_index = np.argmax(self.stability_margin[:, 1])
         return self.stability_margin[max_stability_margin_time_index, 0]
-    
+
     @cached_property
     def max_stability_margin(self):
         """Maximum stability margin."""
@@ -2351,7 +2351,7 @@ def min_stability_margin_time(self):
         """Time of minimum stability margin."""
         min_stability_margin_time_index = np.argmin(self.stability_margin[:, 1])
         return self.stability_margin[min_stability_margin_time_index, 0]
-    
+
     @cached_property
     def min_stability_margin(self):
         """Minimum stability margin."""

From 3eeb5a1bc2ff00c7846f97727ad5da72924427bf Mon Sep 17 00:00:00 2001
From: Gui-FernandesBR <guilherme_fernandes@usp.br>
Date: Sat, 30 Sep 2023 00:10:41 -0300
Subject: [PATCH 14/16] ENH: improve plots not related to stability

---
 rocketpy/mathutils/function.py   | 17 +++++++++++++++--
 rocketpy/plots/flight_plots.py   | 12 ++----------
 rocketpy/plots/rocket_plots.py   | 23 ++++++++++++++++++++---
 rocketpy/prints/flight_prints.py |  9 ++++++++-
 rocketpy/prints/rocket_prints.py | 12 ++++--------
 rocketpy/rocket/rocket.py        |  3 +++
 6 files changed, 52 insertions(+), 24 deletions(-)

diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py
index 371fa492d..1beb3f8ec 100644
--- a/rocketpy/mathutils/function.py
+++ b/rocketpy/mathutils/function.py
@@ -1152,6 +1152,8 @@ def plot2D(
         samples=[30, 30],
         force_data=True,
         disp_type="surface",
+        alpha=0.6,
+        cmap="viridis",
     ):
         """Plot 2-Dimensional Function, from a lower limit to an upper limit,
         by sampling the Function several times in the interval. The title of
@@ -1185,6 +1187,12 @@ def plot2D(
         disp_type : string, optional
             Display type of plotted graph, which can be surface, wireframe,
             contour, or contourf. Default value is surface.
+        alpha : float, optional
+            Transparency of plotted graph, which can be a value between 0 and
+            1. Default value is 0.6.
+        cmap : string, optional
+            Colormap of plotted graph, which can be any of the colormaps
+            available in matplotlib. Default value is viridis.
 
         Returns
         -------
@@ -1221,6 +1229,9 @@ def plot2D(
         mesh = [[mesh_x_flat[i], mesh_y_flat[i]] for i in range(len(mesh_x_flat))]
         # Evaluate function at all mesh nodes and convert it to matrix
         z = np.array(self.get_value(mesh)).reshape(mesh_x.shape)
+        z_min, z_max = z.min(), z.max()
+        color_map = plt.cm.get_cmap(cmap)
+        norm = plt.Normalize(z_min, z_max)
         # Plot function
         if disp_type == "surface":
             surf = axes.plot_surface(
@@ -1229,9 +1240,11 @@ def plot2D(
                 z,
                 rstride=1,
                 cstride=1,
-                # cmap=cm.coolwarm,
+                cmap=color_map,
                 linewidth=0,
-                alpha=0.6,
+                alpha=alpha,
+                vmin=z_min,
+                vmax=z_max,
             )
             figure.colorbar(surf)
         elif disp_type == "wireframe":
diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py
index 4668eb13b..62c751a55 100644
--- a/rocketpy/plots/flight_plots.py
+++ b/rocketpy/plots/flight_plots.py
@@ -682,19 +682,11 @@ def fluid_mechanics_data(self):
 
         ax4 = plt.subplot(414)
         ax4.plot(self.flight.angle_of_attack[:, 0], self.flight.angle_of_attack[:, 1])
-        # Make sure bottom and top limits are different
-        if (
-            self.flight.out_of_rail_time
-            * self.flight.angle_of_attack(self.flight.out_of_rail_time)
-            != 0
-        ):
-            ax4.set_xlim(
-                self.flight.out_of_rail_time, 10 * self.flight.out_of_rail_time + 1
-            )
-            ax4.set_ylim(0, self.flight.angle_of_attack(self.flight.out_of_rail_time))
         ax4.set_title("Angle of Attack")
         ax4.set_xlabel("Time (s)")
         ax4.set_ylabel("Angle of Attack (°)")
+        ax4.set_xlim(self.flight.out_of_rail_time, self.first_event_time)
+        ax4.set_ylim(0, self.flight.angle_of_attack(self.flight.out_of_rail_time) + 15)
         ax4.grid()
 
         plt.subplots_adjust(hspace=0.5)
diff --git a/rocketpy/plots/rocket_plots.py b/rocketpy/plots/rocket_plots.py
index 9aacdb17e..32332a908 100644
--- a/rocketpy/plots/rocket_plots.py
+++ b/rocketpy/plots/rocket_plots.py
@@ -78,6 +78,7 @@ def stability_margin(self):
             upper=[2, self.rocket.motor.burn_out_time],  # Mach 2 and burnout
             samples=[20, 20],
             disp_type="surface",
+            alpha=1,
         )
 
         return None
@@ -130,15 +131,31 @@ def all(self):
         None
         """
 
-        # Show plots
+        # Mass Plots
         print("\nMass Plots")
+        print("-" * 40)
         self.total_mass()
         self.reduced_mass()
+
+        # Aerodynamics Plots
         print("\nAerodynamics Plots")
-        self.static_margin()
-        self.stability_margin()
+        print("-" * 40)
+
+        # Drag Plots
+        print("Drag Plots")
+        print("-" * 20)  # Separator for Drag Plots
         self.power_on_drag()
         self.power_off_drag()
+
+        # Stability Plots
+        print("\nStability Plots")
+        print("-" * 20)  # Separator for Stability Plots
+        self.static_margin()
+        self.stability_margin()
+
+        # Thrust-to-Weight Plot
+        print("\nThrust-to-Weight Plot")
+        print("-" * 40)
         self.thrust_to_weight()
 
         return None
diff --git a/rocketpy/prints/flight_prints.py b/rocketpy/prints/flight_prints.py
index 4be5ca994..5b1f9b914 100644
--- a/rocketpy/prints/flight_prints.py
+++ b/rocketpy/prints/flight_prints.py
@@ -301,12 +301,19 @@ def impact_conditions(self):
             print("\nImpact Conditions\n")
             print("X Impact: {:.3f} m".format(self.flight.x_impact))
             print("Y Impact: {:.3f} m".format(self.flight.y_impact))
+            print("Latitude: {:.7f}°".format(self.flight.latitude(self.flight.t_final)))
+            print(
+                "Longitude: {:.7f}°".format(self.flight.longitude(self.flight.t_final))
+            )
             print("Time of Impact: {:.3f} s".format(self.flight.t_final))
             print("Velocity at Impact: {:.3f} m/s".format(self.flight.impact_velocity))
         elif self.flight.terminate_on_apogee is False:
             print("End of Simulation")
-            print("Time: {:.3f} s".format(self.flight.solution[-1][0]))
+            t_final = self.flight.solution[-1][0]
+            print("Time: {:.3f} s".format(t_final))
             print("Altitude: {:.3f} m".format(self.flight.solution[-1][3]))
+            print("Latitude: {:.3f}°".format(self.flight.latitude(t_final)))
+            print("Longitude: {:.3f}°".format(self.flight.longitude(t_final)))
 
         return None
 
diff --git a/rocketpy/prints/rocket_prints.py b/rocketpy/prints/rocket_prints.py
index 4782de867..e09d999f5 100644
--- a/rocketpy/prints/rocket_prints.py
+++ b/rocketpy/prints/rocket_prints.py
@@ -63,7 +63,7 @@ def inertia_details(self):
             )
         )
         print(
-            "Rocket Inertia (with motor, but without propellant) 23: {:.3f} kg*m2".format(
+            "Rocket Inertia (with motor, but without propellant) 23: {:.3f} kg*m2\n".format(
                 self.rocket.dry_I_23
             )
         )
@@ -82,7 +82,7 @@ def rocket_geometrical_parameters(self):
         print("Rocket Frontal Area: " + "{:.6f}".format(self.rocket.area) + " m2")
         print("\nRocket Distances")
         print(
-            "Rocket Center of Dry Mass - Center of Mass withour Motor: "
+            "Rocket Center of Dry Mass - Center of Mass without Motor: "
             + "{:.3f} m".format(
                 abs(
                     self.rocket.center_of_mass_without_motor
@@ -91,7 +91,7 @@ def rocket_geometrical_parameters(self):
             )
         )
         print(
-            "Rocket Center of Dry Mass - Nozzle Exit Distance: "
+            "Rocket Center of Dry Mass - Nozzle Exit: "
             + "{:.3f} m".format(
                 abs(
                     self.rocket.center_of_dry_mass_position - self.rocket.motor_position
@@ -109,7 +109,7 @@ def rocket_geometrical_parameters(self):
         )
         print(
             "Rocket Center of Mass - Rocket Loaded Center of Mass: "
-            + "{:.3f} m".format(
+            + "{:.3f} m\n".format(
                 abs(
                     self.rocket.center_of_mass(0)
                     - self.rocket.center_of_dry_mass_position
@@ -192,18 +192,14 @@ def all(self):
         """
         # Print inertia details
         self.inertia_details()
-        print()
 
         # Print rocket geometrical parameters
         self.rocket_geometrical_parameters()
-        print()
 
         # Print rocket aerodynamics quantities
         self.rocket_aerodynamics_quantities()
-        print()
 
         # Print parachute data
         self.parachute_data()
-        print()
 
         return None
diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index ea4d3d7c1..a5217e318 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -327,6 +327,7 @@ def evaluate_total_mass(self):
         # Calculate total mass by summing up propellant and dry mass
         self.total_mass = self.mass + self.motor.total_mass
         self.total_mass.set_outputs("Total Mass (Rocket + Propellant) (kg)")
+        self.total_mass.set_title("Total Mass (Rocket + Propellant) (kg) x Time (s)")
 
         # Return total mass
         return self.total_mass
@@ -427,6 +428,7 @@ def evaluate_reduced_mass(self):
         # calculate reduced mass
         self.reduced_mass = motor_mass * mass / (motor_mass + mass)
         self.reduced_mass.set_outputs("Reduced Mass (kg)")
+        self.reduced_mass.set_title("Reduced Mass (kg) x Time (s)")
 
         # Return reduced mass
         return self.reduced_mass
@@ -443,6 +445,7 @@ def evaluate_thrust_to_weight(self):
         self.thrust_to_weight = self.motor.thrust / (9.80665 * self.total_mass)
         self.thrust_to_weight.set_inputs("Time (s)")
         self.thrust_to_weight.set_outputs("Thrust/Weight")
+        self.thrust_to_weight.set_title(None)
 
     def evaluate_center_of_pressure(self):
         """Evaluates rocket center of pressure position relative to user defined

From 0036e172b34bd014301a0b92b12cfcc44122ddcb Mon Sep 17 00:00:00 2001
From: MateusStano <69485049+MateusStano@users.noreply.github.com>
Date: Sun, 1 Oct 2023 16:12:46 -0300
Subject: [PATCH 15/16] Update rocketpy/simulation/flight.py

Co-authored-by: Pedro Henrique Marinho Bressan <87212571+phmbressan@users.noreply.github.com>
---
 rocketpy/simulation/flight.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py
index f04f91533..698bb2916 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -417,7 +417,7 @@ class Flight:
         of frequency in Hz. Can be called or accessed as array.
     Flight.static_margin : Function
         Rocket's static margin during flight in calibers.
-    Flight.stabilityMargin : rocketpy.Function
+    Flight.stability_margin : rocketpy.Function
             Rocket's stability margin during flight, in calibers.
     Flight.stream_velocity_x : Function
         Freestream velocity x (East) component, in m/s, as a function of

From 312f137365f3b93a72570c0600f75fba2546bbb7 Mon Sep 17 00:00:00 2001
From: MateusStano <mateusstano@usp.br>
Date: Fri, 6 Oct 2023 20:49:41 +0200
Subject: [PATCH 16/16] ENH: remove unecessary static margin discritization

---
 rocketpy/rocket/rocket.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py
index ea4d3d7c1..06590975c 100644
--- a/rocketpy/rocket/rocket.py
+++ b/rocketpy/rocket/rocket.py
@@ -511,9 +511,6 @@ def evaluate_static_margin(self):
             / (2 * self.radius)
             * self._csys
         )
-        self.static_margin.set_discrete(
-            lower=0, upper=self.motor.burn_out_time, samples=200
-        )
         return self.static_margin
 
     def evaluate_dry_inertias(self):