ENH: Add the Coriolis Force to the Flight class (#799)

* wind factor bug corrected

the wind factor wasn't applied to the env.wind_velocity properties

* BUG: StochasticModel visualize attributes of a uniform distribution

It showed the nominal and the standard deviation values and it doesn't make sense in a uniform distribution. In a np.random.uniform the 'nominal value' is the lower bound of the distribution, and the 'standard deviation' value is the upper bound. Now, a new condition has been added for the uniform distributions where the mean and semi range are calculated and showed. This way the visualize_attribute function will show the whole range where the random values are uniformly taken in

* variable names corrections

* Corrections requested by the pylint test

* ENH: add multiplication for 2D functions in rocketpy.function

Added the ability to multiply functions with 2D domains in the __mul__ function

* ENH: StochasticAirBrakes class created

The StochasticAirBrakes class has been created. The __init__.py files in the stochastic and rocketpy folders have also been modified accordingly to incorporate this new class

* ENH: set_air_brakes function created

This functions appends an airbrake and controller objects previuosly created to the rocket

* ENH: add StochasticAirBrake to rocketpy.stochastic_rocket

Some functions has been modified and other has been created in order to include the new StochasticAirBrakes feature into the StochasticRocket class. A new function named 'add_air_brakes' has been created to append a StochasticAirBrakes and Controller objects to the StochasticRocket object. A new function '_create_air_brake' has been introduced to create a sample of an AirBrake object through a StochasticAirBrake object. Enventually, the 'create_object' function has been modified to add the sampled AirBrakes to the sampled Rocket

* BUG: StochasticAirBrake object input in _Controller

When defining the _Controller object a StochasticAirBrake was input. This is already corrected and a AirBrake object is now introduced

* ENH: add time_overshoot option to rocketpy.stochastic_flight

Since the new StochasticAirBrake class is defined, we need the 'time_overshoot' option in the Flight class to ensure that the time step defined in the simulation is the controller sampling rate. The MonteCarlo class has had to be modified as well to include this option.

* DOC: StochasticAirBrakes related documentation added

Documentation related to the StochasticAirBrakes implementation has been added in StochasticAirBrakes, StochasticRocket and Rocket classes.

* ENH: pylint recommendations done

* ENH: Reformatted files to pass Ruff linting checks

* ENH: Update rocketpy/stochastic/stochastic_rocket.py

Unnecessary comment

Co-authored-by: Gui-FernandesBR <[email protected]>

* ENH: more intuitive uniform distribution display in StochasticModel

Co-authored-by: MateusStano <[email protected]>

* DOC: improve drag curve factor definition in StochasticAirBrakes

* ENH: Change assert statement to if

Co-authored-by: Gui-FernandesBR <[email protected]>

* DOC: better explanation of __mul__ function

Co-authored-by: MateusStano <[email protected]>

* ENH: delete set_air_brakes function for simplicity

* ENH: inertial foreces added to u_dot_generalized

* ENH: define Earth's angular velocity vector in Environment

* ENH: some corrections to the Flight class

* ENH: modifications in the Flight class

* DOC: improving Environment documentation

* DOC: more improvements in the Environment class

* ENH: format changes done

* DOC: env.earth_rotation_vector improved

Co-authored-by: Gui-FernandesBR <[email protected]>

* ENH: Coriolis acceleration added to u_dot

* BUG: print left

* ENH: ruff changes

* ENH: CHANGELOG updated

* ENH: remove unecessary frame rotation

* ENH: remove rotation from solid prop udot

* ENH: add coriolis to parachute

* TST: fix tests values

* MNT: remove debug functions

* DEV: changelog

---------

Co-authored-by: Gui-FernandesBR <[email protected]>
Co-authored-by: MateusStano <[email protected]>
Co-authored-by: MateusStano <[email protected]>
Co-authored-by: MateusStano <[email protected]>

Author: 4 people

CHANGELOG.md

@@ -32,11 +32,11 @@ Attention: The newest changes should be on top -->
 
 ### Added
 
+-  ENH: Add the Coriolis Force to the Flight class [#799](https://github.com/RocketPy-Team/RocketPy/pull/799)
 
 ### Changed
 
 - ENH: _MotorPrints inheritance - issue #460 [#828](https://github.com/RocketPy-Team/RocketPy/pull/828)
-
 - MNT: fix deprecations and warnings [#829](https://github.com/RocketPy-Team/RocketPy/pull/829)
 
 ### Fixed

rocketpy/environment/environment.py

@@ -247,6 +247,8 @@ class Environment:
         Number of ensemble members. Only defined when using Ensembles.
     Environment.ensemble_member : int
         Current selected ensemble member. Only defined when using Ensembles.
+    Environment.earth_rotation_vector : list[float]
+        Earth's angular velocity vector in the Flight Coordinate System.
 
     Notes
     -----
@@ -352,6 +354,7 @@ def __init__(
         self.set_location(latitude, longitude)
         self.__initialize_earth_geometry(datum)
         self.__initialize_utm_coordinates()
+        self.__set_earth_rotation_vector()
 
         # Set the gravity model
         self.gravity = self.set_gravity_model(gravity)
@@ -583,6 +586,23 @@ def __reset_wind_direction_function(self):
         self.wind_direction.set_outputs("Wind Direction (Deg True)")
         self.wind_direction.set_title("Wind Direction Profile")
 
+    def __set_earth_rotation_vector(self):
+        """Calculates and stores the Earth's angular velocity vector in the Flight
+        Coordinate System, which is essential for evaluating inertial forces.
+        """
+        # Sidereal day
+        T = 86164.1  # seconds
+
+        # Earth's angular velocity magnitude
+        w_earth = 2 * np.pi / T
+
+        # Vector in the Flight Coordinate System
+        lat = np.radians(self.latitude)
+        w_local = [0, w_earth * np.cos(lat), w_earth * np.sin(lat)]
+
+        # Store the attribute
+        self.earth_rotation_vector = w_local
+
     # Validators (used to verify an attribute is being set correctly.)
 
     def __validate_dictionary(self, file, dictionary):

rocketpy/simulation/flight.py

@@ -1676,6 +1676,12 @@ def u_dot(self, t, u, post_processing=False):  # pylint: disable=too-many-locals
         ax, ay, az = K @ Vector(L)
         az -= self.env.gravity.get_value_opt(z)  # Include gravity
 
+        # Coriolis acceleration
+        _, w_earth_y, w_earth_z = self.env.earth_rotation_vector
+        ax -= 2 * (vz * w_earth_y - vy * w_earth_z)
+        ay -= 2 * (vx * w_earth_z)
+        az -= 2 * (-vx * w_earth_y)
+
         # Create u_dot
         u_dot = [
             vx,
@@ -1742,7 +1748,7 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
         _, _, z, vx, vy, vz, e0, e1, e2, e3, omega1, omega2, omega3 = u
 
         # Create necessary vectors
-        # r = Vector([x, y, z])               # CDM position vector
+        # r = Vector([x, y, z])  # CDM position vector
         v = Vector([vx, vy, vz])  # CDM velocity vector
         e = [e0, e1, e2, e3]  # Euler parameters/quaternions
         w = Vector([omega1, omega2, omega3])  # Angular velocity vector
@@ -1901,9 +1907,6 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
         # Angular velocity derivative
         w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM))
 
-        # Velocity vector derivative
-        v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot))
-
         # Euler parameters derivative
         e_dot = [
             0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3),
@@ -1912,6 +1915,10 @@ def u_dot_generalized(self, t, u, post_processing=False):  # pylint: disable=too
             0.5 * (omega3 * e0 + omega2 * e1 - omega1 * e2),
         ]
 
+        # Velocity vector derivative + Coriolis acceleration
+        w_earth = Vector(self.env.earth_rotation_vector)
+        v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v)
+
         # Position vector derivative
         r_dot = [vx, vy, vz]
 
@@ -1991,6 +1998,12 @@ def u_dot_parachute(self, t, u, post_processing=False):
         ay = Dy / (mp + ma)
         az = (Dz - 9.8 * mp) / (mp + ma)
 
+        # Add coriolis acceleration
+        _, w_earth_y, w_earth_z = self.env.earth_rotation_vector
+        ax -= 2 * (vz * w_earth_y - vy * w_earth_z)
+        ay -= 2 * (vx * w_earth_z)
+        az -= 2 * (-vx * w_earth_y)
+
         if post_processing:
             self.__post_processed_variables.append(
                 [t, ax, ay, az, 0, 0, 0, Dx, Dy, Dz, 0, 0, 0, 0]

tests/integration/test_flight.py

@@ -491,8 +491,9 @@ def test_empty_motor_flight(mock_show, example_plain_env, calisto_motorless):  #
 
 def test_freestream_speed_at_apogee(example_plain_env, calisto):
     """
-    Asserts that a rocket at apogee has a free stream speed of 0.0 m/s in all
-    directions given that the environment doesn't have any wind.
+    Asserts that a rocket at apogee has a free stream speed of near 0.0 m/s
+    in all directions given that the environment doesn't have any wind. Any
+    speed values comes from coriolis effect.
     """
     # NOTE: this rocket doesn't move in x or z direction. There's no wind.
     hard_atol = 1e-12
@@ -508,7 +509,9 @@ def test_freestream_speed_at_apogee(example_plain_env, calisto):
     )
 
     assert np.isclose(
-        test_flight.stream_velocity_x(test_flight.apogee_time), 0.0, atol=hard_atol
+        test_flight.stream_velocity_x(test_flight.apogee_time),
+        0.4641492104717301,
+        atol=hard_atol,
     )
     assert np.isclose(
         test_flight.stream_velocity_y(test_flight.apogee_time), 0.0, atol=hard_atol
@@ -518,9 +521,13 @@ def test_freestream_speed_at_apogee(example_plain_env, calisto):
         test_flight.stream_velocity_z(test_flight.apogee_time), 0.0, atol=soft_atol
     )
     assert np.isclose(
-        test_flight.free_stream_speed(test_flight.apogee_time), 0.0, atol=soft_atol
+        test_flight.free_stream_speed(test_flight.apogee_time),
+        0.4641492104717798,
+        atol=hard_atol,
+    )
+    assert np.isclose(
+        test_flight.apogee_freestream_speed, 0.4641492104717798, atol=hard_atol
     )
-    assert np.isclose(test_flight.apogee_freestream_speed, 0.0, atol=soft_atol)
 
 
 def test_rocket_csys_equivalence(
@@ -546,6 +553,7 @@ def test_rocket_csys_equivalence(
     assert np.isclose(
         flight_calisto_robust.x_impact,
         flight_calisto_nose_to_tail_robust.x_impact,
+        atol=1e-3,
     )
     assert np.isclose(
         flight_calisto_robust.y_impact,

tests/unit/test_flight.py

@@ -109,15 +109,15 @@ def test_get_solution_at_time(flight_calisto):
         flight_calisto.get_solution_at_time(flight_calisto.t_final),
         np.array(
             [
-                48.4313533,
-                0.0,
-                985.7665845,
-                -0.00000229951048,
-                0.0,
-                11.2223284,
-                -341.028803,
-                0.999048222,
-                -0.0436193874,
+                48.43719482805657,
+                -14.836008075478597,
+                985.9858934483618,
+                -3.4415459237894554e-05,
+                0.0007572309307800201,
+                11.21695000766671,
+                -341.1460775169661,
+                0.9990482215818578,
+                -0.043619387365336,
                 0.0,
                 0.0,
                 0.0,
@@ -212,7 +212,7 @@ def test_export_sensor_data(flight_calisto_with_sensors):
     [
         ("t_initial", (0.25886, -0.649623, 0)),
         ("out_of_rail_time", (0.792028, -1.987634, 0)),
-        ("apogee_time", (-0.522875, -0.741825, 0)),
+        ("apogee_time", (-0.509420, -0.732933, -2.089120e-14)),
         ("t_final", (0, 0, 0)),
     ],
 )
@@ -251,8 +251,8 @@ def test_aerodynamic_moments(flight_calisto_custom_wind, flight_time, expected_v
     [
         ("t_initial", (1.654150, 0.659142, -0.067103)),
         ("out_of_rail_time", (5.052628, 2.013361, -1.75370)),
-        ("apogee_time", (2.339424, -1.648934, -0.938867)),
-        ("t_final", (0, 0, 159.2210)),
+        ("apogee_time", (2.321838, -1.613641, -0.962108)),
+        ("t_final", (-0.025792, 0.012030, 159.202481)),
     ],
 )
 def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_values):
@@ -292,7 +292,7 @@ def test_aerodynamic_forces(flight_calisto_custom_wind, flight_time, expected_va
         ("out_of_rail_time", (0, 2.248540, 25.700928)),
         (
             "apogee_time",
-            (-14.488364, 15.638049, -0.000191),
+            (-14.826350, 15.670022, -0.000264),
         ),
         ("t_final", (5, 2, -5.660155)),
     ],
@@ -540,7 +540,7 @@ def test_lat_lon_conversion_from_origin(mock_show, example_plain_env, calisto_ro
         heading=0,
     )
 
-    assert abs(test_flight.longitude(test_flight.t_final) - 0) < 1e-12
+    assert abs(test_flight.longitude(test_flight.t_final)) < 1e-4
     assert test_flight.latitude(test_flight.t_final) > 0