lox_bodies/

lib.rs

// SPDX-FileCopyrightText: 2023 Angus Morrison <github@angus-morrison.com>
// SPDX-FileCopyrightText: 2023 Helge Eichhorn <git@helgeeichhorn.de>
//
// SPDX-License-Identifier: MPL-2.0

#![warn(missing_docs)]

//! Definitions for celestial bodies, barycenters, and minor bodies in the solar system.
//!
//! Body properties such as gravitational parameters, radii, and rotational elements are derived
//! from the IAU/NAIF SPICE planetary constants kernel (PCK).

pub use crate::dynamic::DynOrigin;
pub use generated::*;
use lox_core::elements::GravitationalParameter;
use lox_core::f64::consts::{SECONDS_PER_DAY, SECONDS_PER_JULIAN_CENTURY};
use lox_core::units::Distance;
use std::fmt::{Display, Formatter};
use thiserror::Error;

/// Dynamic dispatch variants of origin types.
pub mod dynamic;
#[allow(clippy::approx_constant)]
mod generated;

/// A NAIF ID code identifying a celestial body or barycenter.
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
#[repr(transparent)]
pub struct NaifId(pub i32);

impl Display for NaifId {
    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
        write!(f, "{}", self.0)
    }
}

/// `Origin` is implemented for all bodies and barycenters.
pub trait Origin {
    /// Returns the NAIF ID of the origin.
    fn id(&self) -> NaifId;
    /// Returns the name of the origin.
    fn name(&self) -> &'static str;
}

/// Error returned when a property is not defined for a given origin.
#[derive(Clone, Debug, Error, Eq, PartialEq)]
#[error("undefined property '{prop}' for origin '{origin}'")]
pub struct UndefinedOriginPropertyError {
    origin: String,
    prop: String,
}

/// Equatorial-a, equatorial-b, and polar radii of a triaxial ellipsoid.
pub type Radii = (Distance, Distance, Distance);

/// Fallible accessor for the triaxial ellipsoid radii of a body.
pub trait TryTriaxialEllipsoid: Origin {
    /// Returns the triaxial ellipsoid radii, or an error if undefined.
    fn try_radii(&self) -> Result<Radii, UndefinedOriginPropertyError>;
}

/// Infallible accessor for the triaxial ellipsoid radii of a body.
pub trait TriaxialEllipsoid: Origin {
    /// Returns the triaxial ellipsoid radii.
    fn radii(&self) -> Radii;
}

impl<T: TriaxialEllipsoid> TryTriaxialEllipsoid for T {
    fn try_radii(&self) -> Result<Radii, UndefinedOriginPropertyError> {
        Ok(self.radii())
    }
}

fn flattening(equatorial_radius: Distance, polar_radius: Distance) -> f64 {
    let r_eq = equatorial_radius.as_f64();
    let r_p = polar_radius.as_f64();
    (r_eq - r_p) / r_eq
}

/// Infallible accessor for the spheroid properties of a body.
pub trait Spheroid: TriaxialEllipsoid {
    /// Returns the equatorial radius.
    fn equatorial_radius(&self) -> Distance {
        self.radii().0
    }

    /// Returns the polar radius.
    fn polar_radius(&self) -> Distance {
        self.radii().2
    }

    /// Returns the flattening factor.
    fn flattening(&self) -> f64 {
        flattening(self.equatorial_radius(), self.polar_radius())
    }
}

/// Fallible accessor for the spheroid properties of a body.
pub trait TrySpheroid: TryTriaxialEllipsoid {
    /// Returns the equatorial radius, or an error if undefined.
    fn try_equatorial_radius(&self) -> Result<Distance, UndefinedOriginPropertyError>;

    /// Returns the polar radius, or an error if undefined.
    fn try_polar_radius(&self) -> Result<Distance, UndefinedOriginPropertyError>;

    /// Returns the flattening factor, or an error if undefined.
    fn try_flattening(&self) -> Result<f64, UndefinedOriginPropertyError> {
        self.try_radii().map(|radii| flattening(radii.0, radii.2))
    }
}

impl<T: Spheroid> TrySpheroid for T {
    fn try_equatorial_radius(&self) -> Result<Distance, UndefinedOriginPropertyError> {
        Ok(self.equatorial_radius())
    }

    fn try_polar_radius(&self) -> Result<Distance, UndefinedOriginPropertyError> {
        Ok(self.polar_radius())
    }

    fn try_flattening(&self) -> Result<f64, UndefinedOriginPropertyError> {
        Ok(self.flattening())
    }
}

/// Fallible accessor for the mean radius of a body.
pub trait TryMeanRadius: Origin {
    /// Returns the mean radius, or an error if undefined.
    fn try_mean_radius(&self) -> Result<Distance, UndefinedOriginPropertyError>;
}

/// Infallible accessor for the mean radius of a body.
pub trait MeanRadius: Origin {
    /// Returns the mean radius.
    fn mean_radius(&self) -> Distance;
}

impl<T: MeanRadius> TryMeanRadius for T {
    fn try_mean_radius(&self) -> Result<Distance, UndefinedOriginPropertyError> {
        Ok(self.mean_radius())
    }
}

/// Infallible accessor for the gravitational parameter of a body.
pub trait PointMass: Origin {
    /// Returns the gravitational parameter.
    fn gravitational_parameter(&self) -> GravitationalParameter;
}

/// Fallible accessor for the gravitational parameter of a body.
pub trait TryPointMass: Origin {
    /// Returns the gravitational parameter, or an error if undefined.
    fn try_gravitational_parameter(
        &self,
    ) -> Result<GravitationalParameter, UndefinedOriginPropertyError>;
}

impl<T: PointMass> TryPointMass for T {
    fn try_gravitational_parameter(
        &self,
    ) -> Result<GravitationalParameter, UndefinedOriginPropertyError> {
        Ok(self.gravitational_parameter())
    }
}

#[derive(Clone, Copy, Debug)]
pub(crate) enum RotationalElementType {
    RightAscension,
    Declination,
    Rotation,
}

impl RotationalElementType {
    fn sincos(&self, val: f64) -> f64 {
        match self {
            RotationalElementType::Declination => val.cos(),
            _ => val.sin(),
        }
    }

    fn sincos_dot(&self, val: f64) -> f64 {
        match self {
            RotationalElementType::Declination => val.sin(),
            _ => val.cos(),
        }
    }

    fn sign(&self) -> f64 {
        match self {
            RotationalElementType::Declination => -1.0,
            _ => 1.0,
        }
    }

    fn dt(&self) -> f64 {
        match self {
            RotationalElementType::Rotation => SECONDS_PER_DAY,
            _ => SECONDS_PER_JULIAN_CENTURY,
        }
    }
}

pub(crate) struct RotationalElement<const N: usize> {
    typ: RotationalElementType,
    c0: f64,
    c1: f64,
    c2: f64,
    c: [f64; N],
    theta0: [f64; N],
    theta1: [f64; N],
}

impl<const N: usize> RotationalElement<N> {
    fn trig_term(&self, t: f64) -> f64 {
        self.c
            .iter()
            .zip(self.theta0.iter())
            .zip(self.theta1.iter())
            .map(|((&c, &theta0), &theta1)| {
                c * self
                    .typ
                    .sincos(theta0 + theta1 * t / SECONDS_PER_JULIAN_CENTURY)
            })
            .sum()
    }

    fn trig_term_dot(&self, t: f64) -> f64 {
        self.c
            .iter()
            .zip(self.theta0.iter())
            .zip(self.theta1.iter())
            .map(|((&c, &theta0), &theta1)| {
                c * theta1 / SECONDS_PER_JULIAN_CENTURY
                    * self
                        .typ
                        .sincos_dot(theta0 + theta1 * t / SECONDS_PER_JULIAN_CENTURY)
            })
            .sum()
    }

    fn angle(&self, t: f64) -> f64 {
        self.c0
            + self.c1 * t / self.typ.dt()
            + self.c2 * t.powi(2) / self.typ.dt().powi(2)
            + self.trig_term(t)
    }

    fn angle_dot(&self, t: f64) -> f64 {
        self.c1 / self.typ.dt()
            + 2.0 * self.c2 * t / self.typ.dt().powi(2)
            + self.typ.sign() * self.trig_term_dot(t)
    }
}

/// Right ascension, declination, and prime meridian angles in radians.
pub type Elements = (f64, f64, f64);

/// Infallible accessor for the rotational elements of a body.
pub trait RotationalElements: Origin {
    /// Returns the right ascension, declination, and prime meridian at epoch `t` (seconds since J2000 TDB).
    fn rotational_elements(&self, t: f64) -> Elements;

    /// Returns the rates of the rotational elements at epoch `t`.
    fn rotational_element_rates(&self, t: f64) -> Elements;

    /// Returns the right ascension at epoch `t`.
    fn right_ascension(&self, t: f64) -> f64 {
        self.rotational_elements(t).0
    }

    /// Returns the right ascension rate at epoch `t`.
    fn right_ascension_rate(&self, t: f64) -> f64 {
        self.rotational_element_rates(t).0
    }

    /// Returns the declination at epoch `t`.
    fn declination(&self, t: f64) -> f64 {
        self.rotational_elements(t).1
    }

    /// Returns the declination rate at epoch `t`.
    fn declination_rate(&self, t: f64) -> f64 {
        self.rotational_element_rates(t).1
    }

    /// Returns the prime meridian angle at epoch `t`.
    fn rotation_angle(&self, t: f64) -> f64 {
        self.rotational_elements(t).2
    }

    /// Returns the prime meridian rotation rate at epoch `t`.
    fn rotation_rate(&self, t: f64) -> f64 {
        self.rotational_element_rates(t).2
    }
}

/// Fallible accessor for the rotational elements of a body.
pub trait TryRotationalElements: Origin {
    /// Returns the rotational elements at epoch `t`, or an error if undefined.
    fn try_rotational_elements(&self, t: f64) -> Result<Elements, UndefinedOriginPropertyError>;

    /// Returns the rotational element rates at epoch `t`, or an error if undefined.
    fn try_rotational_element_rates(
        &self,
        t: f64,
    ) -> Result<Elements, UndefinedOriginPropertyError>;

    /// Returns the right ascension at epoch `t`, or an error if undefined.
    fn try_right_ascension(&self, t: f64) -> Result<f64, UndefinedOriginPropertyError> {
        self.try_rotational_elements(t).map(|r| r.0)
    }

    /// Returns the right ascension rate at epoch `t`, or an error if undefined.
    fn try_right_ascension_rate(&self, t: f64) -> Result<f64, UndefinedOriginPropertyError> {
        self.try_rotational_element_rates(t).map(|r| r.0)
    }

    /// Returns the declination at epoch `t`, or an error if undefined.
    fn try_declination(&self, t: f64) -> Result<f64, UndefinedOriginPropertyError> {
        self.try_rotational_elements(t).map(|r| r.1)
    }

    /// Returns the declination rate at epoch `t`, or an error if undefined.
    fn try_declination_rate(&self, t: f64) -> Result<f64, UndefinedOriginPropertyError> {
        self.try_rotational_element_rates(t).map(|r| r.1)
    }

    /// Returns the prime meridian angle at epoch `t`, or an error if undefined.
    fn try_rotation_angle(&self, t: f64) -> Result<f64, UndefinedOriginPropertyError> {
        self.try_rotational_elements(t).map(|r| r.2)
    }

    /// Returns the prime meridian rotation rate at epoch `t`, or an error if undefined.
    fn try_rotation_rate(&self, t: f64) -> Result<f64, UndefinedOriginPropertyError> {
        self.try_rotational_element_rates(t).map(|r| r.2)
    }
}

impl<T: RotationalElements> TryRotationalElements for T {
    fn try_rotational_elements(&self, t: f64) -> Result<Elements, UndefinedOriginPropertyError> {
        Ok(self.rotational_elements(t))
    }

    fn try_rotational_element_rates(
        &self,
        t: f64,
    ) -> Result<Elements, UndefinedOriginPropertyError> {
        Ok(self.rotational_element_rates(t))
    }
}

/// Fallible accessor for the J2 zonal harmonic coefficient.
pub trait TryJ2 {
    /// Returns the J2 coefficient, or an error if undefined.
    fn try_j2(&self) -> Result<f64, UndefinedOriginPropertyError>;
}

/// Infallible accessor for the J2 zonal harmonic coefficient.
pub trait J2 {
    /// Returns the J2 coefficient.
    fn j2(&self) -> f64;
}

impl<T> TryJ2 for T
where
    T: J2,
{
    fn try_j2(&self) -> Result<f64, UndefinedOriginPropertyError> {
        Ok(self.j2())
    }
}

impl J2 for Earth {
    fn j2(&self) -> f64 {
        1.08262668e-3
    }
}

#[cfg(test)]
mod tests {
    use lox_test_utils::assert_approx_eq;

    use super::*;

    // Jupiter is manually redefined here using known data. This avoids a dependency on the
    // correctness of the PCK parser to test RotationalElements, and prevents compiler errors
    // when generated files are malformed or deleted in preparation for regeneration.
    #[derive(Debug, Copy, Clone, Eq, PartialEq)]
    pub struct Jupiter;

    impl Origin for Jupiter {
        fn id(&self) -> NaifId {
            NaifId(599)
        }

        fn name(&self) -> &'static str {
            "Jupiter"
        }
    }
    #[derive(Debug, Copy, Clone, Eq, PartialEq)]
    pub struct Rupert;

    impl Origin for Rupert {
        fn id(&self) -> NaifId {
            NaifId(1099)
        }

        fn name(&self) -> &'static str {
            "Persephone/Rupert"
        }
    }

    #[test]
    fn test_body() {
        let body = Jupiter;
        let id = body.id().0;
        let name = body.name();
        assert_eq!(id, 599);
        assert_eq!(name, "Jupiter");

        let body = Rupert;
        let id = body.id().0;
        let name = body.name();
        assert_eq!(id, 1099);
        assert_eq!(name, "Persephone/Rupert");
    }

    const RIGHT_ASCENSION_JUPITER: RotationalElement<15> = RotationalElement {
        typ: RotationalElementType::RightAscension,
        c0: 4.6784701644349695f64,
        c1: -0.00011342894808711148f64,
        c2: 0f64,
        c: [
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0.0000020420352248333656f64,
            0.000016371188383706813f64,
            0.000024993114888558796f64,
            0.0000005235987755982989f64,
            0.00003752457891787809f64,
        ],
        theta0: [
            1.2796754075622423f64,
            0.42970006184100396f64,
            4.9549897464119015f64,
            6.2098814785958245f64,
            2.092649773141201f64,
            4.010766621082969f64,
            6.147922290150026f64,
            1.9783307071355725f64,
            2.5593508151244846f64,
            0.8594001236820079f64,
            1.734171606432425f64,
            3.0699533280603655f64,
            5.241627996900319f64,
            1.9898901100379935f64,
            0.864134346731335f64,
        ],
        theta1: [
            1596.503281347521f64,
            787.7927551311844f64,
            84.66068602648895f64,
            20.792107379008446f64,
            4.574507969477138f64,
            1.1222467090323538f64,
            41.58421475801689f64,
            105.9414855960558f64,
            3193.006562695042f64,
            1575.5855102623689f64,
            84.65553032387855f64,
            20.80363527871787f64,
            4.582318317879813f64,
            105.94580703128374f64,
            1.1222467090323538f64,
        ],
    };

    const DECLINATION_JUPITER: RotationalElement<15> = RotationalElement {
        typ: RotationalElementType::Declination,
        c0: 1.1256553894213766f64,
        c1: 0.00004211479485062318f64,
        c2: 0f64,
        c: [
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0f64,
            0.0000008726646259971648f64,
            0.000007051130178057092f64,
            0.000010768681484805013f64,
            -0.00000022689280275926283f64,
            0.00001616174887346749f64,
        ],
        theta0: [
            1.2796754075622423f64,
            0.42970006184100396f64,
            4.9549897464119015f64,
            6.2098814785958245f64,
            2.092649773141201f64,
            4.010766621082969f64,
            6.147922290150026f64,
            1.9783307071355725f64,
            2.5593508151244846f64,
            0.8594001236820079f64,
            1.734171606432425f64,
            3.0699533280603655f64,
            5.241627996900319f64,
            1.9898901100379935f64,
            0.864134346731335f64,
        ],
        theta1: [
            1596.503281347521f64,
            787.7927551311844f64,
            84.66068602648895f64,
            20.792107379008446f64,
            4.574507969477138f64,
            1.1222467090323538f64,
            41.58421475801689f64,
            105.9414855960558f64,
            3193.006562695042f64,
            1575.5855102623689f64,
            84.65553032387855f64,
            20.80363527871787f64,
            4.582318317879813f64,
            105.94580703128374f64,
            1.1222467090323538f64,
        ],
    };

    const ROTATION_JUPITER: RotationalElement<0> = RotationalElement {
        typ: RotationalElementType::Rotation,
        c0: 4.973315703557842f64,
        c1: 15.193719457141356f64,
        c2: 0f64,
        c: [],
        theta0: [],
        theta1: [],
    };

    impl RotationalElements for Jupiter {
        fn rotational_elements(&self, t: f64) -> Elements {
            (
                RIGHT_ASCENSION_JUPITER.angle(t),
                DECLINATION_JUPITER.angle(t),
                ROTATION_JUPITER.angle(t),
            )
        }

        fn rotational_element_rates(&self, t: f64) -> Elements {
            (
                RIGHT_ASCENSION_JUPITER.angle_dot(t),
                DECLINATION_JUPITER.angle_dot(t),
                ROTATION_JUPITER.angle_dot(t),
            )
        }
    }

    #[test]
    fn test_rotational_elements_right_ascension() {
        assert_approx_eq!(
            Jupiter.right_ascension(0.0),
            4.678480799964803,
            rtol <= 1e-8
        );
    }

    #[test]
    fn test_rotational_elements_right_ascension_dot() {
        assert_approx_eq!(
            Jupiter.right_ascension_rate(0.0),
            -1.3266588500099516e-13,
            rtol <= 1e-8
        );
    }

    #[test]
    fn test_rotational_elements_declination() {
        assert_approx_eq!(Jupiter.declination(0.0), 1.1256642372977634, rtol <= 1e-8);
    }

    #[test]
    fn test_rotational_elements_declination_dot() {
        assert_approx_eq!(
            Jupiter.declination_rate(0.0),
            3.004482367136341e-15,
            rtol <= 1e-8
        );
    }

    #[test]
    fn test_rotational_elements_prime_meridian() {
        assert_approx_eq!(Jupiter.rotation_angle(0.0), 4.973315703557842, rtol <= 1e-8);
    }

    #[test]
    fn test_rotational_elements_prime_meridian_dot() {
        assert_approx_eq!(
            Jupiter.rotation_rate(0.0),
            0.00017585323445765458,
            rtol <= 1e-8
        );
    }

    #[test]
    fn test_rotational_elements_error() {
        assert!(DynOrigin::Mundilfari.try_rotational_elements(0.0).is_err());
        assert!(
            DynOrigin::Mundilfari
                .try_rotational_element_rates(0.0)
                .is_err()
        );
    }
}