diff --git a/crates/lox-space/examples/iss.rs b/crates/lox-space/examples/iss.rs index 1e137215..1d28f386 100644 --- a/crates/lox-space/examples/iss.rs +++ b/crates/lox-space/examples/iss.rs @@ -14,12 +14,9 @@ fn main() { let epoch = Epoch::from_date_and_time(TimeScale::TDB, date, time); let position = DVec3::new(6068279.27, -1692843.94, -2516619.18) * 1e-3; let velocity = DVec3::new(-660.415582, 5495.938726, -5303.093233) * 1e-3; - let iss = Cartesian::new(epoch, Earth, position, velocity); + let iss = Cartesian::new(epoch, Earth, Icrf, position, velocity); - println!( - "ISS Orbit for Julian Day {}", - iss.epoch().days_since_j2000(), - ); + println!("ISS Orbit for Julian Day {}", iss.time().days_since_j2000(),); println!("============================="); println!("Semi-major axis: {:.3} km", iss.semi_major()); println!("Eccentricity: {:.6}", iss.eccentricity()); diff --git a/crates/lox-space/src/prelude.rs b/crates/lox-space/src/prelude.rs index 20f7e2b7..89b86e94 100644 --- a/crates/lox-space/src/prelude.rs +++ b/crates/lox-space/src/prelude.rs @@ -7,8 +7,7 @@ */ pub use lox_core::bodies::*; - +pub use lox_core::frames::*; pub use lox_core::time::dates::*; pub use lox_core::time::epochs::*; - pub use lox_core::two_body::{Cartesian, DVec3, Keplerian, TwoBody}; diff --git a/crates/lox_core/src/bodies.rs b/crates/lox_core/src/bodies.rs index 7808175d..58d0c0dd 100644 --- a/crates/lox_core/src/bodies.rs +++ b/crates/lox_core/src/bodies.rs @@ -92,6 +92,28 @@ macro_rules! body { } } + impl Display for $i { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.name()) + } + } + }; + ($i:ident, $t:ident, $name:literal, $naif_id:literal) => { + #[derive(Clone, Copy, Debug, Eq, PartialEq)] + pub struct $i; + + impl $t for $i {} + + impl Body for $i { + fn id(&self) -> NaifId { + NaifId($naif_id) + } + + fn name(&self) -> &'static str { + $name + } + } + impl Display for $i { fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { write!(f, "{}", self.name()) @@ -118,16 +140,19 @@ body! { Neptune, Planet, 899 } body! { Pluto, Planet, 999 } // Barycenters. -body! { SolarSystemBarycenter, "Solar System Barycenter", 0 } -body! { MercuryBarycenter, "Mercury Barycenter", 1 } -body! { VenusBarycenter, "Venus Barycenter", 2 } -body! { EarthBarycenter, "Earth Barycenter", 3 } -body! { MarsBarycenter, "Mars Barycenter", 4 } -body! { JupiterBarycenter, "Jupiter Barycenter", 5 } -body! { SaturnBarycenter, "Saturn Barycenter", 6 } -body! { UranusBarycenter, "Uranus Barycenter", 7 } -body! { NeptuneBarycenter, "Neptune Barycenter", 8 } -body! { PlutoBarycenter, "Pluto Barycenter", 9 } +pub trait Barycenter: PointMass + DynClone {} +clone_trait_object!(Barycenter); + +body! { SolarSystemBarycenter, Barycenter, "Solar System Barycenter", 0 } +body! { MercuryBarycenter, Barycenter, "Mercury Barycenter", 1 } +body! { VenusBarycenter, Barycenter, "Venus Barycenter", 2 } +body! { EarthBarycenter, Barycenter, "Earth Barycenter", 3 } +body! { MarsBarycenter, Barycenter, "Mars Barycenter", 4 } +body! { JupiterBarycenter, Barycenter, "Jupiter Barycenter", 5 } +body! { SaturnBarycenter, Barycenter, "Saturn Barycenter", 6 } +body! { UranusBarycenter, Barycenter, "Uranus Barycenter", 7 } +body! { NeptuneBarycenter, Barycenter, "Neptune Barycenter", 8 } +body! { PlutoBarycenter, Barycenter, "Pluto Barycenter", 9 } impl PointMass for SolarSystemBarycenter { fn gravitational_parameter(&self) -> f64 { @@ -545,12 +570,10 @@ pub trait TriAxial: Ellipsoid { fn along_orbit_radius(&self) -> f64; } -pub trait PointMass: Body + DynClone { +pub trait PointMass: Body { fn gravitational_parameter(&self) -> f64; } -clone_trait_object!(PointMass); - pub type PolynomialCoefficients = (f64, f64, f64, &'static [f64]); pub type NutationPrecessionCoefficients = (&'static [f64], &'static [f64]); diff --git a/crates/lox_core/src/frames.rs b/crates/lox_core/src/frames.rs index 3dc8cd24..bdffaece 100644 --- a/crates/lox_core/src/frames.rs +++ b/crates/lox_core/src/frames.rs @@ -15,11 +15,25 @@ pub mod iau; // TODO: Replace with proper `Epoch` type type Epoch = f64; -pub trait ReferenceFrame { - fn is_inertial(&self) -> bool; +pub trait ReferenceFrame {} + +pub trait InertialFrame: ReferenceFrame { + fn is_inertial(&self) -> bool { + true + } + + fn is_rotating(&self) -> bool { + false + } +} + +pub trait RotatingFrame: ReferenceFrame { + fn is_inertial(&self) -> bool { + false + } fn is_rotating(&self) -> bool { - !self.is_inertial() + true } } @@ -32,11 +46,8 @@ impl Display for Icrf { } } -impl ReferenceFrame for Icrf { - fn is_inertial(&self) -> bool { - true - } -} +impl ReferenceFrame for Icrf {} +impl InertialFrame for Icrf {} pub fn rotation_matrix_derivative(m: DMat3, v: DVec3) -> DMat3 { let sx = DVec3::new(0.0, v.z, v.y); diff --git a/crates/lox_core/src/frames/iau.rs b/crates/lox_core/src/frames/iau.rs index 91d700ae..19c81443 100644 --- a/crates/lox_core/src/frames/iau.rs +++ b/crates/lox_core/src/frames/iau.rs @@ -6,22 +6,25 @@ * file, you can obtain one at https://mozilla.org/MPL/2.0/. */ -use std::fmt::Debug; +use std::fmt::{Debug, Display, Formatter}; use glam::{DMat3, DVec3}; use crate::bodies::RotationalElements; -use crate::frames::{Epoch, FromFrame, Icrf, ReferenceFrame, Rotation}; +use crate::frames::{Epoch, FromFrame, Icrf, ReferenceFrame, RotatingFrame, Rotation}; #[derive(Debug, Copy, Clone)] pub struct BodyFixed(pub T); -impl ReferenceFrame for BodyFixed { - fn is_inertial(&self) -> bool { - false +impl Display for BodyFixed { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + write!(f, "IAU_{}", &self.0.name().to_uppercase()) } } +impl ReferenceFrame for BodyFixed {} +impl RotatingFrame for BodyFixed {} + impl FromFrame for BodyFixed { fn rotation_from(&self, _: Icrf, t: Epoch) -> Rotation { let (right_ascension, declination, prime_meridian) = T::rotational_elements(t); @@ -49,6 +52,7 @@ mod tests { #[test] fn test_bodyfixed() { let iau_jupiter = BodyFixed(Jupiter); + assert_eq!(format!("{}", iau_jupiter), "IAU_JUPITER"); assert!(!iau_jupiter.is_inertial()); assert!(iau_jupiter.is_rotating()); diff --git a/crates/lox_core/src/two_body.rs b/crates/lox_core/src/two_body.rs index 0bb7dfb2..c2471e16 100644 --- a/crates/lox_core/src/two_body.rs +++ b/crates/lox_core/src/two_body.rs @@ -9,25 +9,181 @@ pub use glam::DVec3; use crate::bodies::PointMass; +use crate::frames::{InertialFrame, ReferenceFrame}; use crate::time::epochs::Epoch; use crate::two_body::elements::{cartesian_to_keplerian, keplerian_to_cartesian}; pub mod elements; +#[derive(Debug, Clone, PartialEq)] +pub struct CartesianState { + position: DVec3, + velocity: DVec3, +} + +impl CartesianState { + pub fn new(position: DVec3, velocity: DVec3) -> Self { + Self { position, velocity } + } + + pub fn from_coords(x: f64, y: f64, z: f64, vx: f64, vy: f64, vz: f64) -> Self { + let position = DVec3::new(x, y, z); + let velocity = DVec3::new(vx, vy, vz); + Self::new(position, velocity) + } + + pub fn position(&self) -> DVec3 { + self.position + } + + pub fn velocity(&self) -> DVec3 { + self.velocity + } + + pub fn to_keplerian_state(&self, grav_param: f64) -> KeplerianState { + let (semi_major, eccentricity, inclination, ascending_node, periapsis_arg, true_anomaly) = + cartesian_to_keplerian(grav_param, self.position, self.velocity); + KeplerianState { + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + } + } + + pub fn semi_major(&self, grav_param: f64) -> f64 { + self.to_keplerian_state(grav_param).semi_major + } + + pub fn eccentricity(&self, grav_param: f64) -> f64 { + self.to_keplerian_state(grav_param).eccentricity + } + + pub fn inclination(&self, grav_param: f64) -> f64 { + self.to_keplerian_state(grav_param).inclination + } + + pub fn ascending_node(&self, grav_param: f64) -> f64 { + self.to_keplerian_state(grav_param).ascending_node + } + + pub fn periapsis_arg(&self, grav_param: f64) -> f64 { + self.to_keplerian_state(grav_param).periapsis_arg + } + + pub fn true_anomaly(&self, grav_param: f64) -> f64 { + self.to_keplerian_state(grav_param).true_anomaly + } +} + +#[derive(Debug, Clone, PartialEq)] +pub struct KeplerianState { + semi_major: f64, + eccentricity: f64, + inclination: f64, + ascending_node: f64, + periapsis_arg: f64, + true_anomaly: f64, +} + +impl KeplerianState { + pub fn new( + semi_major: f64, + eccentricity: f64, + inclination: f64, + ascending_node: f64, + periapsis_arg: f64, + true_anomaly: f64, + ) -> Self { + Self { + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + } + } + + pub fn position(&self, grav_param: f64) -> DVec3 { + self.to_cartesian_state(grav_param).position + } + + pub fn velocity(&self, grav_param: f64) -> DVec3 { + self.to_cartesian_state(grav_param).velocity + } + + pub fn to_cartesian_state(&self, grav_param: f64) -> CartesianState { + let (position, velocity) = keplerian_to_cartesian( + grav_param, + self.semi_major, + self.eccentricity, + self.inclination, + self.ascending_node, + self.periapsis_arg, + self.true_anomaly, + ); + CartesianState::new(position, velocity) + } + + pub fn semi_major(&self) -> f64 { + self.semi_major + } + + pub fn eccentricity(&self) -> f64 { + self.eccentricity + } + + pub fn inclination(&self) -> f64 { + self.inclination + } + + pub fn ascending_node(&self) -> f64 { + self.ascending_node + } + + pub fn periapsis_arg(&self) -> f64 { + self.periapsis_arg + } + + pub fn true_anomaly(&self) -> f64 { + self.true_anomaly + } +} + pub type Elements = (f64, f64, f64, f64, f64, f64); -pub trait Center { - type Center; +impl From for KeplerianState { + fn from(elements: Elements) -> Self { + let (semi_major, eccentricity, inclination, ascending_node, periapsis_arg, true_anomaly) = + elements; + Self { + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + } + } +} + +pub trait CoordinateSystem { + type Origin: PointMass; + type Frame: ReferenceFrame; - fn center(&self) -> Self::Center; + fn origin(&self) -> Self::Origin; + fn reference_frame(&self) -> Self::Frame; } pub trait TwoBody { - fn epoch(&self) -> Epoch; + fn time(&self) -> Epoch; + fn to_cartesian_state(&self) -> CartesianState; + fn to_keplerian_state(&self) -> KeplerianState; fn position(&self) -> DVec3; fn velocity(&self) -> DVec3; - fn cartesian(&self) -> (DVec3, DVec3); - fn keplerian(&self) -> Elements; fn semi_major(&self) -> f64; fn eccentricity(&self) -> f64; fn inclination(&self) -> f64; @@ -37,193 +193,239 @@ pub trait TwoBody { } #[derive(Debug, Clone)] -pub struct Cartesian { - epoch: Epoch, - center: T, - position: DVec3, - velocity: DVec3, +pub struct Cartesian +where + T: PointMass + Copy, + S: ReferenceFrame + Copy, +{ + time: Epoch, + origin: T, + frame: S, + state: CartesianState, } -impl Cartesian { - pub fn new(epoch: Epoch, center: T, position: DVec3, velocity: DVec3) -> Self { +impl Cartesian +where + T: PointMass + Copy, + S: ReferenceFrame + Copy, +{ + pub fn new(time: Epoch, origin: T, frame: S, position: DVec3, velocity: DVec3) -> Self { Self { - epoch, - center, - position, - velocity, + time, + origin, + frame, + state: CartesianState { position, velocity }, } } } -impl Center for Cartesian { - type Center = T; +impl CoordinateSystem for Cartesian +where + T: PointMass + Copy, + S: ReferenceFrame + Copy, +{ + type Origin = T; + type Frame = S; - fn center(&self) -> Self::Center { - self.center + fn origin(&self) -> Self::Origin { + self.origin + } + + fn reference_frame(&self) -> Self::Frame { + self.frame } } -impl TwoBody for Cartesian { - fn epoch(&self) -> Epoch { - self.epoch +impl TwoBody for Cartesian +where + T: PointMass + Copy, + S: ReferenceFrame + Copy, +{ + fn time(&self) -> Epoch { + self.time } - fn position(&self) -> DVec3 { - self.position + fn to_cartesian_state(&self) -> CartesianState { + self.state.clone() } - fn velocity(&self) -> DVec3 { - self.velocity + fn to_keplerian_state(&self) -> KeplerianState { + let mu = self.origin.gravitational_parameter(); + self.state.to_keplerian_state(mu) } - fn cartesian(&self) -> (DVec3, DVec3) { - (self.position, self.velocity) + fn position(&self) -> DVec3 { + self.state.position } - fn keplerian(&self) -> Elements { - let mu = self.center.gravitational_parameter(); - cartesian_to_keplerian(mu, self.position, self.velocity) + fn velocity(&self) -> DVec3 { + self.state.velocity } fn semi_major(&self) -> f64 { - self.keplerian().0 + self.to_keplerian_state().semi_major } fn eccentricity(&self) -> f64 { - self.keplerian().1 + self.to_keplerian_state().eccentricity } fn inclination(&self) -> f64 { - self.keplerian().2 + self.to_keplerian_state().inclination } fn ascending_node(&self) -> f64 { - self.keplerian().3 + self.to_keplerian_state().ascending_node } fn periapsis_arg(&self) -> f64 { - self.keplerian().4 + self.to_keplerian_state().periapsis_arg } fn true_anomaly(&self) -> f64 { - self.keplerian().5 + self.to_keplerian_state().true_anomaly } } -impl From> for Cartesian { - fn from(value: Keplerian) -> Self { - let epoch = value.epoch; - let center = value.center; - let (pos, vel) = value.cartesian(); - Cartesian::new(epoch, center, pos, vel) +impl From> for Cartesian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + fn from(keplerian: Keplerian) -> Self { + let time = keplerian.time; + let origin = keplerian.origin; + let frame = keplerian.frame; + let state = keplerian.to_cartesian_state(); + Cartesian { + time, + origin, + frame, + state, + } } } #[derive(Debug, Clone)] -pub struct Keplerian { - epoch: Epoch, - center: T, - semi_major: f64, - eccentricity: f64, - inclination: f64, - ascending_node: f64, - periapsis_arg: f64, - true_anomaly: f64, +pub struct Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + time: Epoch, + origin: T, + frame: S, + state: KeplerianState, } -impl Keplerian { - pub fn new(epoch: Epoch, center: T, elements: Elements) -> Self { - let (semi_major, eccentricity, inclination, ascending_node, periapsis_arg, true_anomaly) = - elements; +impl Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + pub fn new(time: Epoch, origin: T, frame: S, elements: Elements) -> Self { Self { - epoch, - center, - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, + time, + origin, + frame, + state: elements.into(), } } } -impl Center for Keplerian { - type Center = T; +impl CoordinateSystem for Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + type Origin = T; + type Frame = S; - fn center(&self) -> Self::Center { - self.center + fn origin(&self) -> Self::Origin { + self.origin + } + + fn reference_frame(&self) -> Self::Frame { + self.frame } } -impl TwoBody for Keplerian { - fn epoch(&self) -> Epoch { - self.epoch +impl TwoBody for Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + fn time(&self) -> Epoch { + self.time } - fn position(&self) -> DVec3 { - self.cartesian().0 + fn to_cartesian_state(&self) -> CartesianState { + let mu = self.origin.gravitational_parameter(); + let (position, velocity) = keplerian_to_cartesian( + mu, + self.state.semi_major, + self.state.eccentricity, + self.state.inclination, + self.state.ascending_node, + self.state.periapsis_arg, + self.state.true_anomaly, + ); + CartesianState { position, velocity } } - fn velocity(&self) -> DVec3 { - self.cartesian().1 + fn to_keplerian_state(&self) -> KeplerianState { + self.state.clone() } - fn cartesian(&self) -> (DVec3, DVec3) { - let mu = self.center.gravitational_parameter(); - keplerian_to_cartesian( - mu, - self.semi_major, - self.eccentricity, - self.inclination, - self.ascending_node, - self.periapsis_arg, - self.true_anomaly, - ) + fn position(&self) -> DVec3 { + self.to_cartesian_state().position } - fn keplerian(&self) -> Elements { - ( - self.semi_major, - self.eccentricity, - self.inclination, - self.ascending_node, - self.periapsis_arg, - self.true_anomaly, - ) + fn velocity(&self) -> DVec3 { + self.to_cartesian_state().velocity } fn semi_major(&self) -> f64 { - self.semi_major + self.state.semi_major } fn eccentricity(&self) -> f64 { - self.eccentricity + self.state.eccentricity } fn inclination(&self) -> f64 { - self.inclination + self.state.inclination } fn ascending_node(&self) -> f64 { - self.ascending_node + self.state.ascending_node } fn periapsis_arg(&self) -> f64 { - self.periapsis_arg + self.state.periapsis_arg } fn true_anomaly(&self) -> f64 { - self.true_anomaly + self.state.true_anomaly } } -impl From> for Keplerian { - fn from(value: Cartesian) -> Self { - let epoch = value.epoch; - let center = value.center; - let elements = value.keplerian(); - Self::new(epoch, center, elements) +impl From> for Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + fn from(cartesian: Cartesian) -> Self { + let time = cartesian.time; + let origin = cartesian.origin; + let frame = cartesian.frame; + let state = cartesian.to_keplerian_state(); + Self { + time, + origin, + frame, + state, + } } } @@ -234,6 +436,7 @@ mod tests { use float_eq::assert_float_eq; use crate::bodies::Earth; + use crate::frames::Icrf; use crate::time::dates::{Date, Time}; use crate::time::epochs::TimeScale; @@ -263,14 +466,17 @@ mod tests { ) .mul(1e-3); - let cartesian = Cartesian::new(epoch, Earth, pos, vel); + let cartesian = Cartesian::new(epoch, Earth, Icrf, pos, vel); assert_eq!( - cartesian.cartesian(), - (cartesian.position(), cartesian.velocity()) + cartesian.to_cartesian_state(), + CartesianState { + position: cartesian.position(), + velocity: cartesian.velocity(), + } ); - assert_eq!(cartesian.epoch(), epoch); - assert_eq!(cartesian.center(), Earth); + assert_eq!(cartesian.time(), epoch); + assert_eq!(cartesian.origin(), Earth); assert_eq!(cartesian.position(), pos); assert_eq!(cartesian.velocity(), vel); assert_float_eq!(cartesian.semi_major(), semi_major, rel <= 1e-6); @@ -283,6 +489,7 @@ mod tests { let keplerian = Keplerian::new( epoch, Earth, + Icrf, ( semi_major, eccentricity, @@ -294,18 +501,18 @@ mod tests { ); assert_eq!( - keplerian.keplerian(), - ( + keplerian.to_keplerian_state(), + KeplerianState { semi_major, eccentricity, inclination, ascending_node, periapsis_arg, true_anomaly - ) + } ); - assert_eq!(keplerian.epoch(), epoch); - assert_eq!(keplerian.center(), Earth); + assert_eq!(keplerian.time(), epoch); + assert_eq!(keplerian.origin(), Earth); assert_float_eq!(keplerian.position().x, pos.x, rel <= 1e-8); assert_float_eq!(keplerian.position().y, pos.y, rel <= 1e-8); assert_float_eq!(keplerian.position().z, pos.z, rel <= 1e-8); @@ -322,26 +529,66 @@ mod tests { let cartesian1 = Cartesian::from(keplerian.clone()); let keplerian1 = Keplerian::from(cartesian.clone()); - assert_float_eq!(cartesian.position.x, cartesian1.position.x, rel <= 1e-8); - assert_float_eq!(cartesian.position.y, cartesian1.position.y, rel <= 1e-8); - assert_float_eq!(cartesian.position.z, cartesian1.position.z, rel <= 1e-8); - assert_float_eq!(cartesian.velocity.x, cartesian1.velocity.x, rel <= 1e-6); - assert_float_eq!(cartesian.velocity.y, cartesian1.velocity.y, rel <= 1e-6); - assert_float_eq!(cartesian.velocity.z, cartesian1.velocity.z, rel <= 1e-6); + assert_float_eq!( + cartesian.state.position.x, + cartesian1.state.position.x, + rel <= 1e-8 + ); + assert_float_eq!( + cartesian.state.position.y, + cartesian1.state.position.y, + rel <= 1e-8 + ); + assert_float_eq!( + cartesian.state.position.z, + cartesian1.state.position.z, + rel <= 1e-8 + ); + assert_float_eq!( + cartesian.state.velocity.x, + cartesian1.state.velocity.x, + rel <= 1e-6 + ); + assert_float_eq!( + cartesian.state.velocity.y, + cartesian1.state.velocity.y, + rel <= 1e-6 + ); + assert_float_eq!( + cartesian.state.velocity.z, + cartesian1.state.velocity.z, + rel <= 1e-6 + ); - assert_float_eq!(keplerian.semi_major, keplerian1.semi_major, rel <= 1e-2); - assert_float_eq!(keplerian.eccentricity, keplerian1.eccentricity, abs <= 1e-6); - assert_float_eq!(keplerian.inclination, keplerian1.inclination, rel <= 1e-6); assert_float_eq!( - keplerian.ascending_node, - keplerian1.ascending_node, + keplerian.state.semi_major, + keplerian1.state.semi_major, + rel <= 1e-2 + ); + assert_float_eq!( + keplerian.state.eccentricity, + keplerian1.state.eccentricity, + abs <= 1e-6 + ); + assert_float_eq!( + keplerian.state.inclination, + keplerian1.state.inclination, + rel <= 1e-6 + ); + assert_float_eq!( + keplerian.state.ascending_node, + keplerian1.state.ascending_node, + rel <= 1e-6 + ); + assert_float_eq!( + keplerian.state.periapsis_arg, + keplerian1.state.periapsis_arg, rel <= 1e-6 ); assert_float_eq!( - keplerian.periapsis_arg, - keplerian1.periapsis_arg, + keplerian.state.true_anomaly, + keplerian1.state.true_anomaly, rel <= 1e-6 ); - assert_float_eq!(keplerian.true_anomaly, keplerian1.true_anomaly, rel <= 1e-6); } } diff --git a/crates/lox_py/src/bodies.rs b/crates/lox_py/src/bodies.rs index 0a423b7d..8998c2d6 100644 --- a/crates/lox_py/src/bodies.rs +++ b/crates/lox_py/src/bodies.rs @@ -6,6 +6,7 @@ * file, you can obtain one at https://mozilla.org/MPL/2.0/. */ +use pyo3::exceptions::PyValueError; use pyo3::prelude::*; use lox_core::bodies::*; @@ -13,7 +14,7 @@ use lox_core::bodies::*; use crate::LoxPyError; #[pyclass(name = "Sun")] -#[derive(Clone)] +#[derive(Clone, Default)] pub struct PySun; #[pymethods] @@ -58,13 +59,13 @@ impl PySun { #[pyclass(name = "Barycenter")] #[derive(Clone)] -pub struct PyBarycenter(Box); +pub struct PyBarycenter(Box); #[pymethods] impl PyBarycenter { #[new] pub fn new(name: &str) -> Result { - let barycenter: Option> = match name { + let barycenter: Option> = match name { "ssb" | "SSB" | "solar system barycenter" | "Solar System Barycenter" => { Some(Box::new(SolarSystemBarycenter)) } @@ -318,6 +319,7 @@ impl PyMinorBody { } } +#[derive(Clone)] pub enum PyBody { Barycenter(PyBarycenter), Sun(PySun), @@ -325,3 +327,71 @@ pub enum PyBody { Satellite(PySatellite), MinorBody(PyMinorBody), } + +impl From for PyObject { + fn from(body: PyBody) -> Self { + Python::with_gil(|py| match body { + PyBody::Barycenter(barycenter) => barycenter.clone().into_py(py), + PyBody::Sun(sun) => sun.clone().into_py(py), + PyBody::Planet(planet) => planet.clone().into_py(py), + PyBody::Satellite(satellite) => satellite.clone().into_py(py), + PyBody::MinorBody(minor_body) => minor_body.clone().into_py(py), + }) + } +} + +impl TryFrom for PyBody { + type Error = PyErr; + + fn try_from(body: PyObject) -> Result { + Python::with_gil(|py| { + if let Ok(body) = body.extract::(py) { + Ok(PyBody::Barycenter(body)) + } else if let Ok(body) = body.extract::(py) { + Ok(PyBody::Sun(body)) + } else if let Ok(body) = body.extract::(py) { + Ok(PyBody::Planet(body)) + } else if let Ok(body) = body.extract::(py) { + Ok(PyBody::Satellite(body)) + } else if let Ok(body) = body.extract::(py) { + Ok(PyBody::MinorBody(body)) + } else { + Err(PyValueError::new_err("Invalid body")) + } + }) + } +} + +impl Body for PyBody { + fn id(&self) -> NaifId { + match &self { + PyBody::Barycenter(barycenter) => NaifId(barycenter.id()), + PyBody::Sun(sun) => NaifId(sun.id()), + PyBody::Planet(planet) => NaifId(planet.id()), + PyBody::Satellite(satellite) => NaifId(satellite.id()), + PyBody::MinorBody(minor_body) => NaifId(minor_body.id()), + } + } + + fn name(&self) -> &'static str { + match &self { + PyBody::Barycenter(barycenter) => barycenter.name(), + PyBody::Sun(sun) => sun.name(), + PyBody::Planet(planet) => planet.name(), + PyBody::Satellite(satellite) => satellite.name(), + PyBody::MinorBody(minor_body) => minor_body.name(), + } + } +} + +impl PointMass for PyBody { + fn gravitational_parameter(&self) -> f64 { + match &self { + PyBody::Barycenter(barycenter) => barycenter.gravitational_parameter(), + PyBody::Sun(sun) => sun.gravitational_parameter(), + PyBody::Planet(planet) => planet.gravitational_parameter(), + PyBody::Satellite(satellite) => satellite.gravitational_parameter(), + PyBody::MinorBody(minor_body) => minor_body.gravitational_parameter(), + } + } +} diff --git a/crates/lox_py/src/frames.rs b/crates/lox_py/src/frames.rs new file mode 100644 index 00000000..22702b14 --- /dev/null +++ b/crates/lox_py/src/frames.rs @@ -0,0 +1,68 @@ +/* + * Copyright (c) 2024. Helge Eichhorn and the LOX contributors + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, you can obtain one at https://mozilla.org/MPL/2.0/. + */ + +use std::fmt::{Display, Formatter}; +use std::str::FromStr; + +use lox_core::bodies::{Earth, Jupiter, Mars, Mercury, Neptune, Pluto, Saturn, Uranus, Venus}; +use lox_core::frames::iau::BodyFixed; +use lox_core::frames::Icrf; + +use crate::LoxPyError; + +// TODO: Add other supported IAU frames +#[derive(Debug, Clone)] +pub enum PyFrame { + Icrf(Icrf), + IauMercury(BodyFixed), + IauVenus(BodyFixed), + IauEarth(BodyFixed), + IauMars(BodyFixed), + IauJupiter(BodyFixed), + IauSaturn(BodyFixed), + IauUranus(BodyFixed), + IauNeptune(BodyFixed), + IauPluto(BodyFixed), +} + +impl Display for PyFrame { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + match &self { + PyFrame::Icrf(frame) => write!(f, "{}", frame), + PyFrame::IauMercury(frame) => write!(f, "{}", frame), + PyFrame::IauVenus(frame) => write!(f, "{}", frame), + PyFrame::IauEarth(frame) => write!(f, "{}", frame), + PyFrame::IauMars(frame) => write!(f, "{}", frame), + PyFrame::IauJupiter(frame) => write!(f, "{}", frame), + PyFrame::IauSaturn(frame) => write!(f, "{}", frame), + PyFrame::IauUranus(frame) => write!(f, "{}", frame), + PyFrame::IauNeptune(frame) => write!(f, "{}", frame), + PyFrame::IauPluto(frame) => write!(f, "{}", frame), + } + } +} + +impl FromStr for PyFrame { + type Err = LoxPyError; + + fn from_str(name: &str) -> Result { + match name { + "icrf" | "ICRF" => Ok(PyFrame::Icrf(Icrf)), + "iau_mercury" | "IAU_MERCURY" => Ok(PyFrame::IauMercury(BodyFixed(Mercury))), + "iau_venus" | "IAU_VENUS" => Ok(PyFrame::IauVenus(BodyFixed(Venus))), + "iau_earth" | "IAU_EARTH" => Ok(PyFrame::IauEarth(BodyFixed(Earth))), + "iau_mars" | "IAU_MARS" => Ok(PyFrame::IauMars(BodyFixed(Mars))), + "iau_jupiter" | "IAU_JUPITER" => Ok(PyFrame::IauJupiter(BodyFixed(Jupiter))), + "iau_saturn" | "IAU_SATURN" => Ok(PyFrame::IauSaturn(BodyFixed(Saturn))), + "iau_uranus" | "IAU_URANUS" => Ok(PyFrame::IauUranus(BodyFixed(Uranus))), + "iau_neptune" | "IAU_NEPTUNE" => Ok(PyFrame::IauNeptune(BodyFixed(Neptune))), + "iau_pluto" | "IAU_PLUTO" => Ok(PyFrame::IauPluto(BodyFixed(Pluto))), + _ => Err(LoxPyError::InvalidFrame(name.to_string())), + } + } +} diff --git a/crates/lox_py/src/lib.rs b/crates/lox_py/src/lib.rs index 96da2d21..6a3641b1 100644 --- a/crates/lox_py/src/lib.rs +++ b/crates/lox_py/src/lib.rs @@ -12,19 +12,22 @@ use thiserror::Error; use crate::bodies::{PyBarycenter, PyMinorBody, PyPlanet, PySatellite, PySun}; use crate::time::{PyEpoch, PyTimeScale}; -use crate::twobody::PyCartesian; +use crate::two_body::PyCartesian; use lox_core::errors::LoxError; mod bodies; +mod frames; mod time; -mod twobody; +mod two_body; #[derive(Error, Debug)] pub enum LoxPyError { - #[error("invalid time scale `{0}`")] + #[error("unknown time scale `{0}`")] InvalidTimeScale(String), #[error("unknown body `{0}`")] InvalidBody(String), + #[error("unknown frame `{0}`")] + InvalidFrame(String), #[error(transparent)] LoxError(#[from] LoxError), #[error(transparent)] @@ -34,8 +37,9 @@ pub enum LoxPyError { impl From for PyErr { fn from(value: LoxPyError) -> Self { match value { - LoxPyError::InvalidTimeScale(_) => PyValueError::new_err(value.to_string()), - LoxPyError::InvalidBody(_) => PyValueError::new_err(value.to_string()), + LoxPyError::InvalidTimeScale(_) + | LoxPyError::InvalidFrame(_) + | LoxPyError::InvalidBody(_) => PyValueError::new_err(value.to_string()), LoxPyError::LoxError(value) => PyValueError::new_err(value.to_string()), LoxPyError::PyError(value) => value, } diff --git a/crates/lox_py/src/two_body.rs b/crates/lox_py/src/two_body.rs new file mode 100644 index 00000000..f7c17f45 --- /dev/null +++ b/crates/lox_py/src/two_body.rs @@ -0,0 +1,258 @@ +/* + * Copyright (c) 2024. Helge Eichhorn and the LOX contributors + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, you can obtain one at https://mozilla.org/MPL/2.0/. + */ + +use numpy::{pyarray, PyArray1}; +use pyo3::prelude::*; +use std::str::FromStr; + +use lox_core::bodies::PointMass; +use lox_core::time::epochs::Epoch; +use lox_core::two_body::{CartesianState, KeplerianState}; + +use crate::bodies::PyBody; +use crate::frames::PyFrame; +use crate::time::PyEpoch; + +#[pyclass(name = "Cartesian")] +pub struct PyCartesian { + time: Epoch, + origin: PyBody, + frame: PyFrame, + state: CartesianState, +} + +#[pymethods] +impl PyCartesian { + #[new] + #[allow(clippy::too_many_arguments)] + fn new( + t: &PyEpoch, + body: PyObject, + frame: &str, + x: f64, + y: f64, + z: f64, + vx: f64, + vy: f64, + vz: f64, + ) -> PyResult { + let origin: PyBody = body.try_into()?; + let frame = PyFrame::from_str(frame)?; + Ok(Self { + time: t.0, + origin, + frame, + state: CartesianState::from_coords(x, y, z, vx, vy, vz), + }) + } + + fn to_keplerian(&self) -> PyKeplerian { + let mu = self.origin.gravitational_parameter(); + let state = self.state.to_keplerian_state(mu); + PyKeplerian { + time: self.time, + origin: self.origin.clone(), + frame: self.frame.clone(), + state, + } + } + + fn time(&self) -> PyEpoch { + PyEpoch(self.time) + } + + fn reference_frame(&self) -> String { + format!("{}", self.frame) + } + + fn origin(&self) -> PyObject { + self.origin.clone().into() + } + + fn position(&self) -> Py> { + let position = self.state.position(); + Python::with_gil(|py| pyarray![py, position.x, position.y, position.z].into_py(py)) + } + + fn velocity(&self) -> Py> { + let velocity = self.state.velocity(); + Python::with_gil(|py| pyarray![py, velocity.x, velocity.y, velocity.z].into_py(py)) + } + + fn cartesian(&self) -> (Py>, Py>) { + (self.position(), self.velocity()) + } + + fn keplerian(&self) -> Py> { + let mu = self.origin.gravitational_parameter(); + let keplerian = self.state.to_keplerian_state(mu); + Python::with_gil(|py| { + pyarray![ + py, + keplerian.true_anomaly(), + keplerian.eccentricity(), + keplerian.inclination(), + keplerian.ascending_node(), + keplerian.periapsis_arg(), + keplerian.true_anomaly() + ] + .into_py(py) + }) + } + + fn semi_major(&self) -> f64 { + let mu = self.origin.gravitational_parameter(); + self.state.semi_major(mu) + } + + fn eccentricity(&self) -> f64 { + let mu = self.origin.gravitational_parameter(); + self.state.eccentricity(mu) + } + + fn inclination(&self) -> f64 { + let mu = self.origin.gravitational_parameter(); + self.state.inclination(mu) + } + + fn ascending_node(&self) -> f64 { + let mu = self.origin.gravitational_parameter(); + self.state.ascending_node(mu) + } + + fn periapsis_arg(&self) -> f64 { + let mu = self.origin.gravitational_parameter(); + self.state.periapsis_arg(mu) + } + + fn true_anomaly(&self) -> f64 { + let mu = self.origin.gravitational_parameter(); + self.state.true_anomaly(mu) + } +} + +#[pyclass(name = "Keplerian")] +pub struct PyKeplerian { + time: Epoch, + origin: PyBody, + frame: PyFrame, + state: KeplerianState, +} + +#[pymethods] +impl PyKeplerian { + #[new] + #[allow(clippy::too_many_arguments)] + fn new( + t: &PyEpoch, + body: PyObject, + frame: &str, + semi_major: f64, + eccentricity: f64, + inclination: f64, + ascending_node: f64, + periapsis_arg: f64, + true_anomaly: f64, + ) -> PyResult { + let origin: PyBody = body.try_into()?; + let frame = PyFrame::from_str(frame)?; + Ok(Self { + time: t.0, + origin, + frame, + state: KeplerianState::new( + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ), + }) + } + + fn to_cartesian(&self) -> PyCartesian { + let mu = self.origin.gravitational_parameter(); + let state = self.state.to_cartesian_state(mu); + PyCartesian { + time: self.time, + origin: self.origin.clone(), + frame: self.frame.clone(), + state, + } + } + + fn time(&self) -> PyEpoch { + PyEpoch(self.time) + } + + fn reference_frame(&self) -> String { + format!("{}", self.frame) + } + + fn origin(&self) -> PyObject { + self.origin.clone().into() + } + + fn position(&self) -> Py> { + let mu = self.origin.gravitational_parameter(); + let cartesian = self.state.to_cartesian_state(mu); + let position = cartesian.position(); + Python::with_gil(|py| pyarray![py, position.x, position.y, position.z].into_py(py)) + } + + fn velocity(&self) -> Py> { + let mu = self.origin.gravitational_parameter(); + let cartesian = self.state.to_cartesian_state(mu); + let velocity = cartesian.velocity(); + Python::with_gil(|py| pyarray![py, velocity.x, velocity.y, velocity.z].into_py(py)) + } + + fn cartesian(&self) -> (Py>, Py>) { + (self.position(), self.velocity()) + } + + fn keplerian(&self) -> Py> { + Python::with_gil(|py| { + pyarray![ + py, + self.state.true_anomaly(), + self.state.eccentricity(), + self.state.inclination(), + self.state.ascending_node(), + self.state.periapsis_arg(), + self.state.true_anomaly() + ] + .into_py(py) + }) + } + + fn semi_major(&self) -> f64 { + self.state.semi_major() + } + + fn eccentricity(&self) -> f64 { + self.state.eccentricity() + } + + fn inclination(&self) -> f64 { + self.state.inclination() + } + + fn ascending_node(&self) -> f64 { + self.state.ascending_node() + } + + fn periapsis_arg(&self) -> f64 { + self.state.periapsis_arg() + } + + fn true_anomaly(&self) -> f64 { + self.state.true_anomaly() + } +} diff --git a/crates/lox_py/src/twobody.rs b/crates/lox_py/src/twobody.rs deleted file mode 100644 index 752cc8e2..00000000 --- a/crates/lox_py/src/twobody.rs +++ /dev/null @@ -1,177 +0,0 @@ -/* - * Copyright (c) 2024. Helge Eichhorn and the LOX contributors - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, you can obtain one at https://mozilla.org/MPL/2.0/. - */ - -use numpy::{pyarray, PyArray1}; -use pyo3::exceptions::PyValueError; -use pyo3::prelude::*; - -use lox_core::time::epochs::Epoch; -use lox_core::two_body::elements::cartesian_to_keplerian; -use lox_core::two_body::{DVec3, Elements, TwoBody}; - -use crate::bodies::{PyBarycenter, PyBody, PyMinorBody, PyPlanet, PySatellite, PySun}; -use crate::time::PyEpoch; - -#[pyclass(name = "Cartesian")] -pub struct PyCartesian { - t: Epoch, - x: f64, - y: f64, - z: f64, - vx: f64, - vy: f64, - vz: f64, - center: PyBody, -} - -impl TwoBody for PyCartesian { - fn epoch(&self) -> Epoch { - self.t - } - - fn position(&self) -> DVec3 { - DVec3::new(self.x, self.y, self.z) - } - - fn velocity(&self) -> DVec3 { - DVec3::new(self.vx, self.vy, self.vz) - } - - fn cartesian(&self) -> (DVec3, DVec3) { - (TwoBody::position(self), TwoBody::velocity(self)) - } - - fn keplerian(&self) -> Elements { - let mu = match &self.center { - PyBody::Barycenter(barycenter) => barycenter.gravitational_parameter(), - PyBody::Sun(sun) => sun.gravitational_parameter(), - PyBody::Planet(planet) => planet.gravitational_parameter(), - PyBody::Satellite(satellite) => satellite.gravitational_parameter(), - PyBody::MinorBody(minor_body) => minor_body.gravitational_parameter(), - }; - cartesian_to_keplerian(mu, TwoBody::position(self), TwoBody::velocity(self)) - } - - fn semi_major(&self) -> f64 { - self.keplerian().0 - } - - fn eccentricity(&self) -> f64 { - self.keplerian().1 - } - - fn inclination(&self) -> f64 { - self.keplerian().2 - } - - fn ascending_node(&self) -> f64 { - self.keplerian().3 - } - - fn periapsis_arg(&self) -> f64 { - self.keplerian().4 - } - - fn true_anomaly(&self) -> f64 { - self.keplerian().5 - } -} - -#[pymethods] -impl PyCartesian { - #[new] - #[allow(clippy::too_many_arguments)] - fn new( - t: &PyEpoch, - x: f64, - y: f64, - z: f64, - vx: f64, - vy: f64, - vz: f64, - body: PyObject, - ) -> PyResult { - let body = Python::with_gil(|py| { - if let Ok(body) = body.extract::(py) { - Some(PyBody::Barycenter(body)) - } else if let Ok(body) = body.extract::(py) { - Some(PyBody::Sun(body)) - } else if let Ok(body) = body.extract::(py) { - Some(PyBody::Planet(body)) - } else if let Ok(body) = body.extract::(py) { - Some(PyBody::Satellite(body)) - } else if let Ok(body) = body.extract::(py) { - Some(PyBody::MinorBody(body)) - } else { - None - } - }); - match body { - None => Err(PyValueError::new_err("foo")), - Some(body) => Ok(Self { - t: t.0, - x, - y, - z, - vx, - vy, - vz, - center: body, - }), - } - } - fn center(&self) -> PyObject { - Python::with_gil(|py| match &self.center { - PyBody::Barycenter(barycenter) => barycenter.clone().into_py(py), - PyBody::Sun(sun) => sun.clone().into_py(py), - PyBody::Planet(planet) => planet.clone().into_py(py), - PyBody::Satellite(satellite) => satellite.clone().into_py(py), - PyBody::MinorBody(minor_body) => minor_body.clone().into_py(py), - }) - } - - fn position(&self) -> Py> { - Python::with_gil(|py| pyarray![py, self.x, self.y, self.z].into_py(py)) - } - - fn velocity(&self) -> Py> { - Python::with_gil(|py| pyarray![py, self.vx, self.vy, self.vz].into_py(py)) - } - - fn cartesian(&self) -> (Py>, Py>) { - (self.position(), self.velocity()) - } - - fn keplerian(&self) -> Elements { - TwoBody::keplerian(self) - } - - fn semi_major(&self) -> f64 { - self.keplerian().0 - } - - fn eccentricity(&self) -> f64 { - self.keplerian().1 - } - - fn inclination(&self) -> f64 { - self.keplerian().2 - } - - fn ascending_node(&self) -> f64 { - self.keplerian().3 - } - - fn periapsis_arg(&self) -> f64 { - self.keplerian().4 - } - - fn true_anomaly(&self) -> f64 { - self.keplerian().5 - } -}