ENH: Add the Coriolis Force to the Flight class

rocketpy/environment/environment.py

@@ -248,6 +248,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
     -----
@@ -353,6 +355,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)
@@ -584,6 +587,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

@@ -1679,6 +1679,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,
@@ -1745,7 +1751,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
@@ -1904,9 +1910,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),
@@ -1915,6 +1918,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]
 
@@ -1994,6 +2001,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