feat: Implement low-precision sun/moon ephemerides and third-body gra…

…vity

src/orbit_dynamics/ephemerides.rs

@@ -0,0 +1,143 @@
+/*!
+Provide low-accuracy ephemerides for various celestial bodies.
+ */
+
+
+use nalgebra::Vector3;
+
+use crate::attitude::RotationMatrix;
+use crate::constants::{AS2RAD, MJD2000};
+use crate::time::{Epoch, TimeSystem};
+
+/// Calculate the position of the Sun in the EME2000 inertial frame using low-precision analytical
+/// methods. For most purposes the EME2000 inertial frame is equivalent to the GCRF frame.
+///
+/// # Arguments
+///
+/// - `epc`: Epoch at which to calculate the Sun's position
+///
+/// # Returns
+///
+/// - `r`: Position of the Sun in the J2000 ecliptic frame. Units: [m]
+///
+/// # References
+///
+/// - O. Montenbruck, and E. Gill, *Satellite Orbits: Models, Methods and Applications*, 2012, p.70-73.
+///
+/// # Example
+///
+/// ```
+/// use brahe::ephemerides::sun_position;
+/// use brahe::time::Epoch;
+/// use brahe::TimeSystem;
+/// use brahe::constants::AU;
+///
+/// let epc = Epoch::from_date(2024, 2, 25, TimeSystem::UTC);
+///
+/// let r = sun_position(epc);
+/// ```
+pub fn sun_position(epc: Epoch) -> Vector3<f64> {
+    // Constants
+    let pi = std::f64::consts::PI;
+    let mjd_tt = epc.mjd_as_time_system(TimeSystem::TT);
+    let epsilon = 23.43929111 * pi * 180.0; // Obliquity of J2000 ecliptic
+    let T = (mjd_tt - MJD2000) / 36525.0; // Julian cent. since J2000
+
+    // Variables
+
+    // Mean anomaly, ecliptic longitude and radius
+    let M = 2.0 * pi * (0.9931267 + 99.9973583 * T).fract(); // [rad]
+    let L = 2.0 * pi * (0.7859444 + M / (2.0 * pi).fract() + (6892.0 * M.sin() + 72.0 * (2.0 * M).sin()) / 1296.0e3); // [rad]
+    let r = 149.619e9 - 2.499e9 * M.cos() - 0.021e9 * (2.0 * M).cos(); // [m]
+
+    // Equatorial position vector
+    RotationMatrix::Rx(-epsilon, false) * Vector3::new(r * L.cos(), r * L.sin(), 0.0)
+}
+
+/// Calculate the position of the Moon in the EME2000 inertial frame using low-precision analytical
+/// methods. For most purposes the EME2000 inertial frame is equivalent to the GCRF frame.
+///
+/// # Arguments
+///
+/// - `epc`: Epoch at which to calculate the Moon's position
+///
+/// # Returns
+///
+/// - `r`: Position of the Moon in the J2000 ecliptic frame. Units: [m]
+///
+/// # References
+///
+/// - O. Montenbruck, and E. Gill, *Satellite Orbits: Models, Methods and Applications*, 2012, p.70-73.
+///
+/// # Example
+///
+/// ```
+/// use brahe::ephemerides::moon_position;
+/// use brahe::time::{Epoch, TimeSystem};
+/// use brahe::constants::AU;
+///
+/// let epc = Epoch::from_date(2024, 2, 25, TimeSystem::UTC);
+///
+/// let r = moon_position(epc);
+/// ```
+pub fn moon_position(epc: Epoch) -> Vector3<f64> {
+    // Constants
+    let pi = std::f64::consts::PI;
+    let mjd_tt = epc.mjd_as_time_system(TimeSystem::TT);
+    let epsilon = 23.43929111 * pi * 180.0; // Obliquity of J2000 ecliptic
+    let T = (mjd_tt - MJD2000) / 36525.0; // Julian cent. since J2000
+
+    // Mean elements of lunar orbit
+    let L_0 = (0.606433 + 1336.851344 * T).fract();       // Mean longitude [rev] w.r.t. J2000 equinox
+    let l = 2.0 * pi * (0.374897 + 1325.552410 * T).fract(); // Moon's mean anomaly [rad]
+    let lp = 2.0 * pi * (0.993133 + 99.997361 * T).fract(); // Sun's mean anomaly [rad]
+    let D = 2.0 * pi * (0.827361 + 1236.853086 * T).fract(); // Diff. long. Moon-Sun [rad]
+    let F = 2.0 * pi * (0.259086 + 1342.227825 * T).fract(); // Argument of latitude
+
+    // Ecliptic longitude (w.r.t. equinox of J2000)
+    let dL = 22640.0 * l.sin() - 4586.0 * (l - 2.0 * D).sin() + 2370.0 * (2.0 * D).sin() + 769.0 * (2.0 * l).sin()
+        - 668.0 * (lp).sin() - 412.0 * (2.0 * F).sin() - 212.0 * (2.0 * l - 2.0 * D).sin() - 206.0 * (l + lp - 2.0 * D).sin()
+        + 192.0 * (l + 2.0 * D).sin() - 165.0 * (lp - 2.0 * D).sin() - 125.0 * D.sin() - 110.0 * (l + lp).sin()
+        + 148.0 * (l - lp).sin() - 55.0 * (2.0 * F - 2.0 * D).sin();
+
+    let L = 2.0 * pi * (L_0 + dL / 1296.0e3).fract(); // [rad]
+
+    // Ecliptic latitude
+    let S = F + (dL + 412.0 * (2.0 * F).sin() + 541.0 * lp.sin()) * AS2RAD;
+    let h = F - 2.0 * D;
+    let N = -526.0 * h.sin() + 44.0 * (l + h).sin() - 31.0 * (-l + h).sin() - 23.0 * (lp + h).sin()
+        + 11.0 * (-lp + h).sin() - 25.0 * (-2.0 * l + F).sin() + 21.0 * (-l + F).sin();
+    let B = (18520.0 * S.sin() + N) * AS2RAD;   // [rad]
+
+    // Distance [m]
+    let r = 385000e3 - 20905e3 * l.cos() - 3699e3 * (2.0 * D - l).cos() - 2956e3 * (2.0 * D).cos()
+        - 570e3 * (2.0 * l).cos() + 246e3 * (2.0 * l - 2.0 * D).cos() - 205e3 * (lp - 2.0 * D).cos()
+        - 171e3 * (l + 2.0 * D).cos() - 152e3 * (l + lp - 2.0 * D).cos();
+
+    // Equatorial coordinates
+    RotationMatrix::Rx(-epsilon, false) * Vector3::new(r * L.cos() * B.cos(), r * L.sin() * B.cos(), r * B.sin())
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::time::Epoch;
+    use crate::TimeSystem;
+
+    use super::*;
+
+    #[test]
+    fn test_sun_position() {
+        let epc = Epoch::from_date(2024, 2, 25, TimeSystem::UTC);
+        let r = sun_position(epc);
+
+        // TODO: Validate algorithm with known values
+    }
+
+    #[test]
+    fn test_moon_position() {
+        let epc = Epoch::from_date(2024, 2, 25, TimeSystem::UTC);
+        let r = moon_position(epc);
+
+        // TODO: Validate algorithm with known values
+    }
+}

src/orbit_dynamics/gravity.rs

@@ -479,9 +479,9 @@ impl std::fmt::Debug for GravityModel {
 /// let r_eci: Vector3<f64> = x_eci.fixed_rows::<3>(0).into();
 ///
 /// // Compute the acceleration due to gravity
-/// let a_grav = brahe::gravity::accel_gravity_spherical_harmonics(&r_eci, &R_i2b, &gravity_model, 20, 20);
+/// let a_grav = brahe::gravity::acceleration_gravity_spherical_harmonics(&r_eci, &R_i2b, &gravity_model, 20, 20);
 /// ```
-pub fn accel_gravity_spherical_harmonics(
+pub fn acceleration_gravity_spherical_harmonics(
     r_eci: &Vector3<f64>,
     R_i2b: &Matrix3<f64>,
     gravity_model: &GravityModel,

src/orbit_dynamics/mod.rs

@@ -1,13 +1,16 @@
 /*!
- Module implementing orbit dynamics models.
-*/
+Module implementing orbit dynamics models.
+ */
+
+pub use ephemerides::*;
+pub use gravity::*;
+pub use relativity::*;
+pub use solar_radiation_pressure::*;
+pub use third_body::*;
 
 pub mod gravity;
 pub mod solar_radiation_pressure;
 pub mod relativity;
 pub mod third_body;
+pub mod ephemerides;
 
-pub use gravity::*;
-pub use solar_radiation_pressure::*;
-pub use third_body::*;
-pub use relativity::*;

src/orbit_dynamics/third_body.rs

@@ -1,4 +1,121 @@
 /*!
- Module for the third body perturbations. Also provides low-precession models for the Sun and Moon
- ephemerides.
-*/
+Module for the third body perturbations. Also provides low-precession models for the Sun and Moon
+ephemerides.
+ */
+
+use nalgebra::Vector3;
+
+use crate::{GM_MOON, GM_SUN};
+use crate::ephemerides::{moon_position, sun_position};
+use crate::orbit_dynamics::gravity::acceleration_point_mass_gravity;
+use crate::time::Epoch;
+
+/// Calculate the acceleration due to the Sun on an object at a given epoch.
+/// The calculation is performed using the point-mass gravity model and the
+/// low-precision analytical ephemerides for the Sun position implemented in
+/// the `ephemerides` module.
+///
+/// Should a more accurate calculation be required, you can utilize the
+/// point-mass gravity model and a higher-precision ephemerides for the Sun.
+///
+/// # Arguments
+///
+/// * `epc` - Epoch at which to calculate the Sun's position
+/// * `r_object` - Position of the object in the GCRF frame. Units: [m]
+///
+/// # Returns
+///
+/// * `a` - Acceleration due to the Sun. Units: [m/s^2]
+///
+/// # Example
+///
+/// ```
+/// use brahe::eop::{set_global_eop_provider, FileEOPProvider, EOPExtrapolation};
+/// use brahe::time::Epoch;
+/// use brahe::third_body::third_body_sun;
+/// use brahe::constants::R_EARTH;
+/// use nalgebra::Vector3;
+///
+/// let eop = FileEOPProvider::from_default_standard(true, EOPExtrapolation::Hold).unwrap();
+/// set_global_eop_provider(eop);
+///
+/// let epc = Epoch::from_date(2024, 2, 25, brahe::TimeSystem::UTC);
+/// let r_object = Vector3::new(R_EARTH + 500e3, 0.0, 0.0);
+///
+/// let a = third_body_sun(epc, &r_object);
+/// ```
+pub fn third_body_sun(epc: Epoch, r_object: &Vector3<f64>) -> Vector3<f64> {
+    acceleration_point_mass_gravity(r_object, &sun_position(epc), GM_SUN)
+}
+
+/// Calculate the acceleration due to the Moon on an object at a given epoch.
+/// The calculation is performed using the point-mass gravity model and the
+/// low-precision analytical ephemerides for the Moon position implemented in
+/// the `ephemerides` module.
+///
+/// Should a more accurate calculation be required, you can utilize the
+/// point-mass gravity model and a higher-precision ephemerides for the Moon.
+///
+/// # Arguments
+///
+/// - `epc` - Epoch at which to calculate the Moon's position
+/// - `r_object` - Position of the object in the GCRF frame. Units: [m]
+///
+/// # Returns
+///
+/// - `a` - Acceleration due to the Moon. Units: [m/s^2]
+///
+pub fn third_body_moon(epc: Epoch, r_object: &Vector3<f64>) -> Vector3<f64> {
+    acceleration_point_mass_gravity(r_object, &moon_position(epc), GM_MOON)
+}
+
+#[cfg(test)]
+mod tests {
+    use nalgebra::Vector6;
+
+    use crate::constants::R_EARTH;
+    use crate::coordinates::*;
+    use crate::TimeSystem;
+
+    use super::*;
+
+    #[test]
+    fn test_third_body_sun() {
+        let epc = Epoch::from_date(2023, 1, 1, TimeSystem::UTC);
+
+        let oe = Vector6::new(
+            R_EARTH + 500e3,
+            0.01,
+            97.3,
+            15.0,
+            30.0,
+            45.0,
+        );
+        let r_object = state_osculating_to_cartesian(oe, true).xyz();
+
+        let a = third_body_sun(epc, &r_object);
+
+        // TODO: Do better validation of the implementation and the expected results
+        assert!(a.norm() < 1.0e-1);
+    }
+
+    #[test]
+    fn test_third_body_moon() {
+        let epc = Epoch::from_date(2023, 1, 1, TimeSystem::UTC);
+
+        let oe = Vector6::new(
+            R_EARTH + 500e3,
+            0.01,
+            97.3,
+            15.0,
+            30.0,
+            45.0,
+        );
+        let r_object = state_osculating_to_cartesian(oe, true).xyz();
+
+        let a = third_body_moon(epc, &r_object);
+
+        // TODO: Do better validation of the implementation and the expected results
+        assert!(a.norm() < 1.0e-4);
+    }
+}