#[cfg(feature = "de440")]
mod de440_backend;
mod runtime_backend;
mod vsop87_backend;
#[cfg(feature = "de440")]
pub use de440_backend::De440Ephemeris;
pub use runtime_backend::RuntimeEphemeris;
pub use vsop87_backend::Vsop87Ephemeris;
use crate::coordinates::{
cartesian::{Position, Velocity},
centers::{Barycentric, Geocentric, Heliocentric},
frames::EclipticMeanJ2000,
};
use crate::time::JulianDate;
use qtty::{AstronomicalUnit, Day, Kilometer};
pub type AuPerDay = qtty::Per<AstronomicalUnit, Day>;
pub trait Ephemeris {
fn sun_barycentric(
jd: JulianDate,
) -> Position<Barycentric, EclipticMeanJ2000, AstronomicalUnit>;
fn earth_barycentric(
jd: JulianDate,
) -> Position<Barycentric, EclipticMeanJ2000, AstronomicalUnit>;
fn earth_heliocentric(
jd: JulianDate,
) -> Position<Heliocentric, EclipticMeanJ2000, AstronomicalUnit>;
fn earth_barycentric_velocity(jd: JulianDate) -> Velocity<EclipticMeanJ2000, AuPerDay>;
fn moon_geocentric(jd: JulianDate) -> Position<Geocentric, EclipticMeanJ2000, Kilometer>;
}
pub trait DynEphemeris {
fn sun_barycentric(
&self,
jd: JulianDate,
) -> Position<Barycentric, EclipticMeanJ2000, AstronomicalUnit>;
fn earth_barycentric(
&self,
jd: JulianDate,
) -> Position<Barycentric, EclipticMeanJ2000, AstronomicalUnit>;
fn earth_heliocentric(
&self,
jd: JulianDate,
) -> Position<Heliocentric, EclipticMeanJ2000, AstronomicalUnit>;
fn earth_barycentric_velocity(&self, jd: JulianDate) -> Velocity<EclipticMeanJ2000, AuPerDay>;
fn moon_geocentric(&self, jd: JulianDate)
-> Position<Geocentric, EclipticMeanJ2000, Kilometer>;
}
impl<T: Ephemeris> DynEphemeris for T {
#[inline]
fn sun_barycentric(
&self,
jd: JulianDate,
) -> Position<Barycentric, EclipticMeanJ2000, AstronomicalUnit> {
<T as Ephemeris>::sun_barycentric(jd)
}
#[inline]
fn earth_barycentric(
&self,
jd: JulianDate,
) -> Position<Barycentric, EclipticMeanJ2000, AstronomicalUnit> {
<T as Ephemeris>::earth_barycentric(jd)
}
#[inline]
fn earth_heliocentric(
&self,
jd: JulianDate,
) -> Position<Heliocentric, EclipticMeanJ2000, AstronomicalUnit> {
<T as Ephemeris>::earth_heliocentric(jd)
}
#[inline]
fn earth_barycentric_velocity(&self, jd: JulianDate) -> Velocity<EclipticMeanJ2000, AuPerDay> {
<T as Ephemeris>::earth_barycentric_velocity(jd)
}
#[inline]
fn moon_geocentric(
&self,
jd: JulianDate,
) -> Position<Geocentric, EclipticMeanJ2000, Kilometer> {
<T as Ephemeris>::moon_geocentric(jd)
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::calculus::ephemeris::Vsop87Ephemeris;
const JD_J2000: f64 = 2451545.0;
fn jd() -> JulianDate {
JulianDate::new(JD_J2000)
}
#[test]
fn dyn_ephemeris_sun_barycentric_via_blanket_impl() {
let eph: &dyn DynEphemeris = &Vsop87Ephemeris;
let pos = eph.sun_barycentric(jd());
let mag =
(pos.x().value().powi(2) + pos.y().value().powi(2) + pos.z().value().powi(2)).sqrt();
assert!(mag < 0.02, "Sun-SSB offset too large: {mag} AU");
assert!(pos.x().value().is_finite());
}
#[test]
fn dyn_ephemeris_earth_barycentric_via_blanket_impl() {
let eph: &dyn DynEphemeris = &Vsop87Ephemeris;
let pos = eph.earth_barycentric(jd());
let mag =
(pos.x().value().powi(2) + pos.y().value().powi(2) + pos.z().value().powi(2)).sqrt();
assert!(
mag > 0.9 && mag < 1.1,
"Earth-SSB distance should be ~1 AU, got {mag}"
);
}
#[test]
fn dyn_ephemeris_earth_heliocentric_via_blanket_impl() {
let eph: &dyn DynEphemeris = &Vsop87Ephemeris;
let pos = eph.earth_heliocentric(jd());
let mag =
(pos.x().value().powi(2) + pos.y().value().powi(2) + pos.z().value().powi(2)).sqrt();
assert!(
mag > 0.98 && mag < 1.02,
"Earth heliocentric should be ~1 AU, got {mag}"
);
}
#[test]
fn dyn_ephemeris_earth_barycentric_velocity_via_blanket_impl() {
let eph: &dyn DynEphemeris = &Vsop87Ephemeris;
let vel = eph.earth_barycentric_velocity(jd());
let speed =
(vel.x().value().powi(2) + vel.y().value().powi(2) + vel.z().value().powi(2)).sqrt();
assert!(
speed > 0.01 && speed < 0.03,
"Earth speed ~0.017 AU/day, got {speed}"
);
}
#[test]
fn dyn_ephemeris_moon_geocentric_via_blanket_impl() {
let eph: &dyn DynEphemeris = &Vsop87Ephemeris;
let pos = eph.moon_geocentric(jd());
let dist =
(pos.x().value().powi(2) + pos.y().value().powi(2) + pos.z().value().powi(2)).sqrt();
assert!(
dist > 350_000.0 && dist < 420_000.0,
"Moon geocentric dist: {dist} km"
);
}
#[test]
fn static_vs_dynamic_dispatch_agree() {
let jd_val = jd();
let pos_static = <Vsop87Ephemeris as Ephemeris>::sun_barycentric(jd_val);
let pos_dyn = {
let eph: &dyn DynEphemeris = &Vsop87Ephemeris;
eph.sun_barycentric(jd_val)
};
assert!((pos_static.x().value() - pos_dyn.x().value()).abs() < 1e-15);
assert!((pos_static.y().value() - pos_dyn.y().value()).abs() < 1e-15);
assert!((pos_static.z().value() - pos_dyn.z().value()).abs() < 1e-15);
}
}