Merge branch 'enh/new-stability-margin' of https://github.com/RocketPy-Team/RocketPy into enh/new-stability-margin

Author: MateusStano

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)
@@ -718,7 +710,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)

rocketpy/plots/rocket_plots.py

@@ -73,7 +73,13 @@ 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",
+            alpha=1,
+        )
 
         return None
 
@@ -125,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

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,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
 
@@ -188,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

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

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
@@ -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):