Merge pull request #377 from RocketPy-Team/enh/new-stability-margin

ENH: Adding Stability Margin with Mach dependency

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":

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)

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

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()

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