From 04429a70d22dac3fbed01fdfd53b7d1fbe6f11ad Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:23:51 -0400 Subject: [PATCH 01/17] ENH: add Rocket.evaluate_nozzle_to_center_of_dry_mass_position() --- rocketpy/rocket/rocket.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index c65f1eca1..8d99ee24e 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -734,6 +734,21 @@ def evaluate_inertias(self): self.I_23, ) + def evaluate_nozzle_to_center_of_dry_mass_position(self): + """Evaluates the distance between the nozzle exit and the rocket's + center of dry mass position. + + Returns + ------- + self.nozzle_to_center_of_dry_mass_position : float + Distance between the nozzle exit and the rocket's center of dry + mass position, in meters. + """ + self.nozzle_to_center_of_dry_mass_position = ( + -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys + ) + return self.nozzle_to_center_of_dry_mass_position + def evaluate_nozzle_gyration_tensor(self): pass From 5ef3f04f858c93ccc3b90efa02f9c56d89f924bb Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:24:15 -0400 Subject: [PATCH 02/17] ENH: Add Rocket.evaluate_nozzle_gyration_tensor() --- rocketpy/rocket/rocket.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 8d99ee24e..0afd631c6 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -750,7 +750,28 @@ def evaluate_nozzle_to_center_of_dry_mass_position(self): return self.nozzle_to_center_of_dry_mass_position 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_center_of_dry_mass_position**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 add_motor(self, motor, position): """Adds a motor to the rocket. From cb27b004ce8aa3f1fe3640c979b75dfc73a52cdc Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:24:38 -0400 Subject: [PATCH 03/17] ENH: Adds Rocket.evaluate_z_coordinate_com_to_cdm() --- rocketpy/rocket/rocket.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 0afd631c6..3657c3fb9 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -773,6 +773,36 @@ def evaluate_nozzle_gyration_tensor(self): ) return self.nozzle_gyration_tensor + def evaluate_z_coordinate_com_to_cdm(self): + """Evaluates the z-coordinate of the center of mass (COM) relative to + the center of dry mass (CDM). + + Notes + ----- + 1. The `z_coordinate_com_to_cdm` plus `center_of_mass` should be equal + to `center_of_dry_mass_position` at every time step. + 2. The `z_coordinate_com_to_cdm` is a function of time and will usually + already be discretized. + + Returns + ------- + self.z_coordinate_com_to_cdm : Function + Function of time expressing the z-coordinate of the center of mass + relative to the center of dry mass. + """ + self.z_coordinate_com_to_cdm = ( + -1 + * ( + (self.center_of_propellant_position - self.center_of_dry_mass_position) + * self._csys + ) + * self.motor.propellant_mass + / self.total_mass + ) + self.z_coordinate_com_to_cdm.set_inputs("Time (s)") + self.z_coordinate_com_to_cdm.set_outputs("Z Coordinate COM to CDM (m)") + self.z_coordinate_com_to_cdm.set_title("Z Coordinate COM to CDM") + return self.z_coordinate_com_to_cdm def add_motor(self, motor, position): """Adds a motor to the rocket. From 4da4d19d52fb0697cb96993eb578e22e0c69804d Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:25:35 -0400 Subject: [PATCH 04/17] MNT: updates Rocket.add_motor() to accommodate new evaluations --- rocketpy/rocket/rocket.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 3657c3fb9..57a88c67c 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -844,6 +844,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_center_of_dry_mass_position() self.evaluate_center_of_mass() self.evaluate_dry_inertias() self.evaluate_inertias() @@ -852,6 +853,8 @@ def add_motor(self, motor, position): self.evaluate_center_of_pressure() self.evaluate_stability_margin() self.evaluate_static_margin() + self.evaluate_z_coordinate_com_to_cdm() + self.evaluate_nozzle_gyration_tensor() def add_surfaces(self, surfaces, positions): """Adds one or more aerodynamic surfaces to the rocket. The aerodynamic From 1a00ca1f971e5157c1e093fb479ba73fd6fe16cd Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:26:19 -0400 Subject: [PATCH 05/17] ENH: Adds methods to get inertia tensors at specific time --- rocketpy/rocket/rocket.py | 61 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 57a88c67c..d3d00612a 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -803,6 +803,67 @@ def evaluate_z_coordinate_com_to_cdm(self): self.z_coordinate_com_to_cdm.set_outputs("Z Coordinate COM to CDM (m)") self.z_coordinate_com_to_cdm.set_title("Z Coordinate COM to CDM") return self.z_coordinate_com_to_cdm + + 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. From 74b38b7fc744d8068985e41fe0dcedd037e8e641 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:26:47 -0400 Subject: [PATCH 06/17] MNT: updates docstring of Rocket class --- rocketpy/rocket/rocket.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index d3d00612a..68bd026e8 100644 --- a/rocketpy/rocket/rocket.py +++ b/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.z_coordinate_com_to_cdm : 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 ` for more information. + Rocket.nozzle_to_center_of_dry_mass_position : 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 From fac8a023d77705d103a6da6ec4a825cbcddd304f Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:27:20 -0400 Subject: [PATCH 07/17] DOC: Adds Inertia tensors section to Rocket class page --- docs/user/rocket.rst | 51 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/docs/user/rocket.rst b/docs/user/rocket.rst index 28e091e8d..6b4a9f257 100644 --- a/docs/user/rocket.rst +++ b/docs/user/rocket.rst @@ -428,3 +428,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. + From 9a2f842a06b40267142ff7ce9ea3c327e14dbdfd Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:29:53 -0400 Subject: [PATCH 08/17] ENH: updates Flight.u_dot_generalized to use pre-calculated attributes from Rocket class --- rocketpy/simulation/flight.py | 81 ++++++++--------------------------- 1 file changed, 18 insertions(+), 63 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index b762c376f..10c454e2e 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1570,70 +1570,26 @@ def u_dot_generalized(self, t, u, post_processing=False): w = Vector([omega1, omega2, omega3]) # Angular velocity vector # Retrieve necessary quantities - rho = self.env.density.get_value_opt(z) + ## Rocket mass total_mass = self.rocket.total_mass.get_value_opt(t) total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t) - total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate(t) + total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate_complex_step(t) ## CM position vector and time derivatives relative to CDM in body frame - r_CM_z = ( - -1 - * ( - ( - self.rocket.center_of_propellant_position - - self.rocket.center_of_dry_mass_position - ) - * self.rocket._csys - ) - * self.rocket.motor.propellant_mass - / total_mass - ) - r_CM = Vector([0, 0, r_CM_z.get_value_opt(t)]) - r_CM_dot = Vector([0, 0, r_CM_z.differentiate(t)]) + r_CM_z = self.rocket.z_coordinate_com_to_cdm + r_CM_t = r_CM_z.get_value_opt(t) + r_CM = Vector([0, 0, r_CM_t]) + r_CM_dot = Vector([0, 0, r_CM_z.differentiate_complex_step(t)]) r_CM_ddot = Vector([0, 0, r_CM_z.differentiate(t, order=2)]) + ## Nozzle position vector + r_NOZ = Vector([0, 0, self.rocket.nozzle_to_center_of_dry_mass_position]) ## Nozzle gyration tensor - r_NOZ = ( - -(self.rocket.nozzle_position - self.rocket.center_of_dry_mass_position) - * self.rocket._csys - ) - S_noz_33 = 0.5 * self.rocket.motor.nozzle_radius**2 - S_noz_11 = S_noz_22 = 0.5 * S_noz_33 + 0.25 * r_NOZ**2 - S_noz_12, S_noz_13, S_noz_23 = 0, 0, 0 - S_nozzle = 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], - ] - ) + S_nozzle = self.rocket.nozzle_gyration_tensor ## Inertia tensor - I_11 = self.rocket.I_11.get_value_opt(t) - I_12 = self.rocket.I_12.get_value_opt(t) - I_13 = self.rocket.I_13.get_value_opt(t) - I_22 = self.rocket.I_22.get_value_opt(t) - I_23 = self.rocket.I_23.get_value_opt(t) - I_33 = self.rocket.I_33.get_value_opt(t) - I = Matrix( - [ - [I_11, I_12, I_13], - [I_12, I_22, I_23], - [I_13, I_23, I_33], - ] - ) + I = self.rocket.get_inertia_tensor_at_time(t) ## Inertia tensor time derivative in the body frame - I_11_dot = self.rocket.I_11.differentiate(t) - I_12_dot = self.rocket.I_12.differentiate(t) - I_13_dot = self.rocket.I_13.differentiate(t) - I_22_dot = self.rocket.I_22.differentiate(t) - I_23_dot = self.rocket.I_23.differentiate(t) - I_33_dot = self.rocket.I_33.differentiate(t) - I_dot = 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], - ] - ) - ## Inertia tensor relative to CM + I_dot = self.rocket.get_inertia_tensor_derivative_at_time(t) + + # Calculate the Inertia tensor relative to CM H = (r_CM.cross_matrix @ -r_CM.cross_matrix) * total_mass I_CM = I - H @@ -1723,8 +1679,8 @@ def u_dot_generalized(self, t, u, post_processing=False): R1 += comp_lift_xb R2 += comp_lift_yb # Add to total moment - M1 -= (comp_cpz + r_CM_z.get_value_opt(t)) * comp_lift_yb - M2 += (comp_cpz + r_CM_z.get_value_opt(t)) * comp_lift_xb + M1 -= (comp_cpz + r_CM_t) * comp_lift_yb + M2 += (comp_cpz + r_CM_t) * comp_lift_xb # Calculates Roll Moment try: clf_delta, cld_omega, cant_angle_rad = aero_surface.roll_parameters @@ -1749,13 +1705,12 @@ def u_dot_generalized(self, t, u, post_processing=False): pass weightB = Kt @ Vector([0, 0, -total_mass * self.env.gravity.get_value_opt(z)]) T00 = total_mass * r_CM - r_NOZ_vector = Vector([0, 0, r_NOZ]) - T03 = 2 * total_mass_dot * (r_NOZ_vector - r_CM) - 2 * total_mass * r_CM_dot + T03 = 2 * total_mass_dot * (r_NOZ - r_CM) - 2 * total_mass * r_CM_dot T04 = ( - self.rocket.motor.thrust.get_value_opt(t) * Vector([0, 0, 1]) + Vector([0, 0, self.rocket.motor.thrust.get_value_opt(t)]) - total_mass * r_CM_ddot - 2 * total_mass_dot * r_CM_dot - + total_mass_ddot * (r_NOZ_vector - r_CM) + + total_mass_ddot * (r_NOZ - r_CM) ) T05 = total_mass_dot * S_nozzle - I_dot From 5c31a0b3cfb070bed528feaaa44ee1a88eed70b0 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 00:31:38 -0400 Subject: [PATCH 09/17] DEV: Adds #595 to the CHANGELOG --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fbe00cfa0..dc7c5a6d6 100644 --- a/CHANGELOG.md +++ b/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) From 40a22c4f24dfbd48471eb62368eb3c9c2dda2f67 Mon Sep 17 00:00:00 2001 From: Lint Action Date: Sun, 5 May 2024 04:34:27 +0000 Subject: [PATCH 10/17] Fix code style issues with Black --- rocketpy/rocket/rocket.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 68bd026e8..0bb7f8caf 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -752,12 +752,12 @@ def evaluate_nozzle_to_center_of_dry_mass_position(self): self.nozzle_to_center_of_dry_mass_position : float Distance between the nozzle exit and the rocket's center of dry mass position, in meters. - """ + """ self.nozzle_to_center_of_dry_mass_position = ( -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys ) return self.nozzle_to_center_of_dry_mass_position - + def evaluate_nozzle_gyration_tensor(self): """Calculates and returns the nozzle gyration tensor relative to the rocket's center of dry mass. The gyration tensor is saved and returned From 68d176df04862920a9af702c8252df6243575458 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 01:09:59 -0400 Subject: [PATCH 11/17] TST: add tests for the new Rocket properties --- tests/unit/test_rocket.py | 86 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index 7d942ee60..555add5c4 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -413,3 +413,89 @@ def test_evaluate_center_of_mass(calisto): A predefined instance of the calisto Rocket with a motor, used as a base for testing. """ assert isinstance(calisto.evaluate_center_of_mass(), Function) + + +def test_evaluate_nozzle_to_center_of_dry_mass_position(calisto): + expected_distance = 1.255 + atol = 1e-3 # Equivalent to 1mm + assert ( + pytest.approx(expected_distance, atol) + == calisto.nozzle_to_center_of_dry_mass_position + ) + # Test if calling the function returns the same result + res = calisto.evaluate_nozzle_to_center_of_dry_mass_position() + assert pytest.approx(expected_distance, atol) == res + + +def test_evaluate_nozzle_gyration_tensor(calisto): + expected_gyration_tensor = np.array( + [[0.3940207, 0, 0], [0, 0.3940207, 0], [0, 0, 0.0005445]] + ) + atol = 1e-3 * 1e-2 * 1e-2 # Equivalent to 1g * 1cm^2 + assert np.allclose( + expected_gyration_tensor, np.array(calisto.nozzle_gyration_tensor), atol=atol + ) + # Test if calling the function returns the same result + res = calisto.evaluate_nozzle_gyration_tensor() + assert np.allclose(expected_gyration_tensor, np.array(res), atol=atol) + + +def test_evaluate_z_coordinate_com_to_cdm(calisto): + atol = 1e-3 # Equivalent to 1mm + assert np.allclose( + (calisto.center_of_dry_mass_position - calisto.center_of_mass).source, + calisto.z_coordinate_com_to_cdm.source, + atol=atol, + ) + + +def test_get_inertia_tensor_at_time(calisto): + # Expected values (for t = 0) + # TODO: compute these values by hand or using CAD. + Ixx = 10.31379 + Iyy = 10.31379 + Izz = 0.039942 + + # Set tolerance threshold + atol = 1e-5 + + # Get inertia tensor at t = 0 + inertia_tensor = calisto.get_inertia_tensor_at_time(0) + + # Check if the values are close to the expected ones + assert pytest.approx(Ixx, atol) == inertia_tensor.x[0] + assert pytest.approx(Iyy, atol) == inertia_tensor.y[1] + assert pytest.approx(Izz, atol) == inertia_tensor.z[2] + # Check if products of inertia are zero + assert pytest.approx(0, atol) == inertia_tensor.x[1] + assert pytest.approx(0, atol) == inertia_tensor.x[2] + assert pytest.approx(0, atol) == inertia_tensor.y[0] + assert pytest.approx(0, atol) == inertia_tensor.y[2] + assert pytest.approx(0, atol) == inertia_tensor.z[0] + assert pytest.approx(0, atol) == inertia_tensor.z[1] + + +def test_get_inertia_tensor_derivative_at_time(calisto): + # Expected values (for t = 2s) + # TODO: compute these values by hand or using CAD. + Ixx_dot = -0.634805230901143 + Iyy_dot = -0.634805230901143 + Izz_dot = -0.000671493662305 + + # Set tolerance threshold + atol = 1e-3 + + # Get inertia tensor at t = 2s + inertia_tensor = calisto.get_inertia_tensor_derivative_at_time(2) + + # Check if the values are close to the expected ones + assert pytest.approx(Ixx_dot, atol) == inertia_tensor.x[0] + assert pytest.approx(Iyy_dot, atol) == inertia_tensor.y[1] + assert pytest.approx(Izz_dot, atol) == inertia_tensor.z[2] + # Check if products of inertia are zero + assert pytest.approx(0, atol) == inertia_tensor.x[1] + assert pytest.approx(0, atol) == inertia_tensor.x[2] + assert pytest.approx(0, atol) == inertia_tensor.y[0] + assert pytest.approx(0, atol) == inertia_tensor.y[2] + assert pytest.approx(0, atol) == inertia_tensor.z[0] + assert pytest.approx(0, atol) == inertia_tensor.z[1] From 080e9f031ed7c2cc57761dd9e99a1b33f0c58d7e Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 01:17:15 -0400 Subject: [PATCH 12/17] ENH: Update differentiate_complex_step method to use smaller dx value for more accurate differentiation --- rocketpy/mathutils/function.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 7994c355f..b9efdba95 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -2427,7 +2427,7 @@ def differentiate(self, x, dx=1e-6, order=1): + self.get_value_opt(x - dx) ) / dx**2 - def differentiate_complex_step(self, x, dx=1e-6): + def differentiate_complex_step(self, x, dx=1e-12): """Differentiate a Function object at a given point using the complex step method. This method can be faster than ``Function.differentiate`` since it requires only one evaluation of the function. However, the From b691ead9e63a256024a378a725be1cab11df4669 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sun, 5 May 2024 01:33:09 -0400 Subject: [PATCH 13/17] ENH: use differentiate.complex_time in aero_surface.py --- rocketpy/rocket/aero_surface.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py index d41240ac9..8f0d37f29 100644 --- a/rocketpy/rocket/aero_surface.py +++ b/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": From 7531ca6a58f034613cf9554edc05031b8be130da Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Wed, 15 May 2024 07:16:22 -0400 Subject: [PATCH 14/17] MNT: rename Rocket.nozzle_to_cdm attribute --- rocketpy/rocket/rocket.py | 18 ++++++++---------- rocketpy/simulation/flight.py | 2 +- tests/unit/test_rocket.py | 9 +++------ 3 files changed, 12 insertions(+), 17 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 0bb7f8caf..7d48d8724 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -150,7 +150,7 @@ class Rocket: defined rocket coordinate system. See :doc:`Positions and Coordinate Systems ` for more information. - Rocket.nozzle_to_center_of_dry_mass_position : float + 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 @@ -743,20 +743,20 @@ def evaluate_inertias(self): self.I_23, ) - def evaluate_nozzle_to_center_of_dry_mass_position(self): + def evaluate_nozzle_to_cdm(self): """Evaluates the distance between the nozzle exit and the rocket's - center of dry mass position. + center of dry mass. Returns ------- - self.nozzle_to_center_of_dry_mass_position : float + self.nozzle_to_cdm : float Distance between the nozzle exit and the rocket's center of dry mass position, in meters. """ - self.nozzle_to_center_of_dry_mass_position = ( + self.nozzle_to_cdm = ( -(self.nozzle_position - self.center_of_dry_mass_position) * self._csys ) - return self.nozzle_to_center_of_dry_mass_position + return self.nozzle_to_cdm def evaluate_nozzle_gyration_tensor(self): """Calculates and returns the nozzle gyration tensor relative to the @@ -769,9 +769,7 @@ def evaluate_nozzle_gyration_tensor(self): 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_center_of_dry_mass_position**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( [ @@ -914,7 +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_center_of_dry_mass_position() + self.evaluate_nozzle_to_cdm() self.evaluate_center_of_mass() self.evaluate_dry_inertias() self.evaluate_inertias() diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 10c454e2e..373fd16d0 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1581,7 +1581,7 @@ def u_dot_generalized(self, t, u, post_processing=False): r_CM_dot = Vector([0, 0, r_CM_z.differentiate_complex_step(t)]) r_CM_ddot = Vector([0, 0, r_CM_z.differentiate(t, order=2)]) ## Nozzle position vector - r_NOZ = Vector([0, 0, self.rocket.nozzle_to_center_of_dry_mass_position]) + r_NOZ = Vector([0, 0, self.rocket.nozzle_to_cdm]) ## Nozzle gyration tensor S_nozzle = self.rocket.nozzle_gyration_tensor ## Inertia tensor diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index 555add5c4..40ae2aec0 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -415,15 +415,12 @@ def test_evaluate_center_of_mass(calisto): assert isinstance(calisto.evaluate_center_of_mass(), Function) -def test_evaluate_nozzle_to_center_of_dry_mass_position(calisto): +def test_evaluate_nozzle_to_cdm(calisto): expected_distance = 1.255 atol = 1e-3 # Equivalent to 1mm - assert ( - pytest.approx(expected_distance, atol) - == calisto.nozzle_to_center_of_dry_mass_position - ) + assert pytest.approx(expected_distance, atol) == calisto.nozzle_to_cdm # Test if calling the function returns the same result - res = calisto.evaluate_nozzle_to_center_of_dry_mass_position() + res = calisto.evaluate_nozzle_to_cdm() assert pytest.approx(expected_distance, atol) == res From 2eedf980a92f4534869366bd50043f284b4bdab2 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Wed, 15 May 2024 07:16:42 -0400 Subject: [PATCH 15/17] DOC: adds Inertia Tensor to the toctree --- docs/user/rocket.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/user/rocket.rst b/docs/user/rocket.rst index 6b4a9f257..fa09f058e 100644 --- a/docs/user/rocket.rst +++ b/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. From 1b038d5f3ac97354a480e9aedf29ce7bcebcbf60 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Wed, 15 May 2024 09:02:27 -0400 Subject: [PATCH 16/17] MNT: use rocket.z_coordinate_com_to_cdm in u_dot --- rocketpy/simulation/flight.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 798a08170..68672d802 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1341,12 +1341,8 @@ def u_dot(self, t, u, post_processing=False): ) * self.rocket._csys ) - # c = -self.rocket.distance_rocket_nozzle - c = ( - -(self.rocket.nozzle_position - self.rocket.center_of_dry_mass_position) - * self.rocket._csys - ) - a = b * Mt / M + c = self.rocket.nozzle_to_cdm + a = self.rocket.z_coordinate_com_to_cdm.get_value_opt(t) rN = self.rocket.motor.nozzle_radius # Prepare transformation matrix a11 = 1 - 2 * (e2**2 + e3**2) From 61ebbc754dc41a7436b9ddc5ccc39310fecc4440 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> Date: Thu, 16 May 2024 20:30:21 +0000 Subject: [PATCH 17/17] MNT: renames z_coordinate_com_to_cdm to com_to_cdm_function --- rocketpy/rocket/rocket.py | 22 +++++++++++----------- rocketpy/simulation/flight.py | 4 ++-- tests/unit/test_rocket.py | 4 ++-- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index c20a35ef2..a0db26eff 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -73,7 +73,7 @@ class Rocket: for more information regarding the coordinate system. Expressed in meters as a function of time. - Rocket.z_coordinate_com_to_cdm : Function + 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 @@ -780,24 +780,24 @@ def evaluate_nozzle_gyration_tensor(self): ) return self.nozzle_gyration_tensor - def evaluate_z_coordinate_com_to_cdm(self): + 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 `z_coordinate_com_to_cdm` plus `center_of_mass` should be equal + 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 `z_coordinate_com_to_cdm` is a function of time and will usually + 2. The `com_to_cdm_function` is a function of time and will usually already be discretized. Returns ------- - self.z_coordinate_com_to_cdm : Function + 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.z_coordinate_com_to_cdm = ( + self.com_to_cdm_function = ( -1 * ( (self.center_of_propellant_position - self.center_of_dry_mass_position) @@ -806,10 +806,10 @@ def evaluate_z_coordinate_com_to_cdm(self): * self.motor.propellant_mass / self.total_mass ) - self.z_coordinate_com_to_cdm.set_inputs("Time (s)") - self.z_coordinate_com_to_cdm.set_outputs("Z Coordinate COM to CDM (m)") - self.z_coordinate_com_to_cdm.set_title("Z Coordinate COM to CDM") - return self.z_coordinate_com_to_cdm + 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 @@ -921,7 +921,7 @@ def add_motor(self, motor, position): self.evaluate_center_of_pressure() self.evaluate_stability_margin() self.evaluate_static_margin() - self.evaluate_z_coordinate_com_to_cdm() + self.evaluate_com_to_cdm_function() self.evaluate_nozzle_gyration_tensor() def add_surfaces(self, surfaces, positions): diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 68672d802..4cbc752fb 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -1342,7 +1342,7 @@ def u_dot(self, t, u, post_processing=False): * self.rocket._csys ) c = self.rocket.nozzle_to_cdm - a = self.rocket.z_coordinate_com_to_cdm.get_value_opt(t) + a = self.rocket.com_to_cdm_function.get_value_opt(t) rN = self.rocket.motor.nozzle_radius # Prepare transformation matrix a11 = 1 - 2 * (e2**2 + e3**2) @@ -1571,7 +1571,7 @@ def u_dot_generalized(self, t, u, post_processing=False): total_mass_dot = self.rocket.total_mass_flow_rate.get_value_opt(t) total_mass_ddot = self.rocket.total_mass_flow_rate.differentiate_complex_step(t) ## CM position vector and time derivatives relative to CDM in body frame - r_CM_z = self.rocket.z_coordinate_com_to_cdm + r_CM_z = self.rocket.com_to_cdm_function r_CM_t = r_CM_z.get_value_opt(t) r_CM = Vector([0, 0, r_CM_t]) r_CM_dot = Vector([0, 0, r_CM_z.differentiate_complex_step(t)]) diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index 40ae2aec0..4d934efef 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -437,11 +437,11 @@ def test_evaluate_nozzle_gyration_tensor(calisto): assert np.allclose(expected_gyration_tensor, np.array(res), atol=atol) -def test_evaluate_z_coordinate_com_to_cdm(calisto): +def test_evaluate_com_to_cdm_function(calisto): atol = 1e-3 # Equivalent to 1mm assert np.allclose( (calisto.center_of_dry_mass_position - calisto.center_of_mass).source, - calisto.z_coordinate_com_to_cdm.source, + calisto.com_to_cdm_function.source, atol=atol, )