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 5b012300f..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)
@@ -714,11 +706,32 @@ 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.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 9764eb467..32332a908 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,24 @@ 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.plot2D(
+            lower=0,
+            upper=[2, self.rocket.motor.burn_out_time],  # Mach 2 and burnout
+            samples=[20, 20],
+            disp_type="surface",
+            alpha=1,
+        )
+
+        return None
+
     def power_on_drag(self):
         """Plots power on drag of the rocket as a function of time.
 
@@ -113,14 +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()
+        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 9614416ad..5b1f9b914 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(
@@ -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
 
@@ -362,6 +369,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 +403,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 +459,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 f8d098c3d..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
@@ -135,34 +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)
-            + " 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
 
@@ -186,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 0ebee9ae1..f466f3762 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,9 @@ 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)
+        self.evaluate_center_of_pressure()
+        self.evaluate_stability_margin()
         self.evaluate_static_margin()
 
         # Initialize plots and prints object
@@ -303,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
@@ -403,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
@@ -419,47 +445,76 @@ 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_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
 
-        # Calculate static margin
-        self.static_margin = (self.center_of_mass - self.cp_position) / (
-            2 * self.radius
+        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
         )
-        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 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.set_source(
+            lambda time: (self.center_of_mass(time) - self.cp_position(0))
+            / (2 * self.radius)
+            * self._csys
         )
-        return None
+        return self.static_margin
 
     def evaluate_dry_inertias(self):
         """Calculates and returns the rocket's dry inertias relative to
@@ -658,6 +713,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 +753,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 fce90c519..698bb2916 100644
--- a/rocketpy/simulation/flight.py
+++ b/rocketpy/simulation/flight.py
@@ -417,6 +417,8 @@ 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.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
         time. Can be called or accessed as array.
@@ -2332,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):
@@ -2564,6 +2589,26 @@ def static_margin(self):
         """Static margin of the rocket."""
         return self.rocket.static_margin
 
+    @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
+        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.
+        """
+        return [(t, self.rocket.stability_margin(m, t)) for t, m in self.mach_number]
+
     # Rail Button Forces
     @cached_property
     def effective_1rl(self):
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(
diff --git a/tests/test_rocket.py b/tests/test_rocket.py
index 0eebdf1f7..185ce12e6 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
@@ -169,8 +170,8 @@ 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.cp_position, 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)
 
 
 @pytest.mark.parametrize(
@@ -188,8 +189,8 @@ 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 calisto.cp_position == pytest.approx(cpz, 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)
     assert pytest.approx(dimensionless_calisto.static_margin(0), 1e-8) == pytest.approx(
@@ -199,11 +200,11 @@ 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)
-    assert pytest.approx(dimensionless_calisto.cp_position / m, 1e-8) == pytest.approx(
-        calisto.cp_position, 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)
 
 
 def test_add_tail_assert_cp_cm_plus_tail(calisto, dimensionless_calisto, m):
@@ -224,8 +225,10 @@ 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 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(
         top_radius=0.0635 * m,
@@ -240,11 +243,11 @@ 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)
-    assert pytest.approx(dimensionless_calisto.cp_position / m, 1e-8) == pytest.approx(
-        calisto.cp_position, 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)
 
 
 @pytest.mark.parametrize(
@@ -280,7 +283,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 +323,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)
 
@@ -355,8 +358,10 @@ 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 calisto.cp_position == pytest.approx(cpz, 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(
         4,
@@ -372,11 +377,11 @@ 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)
-    assert pytest.approx(dimensionless_calisto.cp_position / m, 1e-8) == pytest.approx(
-        calisto.cp_position, 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)
 
 
 def test_add_cm_eccentricity_assert_properties_set(calisto):