Merge pull request #595 from RocketPy-Team/enh/pre-calculate-attribut…

…es-in-rocket-class

ENH: Pre-calculate attributes in Rocket class

CHANGELOG.md

@@ -32,6 +32,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 
 ### Added
 
+- ENH: Pre-calculate attributes in Rocket class [#595](https://github.com/RocketPy-Team/RocketPy/pull/595)
 - ENH: Complex step differentiation [#594](https://github.com/RocketPy-Team/RocketPy/pull/594)
 - ENH: Exponential backoff decorator (fix #449) [#588](https://github.com/RocketPy-Team/RocketPy/pull/588)
 - ENH: Function Validation Rework & Swap `np.searchsorted` to `bisect_left` [#582](https://github.com/RocketPy-Team/RocketPy/pull/582)

docs/user/rocket.rst

@@ -12,6 +12,7 @@ Defining a Rocket in RocketPy is simple and requires a few steps:
 4. Add, if desired, parachutes;
 5. Set, if desired, rail guides;
 6. See results.
+7. Inertia Tensors.
 
 Lets go through each of these steps in detail.
 
@@ -428,3 +429,54 @@ The lets check all the information available about the rocket:
 .. jupyter-execute::
 
     calisto.all_info()
+
+7. Inertia Tensors
+------------------
+
+The inertia tensor of the rocket at a given time can be obtained using the
+``get_inertia_tensor_at_time`` method. This method evaluates each component of
+the inertia tensor at the specified time and returns a
+:class:`rocketpy.mathutils.Matrix` object.
+
+The inertia tensor is a matrix that looks like this:
+
+.. math::
+    :label: inertia_tensor
+
+    \mathbf{I} = \begin{bmatrix}
+    I_{11} & I_{12} & I_{13} \\
+    I_{21} & I_{22} & I_{23} \\
+    I_{31} & I_{32} & I_{33}
+    \end{bmatrix}
+
+For example, to get the inertia tensor of the rocket at time 0.5 seconds, you
+can use the following code:
+
+.. jupyter-execute::
+
+    calisto.get_inertia_tensor_at_time(0.5)
+
+Derivative of the Inertia Tensor
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+You can also get the derivative of the inertia tensor at a given time using the
+``get_inertia_tensor_derivative_at_time`` method. Here's an example:
+
+.. jupyter-execute::
+
+    calisto.get_inertia_tensor_derivative_at_time(0.5)
+
+Implications from these results
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The inertia tensor reveals important information about the rocket's symmetry 
+and ease of rotation:
+
+1. **Axis Symmetry**: If I\ :sub:`11` and I\ :sub:`22` are equal, the rocket is symmetric around the axes perpendicular to the rocket's center axis. In our defined rocket, I\ :sub:`11` and I\ :sub:`22` are indeed equal, indicating that our rocket is axisymmetric.
+
+2. **Zero Products of Inertia**: The off-diagonal elements of the inertia tensor are zero, which means the products of inertia are zero. This indicates that the rocket is symmetric around its center axis.
+
+3. **Ease of Rotation**: The I\ :sub:`33` value is significantly lower than the other two. This suggests that the rocket is easier to rotate around its center axis than around the axes perpendicular to the rocket. This is an important factor when considering the rocket's stability and control.
+
+However, these conclusions are based on the assumption that the inertia tensor is calculated with respect to the rocket's center of mass and aligned with the principal axes of the rocket. If the inertia tensor is calculated with respect to a different point or not aligned with the principal axes, the conclusions may not hold.
+

rocketpy/rocket/aero_surface.py

@@ -367,7 +367,7 @@ def evaluate_geometrical_parameters(self):
         return None
 
     def evaluate_nose_shape(self):
-        """Calculates and saves nose cone's shape as lists and reavulates the
+        """Calculates and saves nose cone's shape as lists and re-evaluates the
         nose cone's length for a given bluffness ratio. The shape is saved as
         two vectors, one for the x coordinates and one for the y coordinates.
 
@@ -381,7 +381,9 @@ def evaluate_nose_shape(self):
 
         # Calculate a function to find the tangential intersection point between the circle and nosecone curve.
         def find_x_intercept(x):
-            return x + self.y_nosecone(x) * self.y_nosecone.differentiate(x)
+            return x + self.y_nosecone(x) * self.y_nosecone.differentiate_complex_step(
+                x
+            )
 
         # Calculate a function to find the radius of the nosecone curve
         def find_radius(x):
@@ -770,7 +772,9 @@ def evaluate_lift_coefficient(self):
             )
 
             # Differentiating at alpha = 0 to get cl_alpha
-            clalpha2D_incompressible = self.airfoil_cl.differentiate(x=1e-3, dx=1e-3)
+            clalpha2D_incompressible = self.airfoil_cl.differentiate_complex_step(
+                x=1e-3, dx=1e-3
+            )
 
             # Convert to radians if needed
             if self.airfoil[1] == "degrees":

rocketpy/rocket/rocket.py

@@ -4,6 +4,7 @@
 
 from rocketpy.control.controller import _Controller
 from rocketpy.mathutils.function import Function
+from rocketpy.mathutils.vector_matrix import Matrix
 from rocketpy.motors.motor import EmptyMotor
 from rocketpy.plots.rocket_plots import _RocketPlots
 from rocketpy.prints.rocket_prints import _RocketPrints
@@ -72,6 +73,9 @@ class Rocket:
         for more information
         regarding the coordinate system.
         Expressed in meters as a function of time.
+    Rocket.com_to_cdm_function : Function
+        Function of time expressing the z-coordinate of the center of mass
+        relative to the center of dry mass.
     Rocket.reduced_mass : Function
         Function of time expressing the reduced mass of the rocket,
         defined as the product of the propellant mass and the mass
@@ -146,6 +150,11 @@ class Rocket:
         defined rocket coordinate system.
         See :doc:`Positions and Coordinate Systems </user/positions>`
         for more information.
+    Rocket.nozzle_to_cdm : float
+        Distance between the nozzle exit and the rocket's center of dry mass
+        position, in meters.
+    Rocket.nozzle_gyration_tensor: Matrix
+        Matrix representing the nozzle gyration tensor.
     Rocket.center_of_propellant_position : Function
         Position of the propellant's center of mass relative to the user defined
         rocket reference system. See
@@ -734,8 +743,133 @@ def evaluate_inertias(self):
             self.I_23,
         )
 
+    def evaluate_nozzle_to_cdm(self):
+        """Evaluates the distance between the nozzle exit and the rocket's
+        center of dry mass.
+
+        Returns
+        -------
+        self.nozzle_to_cdm : float
+            Distance between the nozzle exit and the rocket's center of dry
+            mass position, in meters.
+        """
+        self.nozzle_to_cdm = (
+            -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys
+        )
+        return self.nozzle_to_cdm
+
     def evaluate_nozzle_gyration_tensor(self):
-        pass
+        """Calculates and returns the nozzle gyration tensor relative to the
+        rocket's center of dry mass. The gyration tensor is saved and returned
+        in units of kg*m².
+
+        Returns
+        -------
+        self.nozzle_gyration_tensor : Matrix
+            Matrix containing the nozzle gyration tensor.
+        """
+        S_noz_33 = 0.5 * self.motor.nozzle_radius**2
+        S_noz_11 = S_noz_22 = 0.5 * S_noz_33 + 0.25 * self.nozzle_to_cdm**2
+        S_noz_12, S_noz_13, S_noz_23 = 0, 0, 0  # Due to axis symmetry
+        self.nozzle_gyration_tensor = Matrix(
+            [
+                [S_noz_11, S_noz_12, S_noz_13],
+                [S_noz_12, S_noz_22, S_noz_23],
+                [S_noz_13, S_noz_23, S_noz_33],
+            ]
+        )
+        return self.nozzle_gyration_tensor
+
+    def evaluate_com_to_cdm_function(self):
+        """Evaluates the z-coordinate of the center of mass (COM) relative to
+        the center of dry mass (CDM).
+
+        Notes
+        -----
+        1. The `com_to_cdm_function` plus `center_of_mass` should be equal
+        to `center_of_dry_mass_position` at every time step.
+        2. The `com_to_cdm_function` is a function of time and will usually
+        already be discretized.
+
+        Returns
+        -------
+        self.com_to_cdm_function : Function
+            Function of time expressing the z-coordinate of the center of mass
+            relative to the center of dry mass.
+        """
+        self.com_to_cdm_function = (
+            -1
+            * (
+                (self.center_of_propellant_position - self.center_of_dry_mass_position)
+                * self._csys
+            )
+            * self.motor.propellant_mass
+            / self.total_mass
+        )
+        self.com_to_cdm_function.set_inputs("Time (s)")
+        self.com_to_cdm_function.set_outputs("Z Coordinate COM to CDM (m)")
+        self.com_to_cdm_function.set_title("Z Coordinate COM to CDM")
+        return self.com_to_cdm_function
+
+    def get_inertia_tensor_at_time(self, t):
+        """Returns a Matrix representing the inertia tensor of the rocket with
+        respect to the rocket's center of mass at a given time. It evaluates
+        each inertia tensor component at the given time and returns a Matrix
+        with the computed values.
+
+        Parameters
+        ----------
+        t : float
+            Time at which the inertia tensor is to be evaluated.
+
+        Returns
+        -------
+        Matrix
+            Inertia tensor of the rocket at time t.
+        """
+        I_11 = self.I_11.get_value_opt(t)
+        I_12 = self.I_12.get_value_opt(t)
+        I_13 = self.I_13.get_value_opt(t)
+        I_22 = self.I_22.get_value_opt(t)
+        I_23 = self.I_23.get_value_opt(t)
+        I_33 = self.I_33.get_value_opt(t)
+        return Matrix(
+            [
+                [I_11, I_12, I_13],
+                [I_12, I_22, I_23],
+                [I_13, I_23, I_33],
+            ]
+        )
+
+    def get_inertia_tensor_derivative_at_time(self, t):
+        """Returns a Matrix representing the time derivative of the inertia
+        tensor of the rocket with respect to the rocket's center of mass at a
+        given time. It evaluates each inertia tensor component's derivative at
+        the given time and returns a Matrix with the computed values.
+
+        Parameters
+        ----------
+        t : float
+            Time at which the inertia tensor derivative is to be evaluated.
+
+        Returns
+        -------
+        Matrix
+            Inertia tensor time derivative of the rocket at time t.
+        """
+        I_11_dot = self.I_11.differentiate_complex_step(t)
+        I_12_dot = self.I_12.differentiate_complex_step(t)
+        I_13_dot = self.I_13.differentiate_complex_step(t)
+        I_22_dot = self.I_22.differentiate_complex_step(t)
+        I_23_dot = self.I_23.differentiate_complex_step(t)
+        I_33_dot = self.I_33.differentiate_complex_step(t)
+        return Matrix(
+            [
+                [I_11_dot, I_12_dot, I_13_dot],
+                [I_12_dot, I_22_dot, I_23_dot],
+                [I_13_dot, I_23_dot, I_33_dot],
+            ]
+        )
 
     def add_motor(self, motor, position):
         """Adds a motor to the rocket.
@@ -778,6 +912,7 @@ def add_motor(self, motor, position):
         self.evaluate_dry_mass()
         self.evaluate_total_mass()
         self.evaluate_center_of_dry_mass()
+        self.evaluate_nozzle_to_cdm()
         self.evaluate_center_of_mass()
         self.evaluate_dry_inertias()
         self.evaluate_inertias()
@@ -786,6 +921,8 @@ def add_motor(self, motor, position):
         self.evaluate_center_of_pressure()
         self.evaluate_stability_margin()
         self.evaluate_static_margin()
+        self.evaluate_com_to_cdm_function()
+        self.evaluate_nozzle_gyration_tensor()
 
     def add_surfaces(self, surfaces, positions):
         """Adds one or more aerodynamic surfaces to the rocket. The aerodynamic