diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 467b8bf9..bd564d79 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -16,6 +16,29 @@ jobs: - name: Run tests run: cargo test + python: + name: Python smoke tests + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: dtolnay/rust-toolchain@stable + - uses: Swatinem/rust-cache@v2 + - uses: actions/setup-python@v5 + with: + python-version: '3.11' + - name: Generate virtual environment + run: python -m venv .venv + - name: Build Python wrapper + uses: PyO3/maturin-action@v1 + with: + command: develop + working-directory: ./crates/lox_py + args: --extras dev + - name: Run test + run: | + source .venv/bin/activate + pytest crates/lox_py/tests + fmt: name: Rustfmt runs-on: ubuntu-latest diff --git a/.idea/lox-space.iml b/.idea/lox-space.iml index 79dd551d..b2060bea 100644 --- a/.idea/lox-space.iml +++ b/.idea/lox-space.iml @@ -2,7 +2,7 @@ - + @@ -18,11 +18,12 @@ + - + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 00000000..6e1ff8df --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,6 @@ + + + + + \ No newline at end of file diff --git a/Cargo.lock b/Cargo.lock index 20b56ee5..0b872a6c 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -90,28 +90,35 @@ checksum = "baf0a07a401f374238ab8e2f11a104d2851bf9ce711ec69804834de8af45c7af" [[package]] name = "divan" -version = "0.1.2" +version = "0.1.11" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fab20f5802e0b897093184f5dc5d23447fa715604f238dc798f6da188b230019" +checksum = "5398159ee27f2b123d89b856bad61725442f37df5fb98c30cd570c318d594aee" dependencies = [ + "cfg-if", "clap", "condtype", "divan-macros", - "linkme", + "libc", "regex-lite", ] [[package]] name = "divan-macros" -version = "0.1.2" +version = "0.1.11" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1c5d6354551e0b5c451a948814fc47fe745a14eac7835c087d60162661019db4" +checksum = "5092f66eb3563a01e85552731ae82c04c934ff4efd7ad1a0deae7b948f4b3ec4" dependencies = [ "proc-macro2", "quote", "syn", ] +[[package]] +name = "dyn-clone" +version = "1.0.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "545b22097d44f8a9581187cdf93de7a71e4722bf51200cfaba810865b49a495d" + [[package]] name = "errno" version = "0.3.5" @@ -257,9 +264,9 @@ dependencies = [ [[package]] name = "glam" -version = "0.24.2" +version = "0.25.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b5418c17512bdf42730f9032c74e1ae39afc408745ebb2acf72fbc4691c17945" +checksum = "151665d9be52f9bb40fc7966565d39666f2d1e69233571b71b87791c7e0528b3" [[package]] name = "glob" @@ -297,26 +304,6 @@ version = "0.2.8" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058" -[[package]] -name = "linkme" -version = "0.3.17" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "91ed2ee9464ff9707af8e9ad834cffa4802f072caad90639c583dd3c62e6e608" -dependencies = [ - "linkme-impl", -] - -[[package]] -name = "linkme-impl" -version = "0.3.17" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ba125974b109d512fccbc6c0244e7580143e460895dfd6ea7f8bbb692fd94396" -dependencies = [ - "proc-macro2", - "quote", - "syn", -] - [[package]] name = "linux-raw-sys" version = "0.4.10" @@ -345,6 +332,7 @@ name = "lox_core" version = "0.1.0" dependencies = [ "divan", + "dyn-clone", "fast_polynomial", "float_eq", "glam", @@ -366,8 +354,10 @@ dependencies = [ name = "lox_py" version = "0.1.0" dependencies = [ + "float_eq", "lox_core", "pyo3", + "rstest", "thiserror", ] @@ -537,9 +527,9 @@ dependencies = [ [[package]] name = "proptest" -version = "1.3.1" +version = "1.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7c003ac8c77cb07bb74f5f198bce836a689bcd5a42574612bf14d17bfd08c20e" +checksum = "31b476131c3c86cb68032fdc5cb6d5a1045e3e42d96b69fa599fd77701e1f5bf" dependencies = [ "bit-set", "bit-vec", @@ -549,7 +539,7 @@ dependencies = [ "rand", "rand_chacha", "rand_xorshift", - "regex-syntax 0.7.5", + "regex-syntax", "rusty-fork", "tempfile", "unarray", @@ -557,9 +547,9 @@ dependencies = [ [[package]] name = "pyo3" -version = "0.20.0" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "04e8453b658fe480c3e70c8ed4e3d3ec33eb74988bd186561b0cc66b85c3bc4b" +checksum = "9a89dc7a5850d0e983be1ec2a463a171d20990487c3cfcd68b5363f1ee3d6fe0" dependencies = [ "cfg-if", "indoc", @@ -574,9 +564,9 @@ dependencies = [ [[package]] name = "pyo3-build-config" -version = "0.20.0" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a96fe70b176a89cff78f2fa7b3c930081e163d5379b4dcdf993e3ae29ca662e5" +checksum = "07426f0d8fe5a601f26293f300afd1a7b1ed5e78b2a705870c5f30893c5163be" dependencies = [ "once_cell", "target-lexicon", @@ -584,9 +574,9 @@ dependencies = [ [[package]] name = "pyo3-ffi" -version = "0.20.0" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "214929900fd25e6604661ed9cf349727c8920d47deff196c4e28165a6ef2a96b" +checksum = "dbb7dec17e17766b46bca4f1a4215a85006b4c2ecde122076c562dd058da6cf1" dependencies = [ "libc", "pyo3-build-config", @@ -594,9 +584,9 @@ dependencies = [ [[package]] name = "pyo3-macros" -version = "0.20.0" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dac53072f717aa1bfa4db832b39de8c875b7c7af4f4a6fe93cdbf9264cf8383b" +checksum = "05f738b4e40d50b5711957f142878cfa0f28e054aa0ebdfc3fd137a843f74ed3" dependencies = [ "proc-macro2", "pyo3-macros-backend", @@ -606,9 +596,9 @@ dependencies = [ [[package]] name = "pyo3-macros-backend" -version = "0.20.0" +version = "0.20.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7774b5a8282bd4f25f803b1f0d945120be959a36c72e08e7cd031c792fdfd424" +checksum = "0fc910d4851847827daf9d6cdd4a823fbdaab5b8818325c5e97a86da79e8881f" dependencies = [ "heck", "proc-macro2", @@ -688,7 +678,7 @@ dependencies = [ "aho-corasick", "memchr", "regex-automata", - "regex-syntax 0.8.2", + "regex-syntax", ] [[package]] @@ -699,7 +689,7 @@ checksum = "5011c7e263a695dc8ca064cddb722af1be54e517a280b12a5356f98366899e5d" dependencies = [ "aho-corasick", "memchr", - "regex-syntax 0.8.2", + "regex-syntax", ] [[package]] @@ -708,12 +698,6 @@ version = "0.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "30b661b2f27137bdbc16f00eda72866a92bb28af1753ffbd56744fb6e2e9cd8e" -[[package]] -name = "regex-syntax" -version = "0.7.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dbb5fb1acd8a1a18b3dd5be62d25485eb770e05afb408a9627d14d451bae12da" - [[package]] name = "regex-syntax" version = "0.8.2" diff --git a/Cargo.toml b/Cargo.toml index 9d6637f8..333d4c50 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -15,5 +15,8 @@ authors = ["Helge Eichhorn and the lox-space contributors"] [workspace.dependencies] lox_core = { path = "./crates/lox_core" } +float_eq = "1.0.1" +proptest = "1.4.0" +rstest = "0.18.2" thiserror = "1.0" diff --git a/crates/lox-space/examples/iss.rs b/crates/lox-space/examples/iss.rs index 1e137215..66767908 100644 --- a/crates/lox-space/examples/iss.rs +++ b/crates/lox-space/examples/iss.rs @@ -14,20 +14,21 @@ 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 = Cartesian::new(epoch, Earth, Icrf, position, velocity); + let iss = Keplerian::from(iss_cartesian); - 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!("Semi-major axis: {:.3} km", iss.semi_major_axis()); println!("Eccentricity: {:.6}", iss.eccentricity()); println!("Inclination: {:.3}°", iss.inclination().to_degrees()); println!( "Longitude of ascending node: {:.3}°", iss.ascending_node().to_degrees() ); - println!("Argument of perigee: {}°", iss.periapsis_arg().to_degrees()); + println!( + "Argument of perigee: {}°", + iss.periapsis_argument().to_degrees() + ); println!("True anomaly: {:.3}°", iss.true_anomaly().to_degrees()); } diff --git a/crates/lox-space/src/prelude.rs b/crates/lox-space/src/prelude.rs index 20f7e2b7..b9f9c620 100644 --- a/crates/lox-space/src/prelude.rs +++ b/crates/lox-space/src/prelude.rs @@ -7,8 +7,8 @@ */ pub use lox_core::bodies::*; - +pub use lox_core::coords::two_body::{Cartesian, Keplerian}; +pub use lox_core::coords::DVec3; +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/Cargo.toml b/crates/lox_core/Cargo.toml index c3bdf10a..f66026a8 100644 --- a/crates/lox_core/Cargo.toml +++ b/crates/lox_core/Cargo.toml @@ -4,17 +4,18 @@ version = "0.1.0" edition = "2021" [dependencies] -float_eq = "1.0.1" -glam = "0.24.2" -num = "0.4.0" +float_eq.workspace = true thiserror.workspace = true +glam = "0.25.0" +num = "0.4.1" nom = "7.1.3" fast_polynomial = "0.1.0" +dyn-clone = "1.0.16" [dev-dependencies] -proptest = "1.1.0" -rstest = "0.18.2" -divan = "0.1.2" +proptest.workspace = true +rstest.workspace = true +divan = "0.1.11" [[bench]] name = "iau_frames" diff --git a/crates/lox_core/src/bodies.rs b/crates/lox_core/src/bodies.rs index d8420438..58d0c0dd 100644 --- a/crates/lox_core/src/bodies.rs +++ b/crates/lox_core/src/bodies.rs @@ -6,6 +6,7 @@ * file, you can obtain one at https://mozilla.org/MPL/2.0/. */ +use dyn_clone::{clone_trait_object, DynClone}; use std::f64::consts::PI; use std::fmt::{Display, Formatter}; @@ -91,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()) @@ -103,7 +126,8 @@ macro_rules! body { body! { Sun, 10 } // Planets. -pub trait Planet: PointMass + Spheroid {} +pub trait Planet: PointMass + Spheroid + DynClone {} +clone_trait_object!(Planet); body! { Mercury, Planet, 199 } body! { Venus, Planet, 299 } @@ -116,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 { @@ -134,7 +161,8 @@ impl PointMass for SolarSystemBarycenter { } // Satellites. -pub trait Satellite: PointMass + TriAxial {} +pub trait Satellite: PointMass + TriAxial + DynClone {} +clone_trait_object!(Satellite); body! { Moon, Satellite, 301 } body! { Phobos, Satellite, 401 } @@ -290,7 +318,8 @@ body! { Kerberos, 904 } body! { Styx, 905 } // Minor bodies. -pub trait MinorBody: PointMass + TriAxial {} +pub trait MinorBody: PointMass + TriAxial + DynClone {} +clone_trait_object!(MinorBody); body! {Gaspra, 9511010 } body! {Ida, 2431010 } diff --git a/crates/lox_core/src/coords.rs b/crates/lox_core/src/coords.rs new file mode 100644 index 00000000..a59cce1b --- /dev/null +++ b/crates/lox_core/src/coords.rs @@ -0,0 +1,24 @@ +/* + * Copyright (c) 2023. 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/. + */ + +pub use glam::DVec3; + +use crate::bodies::PointMass; +use crate::frames::ReferenceFrame; + +pub mod anomalies; +pub mod states; +pub mod two_body; + +pub trait CoordinateSystem { + type Origin: PointMass; + type Frame: ReferenceFrame; + + fn origin(&self) -> Self::Origin; + fn reference_frame(&self) -> Self::Frame; +} diff --git a/crates/lox_core/src/coords/anomalies.rs b/crates/lox_core/src/coords/anomalies.rs new file mode 100644 index 00000000..661f5ed6 --- /dev/null +++ b/crates/lox_core/src/coords/anomalies.rs @@ -0,0 +1,44 @@ +/* + * Copyright (c) 2023. 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/. + */ + +pub fn hyperbolic_to_true(hyperbolic_anomaly: f64, eccentricity: f64) -> f64 { + 2.0 * (((1.0 + eccentricity) / (eccentricity - 1.0)).sqrt() * (hyperbolic_anomaly / 2.0).tanh()) + .atan() +} + +pub fn eccentric_to_true(eccentric_anomaly: f64, eccentricity: f64) -> f64 { + 2.0 * (((1.0 + eccentricity) / (1.0 - eccentricity)).sqrt() * (eccentric_anomaly / 2.0).tan()) + .atan() +} + +#[cfg(test)] +mod tests { + use std::f64::consts::PI; + + use float_eq::assert_float_eq; + + use super::*; + + #[test] + fn test_hyperbolic() { + assert_float_eq!( + hyperbolic_to_true(PI / 2.0, 1.2), + 2.2797028138935547, + rel <= 1e-8 + ); + } + + #[test] + fn test_eccentric() { + assert_float_eq!( + eccentric_to_true(PI / 2.0, 0.2), + 1.7721542475852272, + rel <= 1e-8 + ); + } +} diff --git a/crates/lox_core/src/coords/states.rs b/crates/lox_core/src/coords/states.rs new file mode 100644 index 00000000..94900e47 --- /dev/null +++ b/crates/lox_core/src/coords/states.rs @@ -0,0 +1,464 @@ +/* + * 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 float_eq::float_eq; +use glam::{DMat3, DVec3}; + +use crate::math::{mod_two_pi, normalize_two_pi}; +use crate::time::epochs::Epoch; + +pub trait TwoBodyState { + fn time(&self) -> Epoch; + fn to_cartesian_state(&self, grav_param: f64) -> CartesianState; + fn to_keplerian_state(&self, grav_param: f64) -> KeplerianState; +} + +#[derive(Debug, Copy, Clone, PartialEq)] +pub struct CartesianState { + time: Epoch, + position: DVec3, + velocity: DVec3, +} + +impl CartesianState { + pub fn new(time: Epoch, position: DVec3, velocity: DVec3) -> Self { + Self { + time, + position, + velocity, + } + } + + pub fn position(&self) -> DVec3 { + self.position + } + + pub fn velocity(&self) -> DVec3 { + self.velocity + } +} + +impl TwoBodyState for CartesianState { + fn time(&self) -> Epoch { + self.time + } + + fn to_cartesian_state(&self, _grav_param: f64) -> CartesianState { + *self + } + + fn to_keplerian_state(&self, grav_param: f64) -> KeplerianState { + let r = self.position.length(); + let v = self.velocity.length(); + let h = self.position.cross(self.velocity); + let hm = h.length(); + let node = DVec3::Z.cross(h); + let e = eccentricity_vector(grav_param, self.position, self.velocity); + let eccentricity = e.length(); + let inclination = h.angle_between(DVec3::Z); + + let equatorial = is_equatorial(inclination); + let circular = is_circular(eccentricity); + + let semi_major = if circular { + hm.powi(2) / grav_param + } else { + -grav_param / (2.0 * (v.powi(2) / 2.0 - grav_param / r)) + }; + + let ascending_node; + let periapsis_arg; + let true_anomaly; + if equatorial && !circular { + ascending_node = 0.0; + periapsis_arg = azimuth(e); + true_anomaly = (h.dot(e.cross(self.position)) / hm).atan2(self.position.dot(e)); + } else if !equatorial && circular { + ascending_node = azimuth(node); + periapsis_arg = 0.0; + true_anomaly = (self.position.dot(h.cross(node)) / hm).atan2(self.position.dot(node)); + } else if equatorial && circular { + ascending_node = 0.0; + periapsis_arg = 0.0; + true_anomaly = azimuth(self.position); + } else { + if semi_major > 0.0 { + let e_se = self.position.dot(self.velocity) / (grav_param * semi_major).sqrt(); + let e_ce = (r * v.powi(2)) / grav_param - 1.0; + true_anomaly = + crate::coords::anomalies::eccentric_to_true(e_se.atan2(e_ce), eccentricity); + } else { + let e_sh = self.position.dot(self.velocity) / (-grav_param * semi_major).sqrt(); + let e_ch = (r * v.powi(2)) / grav_param - 1.0; + true_anomaly = crate::coords::anomalies::hyperbolic_to_true( + ((e_ch + e_sh) / (e_ch - e_sh)).ln() / 2.0, + eccentricity, + ); + } + ascending_node = azimuth(node); + let px = self.position.dot(node); + let py = self.position.dot(h.cross(node)) / hm; + periapsis_arg = py.atan2(px) - true_anomaly; + } + + KeplerianState::new( + self.time, + semi_major, + eccentricity, + inclination, + mod_two_pi(ascending_node), + mod_two_pi(periapsis_arg), + normalize_two_pi(true_anomaly, 0.0), + ) + } +} + +#[derive(Debug, Copy, Clone, PartialEq)] +pub struct KeplerianState { + time: Epoch, + semi_major: f64, + eccentricity: f64, + inclination: f64, + ascending_node: f64, + periapsis_argument: f64, + true_anomaly: f64, +} + +impl KeplerianState { + pub fn new( + time: Epoch, + semi_major: f64, + eccentricity: f64, + inclination: f64, + ascending_node: f64, + periapsis_argument: f64, + true_anomaly: f64, + ) -> Self { + Self { + time, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_argument, + true_anomaly, + } + } + + pub fn semi_major_axis(&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_argument(&self) -> f64 { + self.periapsis_argument + } + + pub fn true_anomaly(&self) -> f64 { + self.true_anomaly + } + + pub fn semi_latus_rectum(&self) -> f64 { + if float_eq!(self.eccentricity, 0.0, abs <= 1e-3) { + self.semi_major + } else { + self.semi_major * (1.0 - self.eccentricity.powi(2)) + } + } + + pub fn to_perifocal(&self, grav_param: f64) -> (DVec3, DVec3) { + let semi_latus = self.semi_latus_rectum(); + let (sin_nu, cos_nu) = self.true_anomaly.sin_cos(); + let sqrt_mu_p = (grav_param / semi_latus).sqrt(); + + let pos = + DVec3::new(cos_nu, sin_nu, 0.0) * (semi_latus / (1.0 + self.eccentricity * cos_nu)); + let vel = DVec3::new(-sin_nu, self.eccentricity + cos_nu, 0.0) * sqrt_mu_p; + + (pos, vel) + } +} + +impl TwoBodyState for KeplerianState { + fn time(&self) -> Epoch { + self.time + } + + fn to_cartesian_state(&self, grav_param: f64) -> CartesianState { + let (pos, vel) = self.to_perifocal(grav_param); + let rot = DMat3::from_rotation_z(self.ascending_node) + * DMat3::from_rotation_x(self.inclination) + * DMat3::from_rotation_z(self.periapsis_argument); + CartesianState::new(self.time, rot * pos, rot * vel) + } + + fn to_keplerian_state(&self, _grav_param: f64) -> KeplerianState { + *self + } +} + +fn azimuth(v: DVec3) -> f64 { + v.y.atan2(v.x) +} + +fn eccentricity_vector(grav_param: f64, pos: DVec3, vel: DVec3) -> DVec3 { + (pos * (vel.dot(vel) - grav_param / pos.length()) - vel * pos.dot(vel)) / grav_param +} + +fn is_equatorial(inclination: f64) -> bool { + float_eq!(inclination.abs(), 0.0, abs <= 1e-3) +} + +fn is_circular(eccentricity: f64) -> bool { + float_eq!(eccentricity, 0.0, abs <= 1e-3) +} + +#[cfg(test)] +mod tests { + use crate::time::epochs::TimeScale; + use float_eq::assert_float_eq; + use glam::DVec3; + + use super::*; + + #[test] + fn test_elliptic() { + let time = Epoch::j2000(TimeScale::TDB); + let grav_param = 3.9860047e14; + let semi_major = 24464560.0; + let eccentricity = 0.7311; + let inclination = 0.122138; + let ascending_node = 1.00681; + let periapsis_arg = 3.10686; + let true_anomaly = 0.44369564302687126; + let pos = DVec3::new( + -0.107622532467967e7, + -0.676589636432773e7, + -0.332308783350379e6, + ); + let vel = DVec3::new( + 0.935685775154103e4, + -0.331234775037644e4, + -0.118801577532701e4, + ); + + let cartesian = CartesianState::new(time, pos, vel); + assert_eq!(cartesian.to_cartesian_state(grav_param), cartesian); + + let keplerian = KeplerianState::new( + time, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ); + assert_eq!(keplerian.to_keplerian_state(grav_param), keplerian); + + let cartesian1 = keplerian.to_cartesian_state(grav_param); + let keplerian1 = cartesian.to_keplerian_state(grav_param); + + assert_eq!(cartesian1.time(), time); + assert_eq!(keplerian1.time(), time); + + assert_float_eq!(pos.x, cartesian1.position.x, rel <= 1e-8); + assert_float_eq!(pos.y, cartesian1.position.y, rel <= 1e-8); + assert_float_eq!(pos.z, cartesian1.position.z, rel <= 1e-8); + assert_float_eq!(vel.x, cartesian1.velocity.x, rel <= 1e-8); + assert_float_eq!(vel.y, cartesian1.velocity.y, rel <= 1e-8); + assert_float_eq!(vel.z, cartesian1.velocity.z, rel <= 1e-8); + + assert_float_eq!(semi_major, keplerian1.semi_major, rel <= 1e-8); + assert_float_eq!(eccentricity, keplerian1.eccentricity, rel <= 1e-8); + assert_float_eq!(inclination, keplerian1.inclination, rel <= 1e-8); + assert_float_eq!(ascending_node, keplerian1.ascending_node, rel <= 1e-8); + assert_float_eq!(periapsis_arg, keplerian1.periapsis_argument, rel <= 1e-8); + assert_float_eq!(true_anomaly, keplerian1.true_anomaly, rel <= 1e-8); + } + + #[test] + fn test_circular() { + let time = Epoch::j2000(TimeScale::TDB); + let grav_param = 3.986004418e14; + let semi_major = 6778136.6; + let eccentricity = 0.0; + let inclination = 15f64.to_radians(); + let ascending_node = 20f64.to_radians(); + let periapsis_arg = 0.0; + let true_anomaly = 30f64.to_radians(); + let pos = DVec3::new(4396398.60746266, 5083838.45333733, 877155.42119322); + let vel = DVec3::new(-5797.06004014, 4716.60916063, 1718.86034246); + let cartesian = CartesianState::new(time, pos, vel); + let keplerian = KeplerianState::new( + time, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ); + + let cartesian1 = keplerian.to_cartesian_state(grav_param); + let keplerian1 = cartesian.to_keplerian_state(grav_param); + + assert_float_eq!(pos.x, cartesian1.position.x, rel <= 1e-8); + assert_float_eq!(pos.y, cartesian1.position.y, rel <= 1e-8); + assert_float_eq!(pos.z, cartesian1.position.z, rel <= 1e-8); + assert_float_eq!(vel.x, cartesian1.velocity.x, rel <= 1e-8); + assert_float_eq!(vel.y, cartesian1.velocity.y, rel <= 1e-8); + assert_float_eq!(vel.z, cartesian1.velocity.z, rel <= 1e-8); + + assert_float_eq!(semi_major, keplerian1.semi_major, rel <= 1e-8); + assert_float_eq!(eccentricity, keplerian1.eccentricity, abs <= 1e-8); + assert_float_eq!(inclination, keplerian1.inclination, rel <= 1e-8); + assert_float_eq!(ascending_node, keplerian1.ascending_node, rel <= 1e-8); + assert_float_eq!(periapsis_arg, keplerian1.periapsis_argument, rel <= 1e-8); + assert_float_eq!(true_anomaly, keplerian1.true_anomaly, rel <= 1e-8); + } + + #[test] + fn test_circular_orekit() { + let time = Epoch::j2000(TimeScale::TDB); + let grav_param = 3.9860047e14; + let semi_major = 24464560.0; + let eccentricity = 0.0; + let inclination = 0.122138; + let ascending_node = 1.00681; + let periapsis_arg = 0.0; + let true_anomaly = 0.048363; + let keplerian = KeplerianState::new( + time, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ); + + let keplerian1 = keplerian + .to_cartesian_state(grav_param) + .to_keplerian_state(grav_param); + + assert_float_eq!(semi_major, keplerian1.semi_major, rel <= 1e-8); + assert_float_eq!(eccentricity, keplerian1.eccentricity, abs <= 1e-8); + assert_float_eq!(inclination, keplerian1.inclination, rel <= 1e-8); + assert_float_eq!(ascending_node, keplerian1.ascending_node, rel <= 1e-8); + assert_float_eq!(periapsis_arg, keplerian1.periapsis_argument, rel <= 1e-8); + assert_float_eq!(true_anomaly, keplerian1.true_anomaly, rel <= 1e-8); + } + + #[test] + fn test_hyperbolic_orekit() { + let time = Epoch::j2000(TimeScale::TDB); + let grav_param = 3.9860047e14; + let semi_major = -24464560.0; + let eccentricity = 1.7311; + let inclination = 0.122138; + let ascending_node = 1.00681; + let periapsis_arg = 3.10686; + let true_anomaly = 0.12741601769795755; + let keplerian = KeplerianState::new( + time, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ); + + let keplerian1 = keplerian + .to_cartesian_state(grav_param) + .to_keplerian_state(grav_param); + + assert_float_eq!(semi_major, keplerian1.semi_major, rel <= 1e-8); + assert_float_eq!(eccentricity, keplerian1.eccentricity, rel <= 1e-8); + assert_float_eq!(inclination, keplerian1.inclination, rel <= 1e-8); + assert_float_eq!(ascending_node, keplerian1.ascending_node, rel <= 1e-8); + assert_float_eq!(periapsis_arg, keplerian1.periapsis_argument, rel <= 1e-8); + assert_float_eq!(true_anomaly, keplerian1.true_anomaly, rel <= 1e-8); + } + + #[test] + fn test_equatorial() { + let time = Epoch::j2000(TimeScale::TDB); + let grav_param = 3.9860047e14; + let semi_major = 24464560.0; + let eccentricity = 0.7311; + let inclination = 0.0; + let ascending_node = 0.0; + let periapsis_arg = 3.10686; + let true_anomaly = 0.44369564302687126; + let keplerian = KeplerianState::new( + time, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ); + + let keplerian1 = keplerian + .to_cartesian_state(grav_param) + .to_keplerian_state(grav_param); + + assert_float_eq!(semi_major, keplerian1.semi_major, rel <= 1e-8); + assert_float_eq!(eccentricity, keplerian1.eccentricity, rel <= 1e-8); + assert_float_eq!(inclination, keplerian1.inclination, rel <= 1e-8); + assert_float_eq!(ascending_node, keplerian1.ascending_node, rel <= 1e-8); + assert_float_eq!(periapsis_arg, keplerian1.periapsis_argument, rel <= 1e-8); + assert_float_eq!(true_anomaly, keplerian1.true_anomaly, rel <= 1e-8); + } + + #[test] + fn test_circular_equatorial() { + let time = Epoch::j2000(TimeScale::TDB); + let grav_param = 3.9860047e14; + let semi_major = 24464560.0; + let eccentricity = 0.0; + let inclination = 0.0; + let ascending_node = 0.0; + let periapsis_arg = 0.0; + let true_anomaly = 0.44369564302687126; + let keplerian = KeplerianState::new( + time, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ); + + let keplerian1 = keplerian + .to_cartesian_state(grav_param) + .to_keplerian_state(grav_param); + + assert_float_eq!(semi_major, keplerian1.semi_major, rel <= 1e-8); + assert_float_eq!(eccentricity, keplerian1.eccentricity, abs <= 1e-8); + assert_float_eq!(inclination, keplerian1.inclination, rel <= 1e-8); + assert_float_eq!(ascending_node, keplerian1.ascending_node, rel <= 1e-8); + assert_float_eq!(periapsis_arg, keplerian1.periapsis_argument, rel <= 1e-8); + assert_float_eq!(true_anomaly, keplerian1.true_anomaly, rel <= 1e-8); + } +} diff --git a/crates/lox_core/src/coords/two_body.rs b/crates/lox_core/src/coords/two_body.rs new file mode 100644 index 00000000..8a15cf5a --- /dev/null +++ b/crates/lox_core/src/coords/two_body.rs @@ -0,0 +1,338 @@ +/* + * 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 glam::DVec3; + +use crate::bodies::PointMass; +use crate::coords::states::{CartesianState, KeplerianState, TwoBodyState}; +use crate::coords::CoordinateSystem; +use crate::frames::{InertialFrame, ReferenceFrame}; +use crate::time::epochs::Epoch; + +pub trait TwoBody +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + fn to_cartesian(&self) -> Cartesian; + + fn to_keplerian(&self) -> Keplerian; +} + +#[derive(Debug, Copy, Clone, PartialEq)] +pub struct Cartesian +where + T: PointMass + Copy, + S: ReferenceFrame + Copy, +{ + state: CartesianState, + origin: T, + frame: S, +} + +impl Cartesian +where + T: PointMass + Copy, + S: ReferenceFrame + Copy, +{ + pub fn new(time: Epoch, origin: T, frame: S, position: DVec3, velocity: DVec3) -> Self { + let state = CartesianState::new(time, position, velocity); + Self { + state, + origin, + frame, + } + } + + pub fn time(&self) -> Epoch { + self.state.time() + } + + pub fn position(&self) -> DVec3 { + self.state.position() + } + + pub fn velocity(&self) -> DVec3 { + self.state.position() + } +} + +impl TwoBody for Cartesian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + fn to_cartesian(&self) -> Cartesian { + *self + } + + fn to_keplerian(&self) -> Keplerian { + Keplerian::from(*self) + } +} + +impl CoordinateSystem for Cartesian +where + T: PointMass + Copy, + S: ReferenceFrame + Copy, +{ + type Origin = T; + type Frame = S; + + fn origin(&self) -> Self::Origin { + self.origin + } + + fn reference_frame(&self) -> Self::Frame { + self.frame + } +} + +impl From> for Cartesian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + fn from(keplerian: Keplerian) -> Self { + let grav_param = keplerian.origin.gravitational_parameter(); + let state = keplerian.state.to_cartesian_state(grav_param); + Cartesian { + state, + origin: keplerian.origin, + frame: keplerian.frame, + } + } +} + +#[derive(Debug, Copy, Clone, PartialEq)] +pub struct Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + state: KeplerianState, + origin: T, + frame: S, +} + +impl Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + #[allow(clippy::too_many_arguments)] + pub fn new( + time: Epoch, + origin: T, + frame: S, + semi_major: f64, + eccentricity: f64, + inclination: f64, + ascending_node: f64, + periapsis_arg: f64, + true_anomaly: f64, + ) -> Self { + let state = KeplerianState::new( + time, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ); + Self { + state, + origin, + frame, + } + } + + pub fn time(&self) -> Epoch { + self.state.time() + } + + pub fn semi_major_axis(&self) -> f64 { + self.state.semi_major_axis() + } + + pub fn eccentricity(&self) -> f64 { + self.state.eccentricity() + } + + pub fn inclination(&self) -> f64 { + self.state.inclination() + } + + pub fn ascending_node(&self) -> f64 { + self.state.ascending_node() + } + + pub fn periapsis_argument(&self) -> f64 { + self.state.periapsis_argument() + } + + pub fn true_anomaly(&self) -> f64 { + self.state.true_anomaly() + } +} + +impl TwoBody for Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + fn to_cartesian(&self) -> Cartesian { + Cartesian::from(*self) + } + + fn to_keplerian(&self) -> Keplerian { + *self + } +} + +impl CoordinateSystem for Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + type Origin = T; + type Frame = S; + + fn origin(&self) -> Self::Origin { + self.origin + } + + fn reference_frame(&self) -> Self::Frame { + self.frame + } +} + +impl From> for Keplerian +where + T: PointMass + Copy, + S: InertialFrame + Copy, +{ + fn from(cartesian: Cartesian) -> Self { + let grav_param = cartesian.origin.gravitational_parameter(); + let state = cartesian.state.to_keplerian_state(grav_param); + Self { + state, + origin: cartesian.origin, + frame: cartesian.frame, + } + } +} + +#[cfg(test)] +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; + + use super::*; + + #[test] + fn test_cartesian() { + let date = Date::new(2023, 3, 25).expect("Date should be valid"); + let time = Time::new(21, 8, 0).expect("Time should be valid"); + let epoch = Epoch::from_date_and_time(TimeScale::TDB, date, time); + let pos = DVec3::new( + -0.107622532467967e7, + -0.676589636432773e7, + -0.332308783350379e6, + ) * 1e-3; + let vel = DVec3::new( + 0.935685775154103e4, + -0.331234775037644e4, + -0.118801577532701e4, + ) * 1e-3; + + let cartesian = Cartesian::new(epoch, Earth, Icrf, pos, vel); + assert_eq!(cartesian.to_cartesian(), cartesian); + + let cartesian1 = cartesian.to_keplerian().to_cartesian(); + + assert_eq!(cartesian1.time(), epoch); + assert_eq!(cartesian1.origin(), Earth); + assert_eq!(cartesian1.reference_frame(), Icrf); + + 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); + } + + #[test] + fn test_keplerian() { + let date = Date::new(2023, 3, 25).expect("Date should be valid"); + let time = Time::new(21, 8, 0).expect("Time should be valid"); + let epoch = Epoch::from_date_and_time(TimeScale::TDB, date, time); + let semi_major = 24464560.0e-3; + let eccentricity = 0.7311; + let inclination = 0.122138; + let ascending_node = 1.00681; + let periapsis_arg = 3.10686; + let true_anomaly = 0.44369564302687126; + + let keplerian = Keplerian::new( + epoch, + Earth, + Icrf, + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ); + assert_eq!(keplerian.to_keplerian(), keplerian); + + let keplerian1 = keplerian.to_cartesian().to_keplerian(); + + assert_eq!(keplerian1.time(), epoch); + assert_eq!(keplerian1.origin(), Earth); + assert_eq!(keplerian1.reference_frame(), Icrf); + + assert_float_eq!( + keplerian.semi_major_axis(), + keplerian1.semi_major_axis(), + rel <= 1e-6 + ); + 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(), + rel <= 1e-6 + ); + assert_float_eq!( + keplerian.periapsis_argument(), + keplerian1.periapsis_argument(), + rel <= 1e-6 + ); + assert_float_eq!( + keplerian.true_anomaly(), + keplerian1.true_anomaly(), + rel <= 1e-6 + ); + } +} diff --git a/crates/lox_core/src/frames.rs b/crates/lox_core/src/frames.rs index 3dc8cd24..434d6b09 100644 --- a/crates/lox_core/src/frames.rs +++ b/crates/lox_core/src/frames.rs @@ -15,15 +15,29 @@ 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 } } -#[derive(Debug, Copy, Clone)] +#[derive(Debug, Copy, Clone, Eq, PartialEq)] pub struct Icrf; impl Display for Icrf { @@ -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); @@ -111,3 +122,15 @@ where self.rotation_from(frame, t).transpose() } } + +#[cfg(test)] +mod test { + use super::*; + + #[test] + fn test_icrf() { + assert!(Icrf.is_inertial()); + assert!(!Icrf.is_rotating()); + assert_eq!(format!("{}", Icrf), "ICRF"); + } +} diff --git a/crates/lox_core/src/frames/iau.rs b/crates/lox_core/src/frames/iau.rs index 91d700ae..2a788dd6 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)] +#[derive(Debug, Copy, Clone, Eq, PartialEq)] 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/lib.rs b/crates/lox_core/src/lib.rs index 0a10ad1e..4e170f47 100644 --- a/crates/lox_core/src/lib.rs +++ b/crates/lox_core/src/lib.rs @@ -8,11 +8,11 @@ pub mod bodies; pub mod communications; +pub mod coords; pub mod ephemeris; pub mod errors; pub mod frames; pub mod time; -pub mod two_body; pub mod types; mod earth; diff --git a/crates/lox_core/src/math.rs b/crates/lox_core/src/math.rs index 365146f5..13c9b5c3 100644 --- a/crates/lox_core/src/math.rs +++ b/crates/lox_core/src/math.rs @@ -24,6 +24,16 @@ pub(crate) fn arcsec_to_rad(arcsec: Arcsec) -> Radians { arcsec * RADIANS_IN_ARCSECOND } +/// Modulus after division by 2π, returning in the range [0,2π). +pub fn mod_two_pi(a: f64) -> f64 { + let w = a % (2.0 * PI); + if w < 0.0 { + w + 2.0 * PI + } else { + w + } +} + #[cfg(test)] mod tests { use float_eq::assert_float_eq; diff --git a/crates/lox_core/src/time/epochs.rs b/crates/lox_core/src/time/epochs.rs index 3e4b337d..c8de3c80 100644 --- a/crates/lox_core/src/time/epochs.rs +++ b/crates/lox_core/src/time/epochs.rs @@ -16,7 +16,7 @@ use crate::time::constants::i64::{SECONDS_PER_DAY, SECONDS_PER_HOUR, SECONDS_PER use crate::time::dates::Calendar::ProlepticJulian; use crate::time::dates::{Date, DateTime, Time}; -#[derive(Debug, Copy, Clone)] +#[derive(Debug, Copy, Clone, Eq, PartialEq)] pub enum TimeScale { TAI, TCB, diff --git a/crates/lox_core/src/two_body.rs b/crates/lox_core/src/two_body.rs deleted file mode 100644 index 9b6ed890..00000000 --- a/crates/lox_core/src/two_body.rs +++ /dev/null @@ -1,339 +0,0 @@ -/* - * Copyright (c) 2023. 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/. - */ - -pub use glam::DVec3; - -use crate::bodies::PointMass; -use crate::time::epochs::Epoch; -use crate::two_body::elements::{cartesian_to_keplerian, keplerian_to_cartesian}; - -pub mod elements; - -type Elements = (f64, f64, f64, f64, f64, f64); - -pub trait TwoBody { - type Center; - fn epoch(&self) -> Epoch; - fn center(&self) -> Self::Center; - 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; - fn ascending_node(&self) -> f64; - fn periapsis_arg(&self) -> f64; - fn true_anomaly(&self) -> f64; -} - -#[derive(Debug, Clone)] -pub struct Cartesian { - epoch: Epoch, - center: T, - position: DVec3, - velocity: DVec3, -} - -impl Cartesian { - pub fn new(epoch: Epoch, center: T, position: DVec3, velocity: DVec3) -> Self { - Self { - epoch, - center, - position, - velocity, - } - } -} - -impl TwoBody for Cartesian { - type Center = T; - - fn epoch(&self) -> Epoch { - self.epoch - } - - fn center(&self) -> Self::Center { - self.center - } - - fn position(&self) -> DVec3 { - self.position - } - - fn velocity(&self) -> DVec3 { - self.velocity - } - - fn cartesian(&self) -> (DVec3, DVec3) { - (self.position, self.velocity) - } - - fn keplerian(&self) -> Elements { - let mu = self.center.gravitational_parameter(); - cartesian_to_keplerian(mu, self.position, self.velocity) - } - - 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 - } -} - -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) - } -} - -#[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, -} - -impl Keplerian { - pub fn new(epoch: Epoch, center: T, elements: Elements) -> Self { - let (semi_major, eccentricity, inclination, ascending_node, periapsis_arg, true_anomaly) = - elements; - Self { - epoch, - center, - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, - } - } -} - -impl TwoBody for Keplerian { - type Center = T; - - fn epoch(&self) -> Epoch { - self.epoch - } - - fn center(&self) -> Self::Center { - self.center - } - - fn position(&self) -> DVec3 { - self.cartesian().0 - } - - fn velocity(&self) -> DVec3 { - self.cartesian().1 - } - - 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 keplerian(&self) -> Elements { - ( - self.semi_major, - self.eccentricity, - self.inclination, - self.ascending_node, - self.periapsis_arg, - self.true_anomaly, - ) - } - - fn semi_major(&self) -> f64 { - self.semi_major - } - - fn eccentricity(&self) -> f64 { - self.eccentricity - } - - fn inclination(&self) -> f64 { - self.inclination - } - - fn ascending_node(&self) -> f64 { - self.ascending_node - } - - fn periapsis_arg(&self) -> f64 { - self.periapsis_arg - } - - fn true_anomaly(&self) -> f64 { - self.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) - } -} - -#[cfg(test)] -mod tests { - use std::ops::Mul; - - use float_eq::assert_float_eq; - - use crate::bodies::Earth; - use crate::time::dates::{Date, Time}; - use crate::time::epochs::TimeScale; - - use super::*; - - #[test] - fn test_two_body() { - let date = Date::new(2023, 3, 25).expect("Date should be valid"); - let time = Time::new(21, 8, 0).expect("Time should be valid"); - let epoch = Epoch::from_date_and_time(TimeScale::TDB, date, time); - let semi_major = 24464560.0e-3; - let eccentricity = 0.7311; - let inclination = 0.122138; - let ascending_node = 1.00681; - let periapsis_arg = 3.10686; - let true_anomaly = 0.44369564302687126; - let pos = DVec3::new( - -0.107622532467967e7, - -0.676589636432773e7, - -0.332308783350379e6, - ) - .mul(1e-3); - let vel = DVec3::new( - 0.935685775154103e4, - -0.331234775037644e4, - -0.118801577532701e4, - ) - .mul(1e-3); - - let cartesian = Cartesian::new(epoch, Earth, pos, vel); - - assert_eq!( - cartesian.cartesian(), - (cartesian.position(), cartesian.velocity()) - ); - assert_eq!(cartesian.epoch(), epoch); - assert_eq!(cartesian.center(), Earth); - assert_eq!(cartesian.position(), pos); - assert_eq!(cartesian.velocity(), vel); - assert_float_eq!(cartesian.semi_major(), semi_major, rel <= 1e-6); - assert_float_eq!(cartesian.eccentricity(), eccentricity, rel <= 1e-6); - assert_float_eq!(cartesian.inclination(), inclination, rel <= 1e-6); - assert_float_eq!(cartesian.ascending_node(), ascending_node, rel <= 1e-6); - assert_float_eq!(cartesian.periapsis_arg(), periapsis_arg, rel <= 1e-6); - assert_float_eq!(cartesian.true_anomaly(), true_anomaly, rel <= 1e-6); - - let keplerian = Keplerian::new( - epoch, - Earth, - ( - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, - ), - ); - - assert_eq!( - keplerian.keplerian(), - ( - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly - ) - ); - assert_eq!(keplerian.epoch(), epoch); - assert_eq!(keplerian.center(), 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); - assert_float_eq!(keplerian.velocity().x, vel.x, rel <= 1e-6); - assert_float_eq!(keplerian.velocity().y, vel.y, rel <= 1e-6); - assert_float_eq!(keplerian.velocity().z, vel.z, rel <= 1e-6); - assert_float_eq!(keplerian.semi_major(), semi_major, rel <= 1e-6); - assert_float_eq!(keplerian.eccentricity(), eccentricity, rel <= 1e-6); - assert_float_eq!(keplerian.inclination(), inclination, rel <= 1e-6); - assert_float_eq!(keplerian.ascending_node(), ascending_node, rel <= 1e-6); - assert_float_eq!(keplerian.periapsis_arg(), periapsis_arg, rel <= 1e-6); - assert_float_eq!(keplerian.true_anomaly(), true_anomaly, rel <= 1e-6); - - 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!(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, - rel <= 1e-6 - ); - assert_float_eq!( - keplerian.periapsis_arg, - keplerian1.periapsis_arg, - rel <= 1e-6 - ); - assert_float_eq!(keplerian.true_anomaly, keplerian1.true_anomaly, rel <= 1e-6); - } -} diff --git a/crates/lox_core/src/two_body/elements.rs b/crates/lox_core/src/two_body/elements.rs deleted file mode 100644 index 5287a7bb..00000000 --- a/crates/lox_core/src/two_body/elements.rs +++ /dev/null @@ -1,421 +0,0 @@ -/* - * Copyright (c) 2023. 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 crate::math::normalize_two_pi; -use float_eq::float_eq; -use glam::{DMat3, DVec3}; -use std::f64::consts::PI; - -pub fn keplerian_to_perifocal( - semi_latus: f64, - eccentricity: f64, - true_anomaly: f64, - grav_param: f64, -) -> (DVec3, DVec3) { - let (sin_nu, cos_nu) = true_anomaly.sin_cos(); - let sqrt_mu_p = (grav_param / semi_latus).sqrt(); - - let pos = DVec3::new(cos_nu, sin_nu, 0.0) * (semi_latus / (1.0 + eccentricity * cos_nu)); - let vel = DVec3::new(-sin_nu, eccentricity + cos_nu, 0.0) * sqrt_mu_p; - - (pos, vel) -} - -pub fn keplerian_to_cartesian( - grav_param: f64, - semi_major: f64, - eccentricity: f64, - inclination: f64, - ascending_node: f64, - periapsis_arg: f64, - true_anomaly: f64, -) -> (DVec3, DVec3) { - let semi_latus = semi_latus_rectum(semi_major, eccentricity); - let (pos, vel) = keplerian_to_perifocal(semi_latus, eccentricity, true_anomaly, grav_param); - let rot = DMat3::from_rotation_z(ascending_node) - * DMat3::from_rotation_x(inclination) - * DMat3::from_rotation_z(periapsis_arg); - (rot * pos, rot * vel) -} - -pub fn semi_latus_rectum(semi_major: f64, eccentricity: f64) -> f64 { - if float_eq!(eccentricity, 0.0, abs <= 1e-3) { - semi_major - } else { - semi_major * (1.0 - eccentricity.powi(2)) - } -} - -type Keplerian = (f64, f64, f64, f64, f64, f64); - -pub fn cartesian_to_keplerian(grav_param: f64, pos: DVec3, vel: DVec3) -> Keplerian { - let r = pos.length(); - let v = vel.length(); - let h = pos.cross(vel); - let hm = h.length(); - let node = DVec3::Z.cross(h); - let e = eccentricity_vector(grav_param, pos, vel); - let eccentricity = e.length(); - let inclination = h.angle_between(DVec3::Z); - - let equatorial = is_equatorial(inclination); - let circular = is_circular(eccentricity); - - let semi_major = if circular { - hm.powi(2) / grav_param - } else { - -grav_param / (2.0 * (v.powi(2) / 2.0 - grav_param / r)) - }; - - let ascending_node; - let periapsis_arg; - let true_anomaly; - if equatorial && !circular { - ascending_node = 0.0; - periapsis_arg = azimuth(e); - true_anomaly = (h.dot(e.cross(pos)) / hm).atan2(pos.dot(e)); - } else if !equatorial && circular { - ascending_node = azimuth(node); - periapsis_arg = 0.0; - true_anomaly = (pos.dot(h.cross(node)) / hm).atan2(pos.dot(node)); - } else if equatorial && circular { - ascending_node = 0.0; - periapsis_arg = 0.0; - true_anomaly = azimuth(pos); - } else { - if semi_major > 0.0 { - let e_se = pos.dot(vel) / (grav_param * semi_major).sqrt(); - let e_ce = (r * v.powi(2)) / grav_param - 1.0; - true_anomaly = eccentric_to_true(e_se.atan2(e_ce), eccentricity); - } else { - let e_sh = pos.dot(vel) / (-grav_param * semi_major).sqrt(); - let e_ch = (r * v.powi(2)) / grav_param - 1.0; - true_anomaly = - hyperbolic_to_true(((e_ch + e_sh) / (e_ch - e_sh)).ln() / 2.0, eccentricity); - } - ascending_node = azimuth(node); - let px = pos.dot(node); - let py = pos.dot(h.cross(node)) / hm; - periapsis_arg = py.atan2(px) - true_anomaly; - } - - ( - semi_major, - eccentricity, - inclination, - mod_two_pi(ascending_node), - mod_two_pi(periapsis_arg), - normalize_two_pi(true_anomaly, 0.0), - ) -} - -fn mod_two_pi(a: f64) -> f64 { - let w = a % (2.0 * PI); - if w < 0.0 { - w + 2.0 * PI - } else { - w - } -} - -fn hyperbolic_to_true(e: f64, ecc: f64) -> f64 { - 2.0 * (((1.0 + ecc) / (ecc - 1.0)).sqrt() * (e / 2.0).tanh()).atan() -} - -fn eccentric_to_true(e: f64, ecc: f64) -> f64 { - 2.0 * (((1.0 + ecc) / (1.0 - ecc)).sqrt() * (e / 2.0).tan()).atan() -} - -fn azimuth(v: DVec3) -> f64 { - v.y.atan2(v.x) -} - -fn eccentricity_vector(grav_param: f64, pos: DVec3, vel: DVec3) -> DVec3 { - (pos * (vel.dot(vel) - grav_param / pos.length()) - vel * pos.dot(vel)) / grav_param -} - -fn is_equatorial(inclination: f64) -> bool { - float_eq!(inclination.abs(), 0.0, abs <= 1e-3) -} - -fn is_circular(eccentricity: f64) -> bool { - float_eq!(eccentricity, 0.0, abs <= 1e-3) -} - -#[cfg(test)] -mod tests { - use super::*; - use float_eq::assert_float_eq; - use glam::DVec3; - - #[test] - fn test_perifocal() { - let semi_latus = 1.13880762905224e7; - let eccentricity = 0.7311; - let true_anomaly = 0.44369564302687126; - let grav_param = 3.9860047e14; - - let pos = DVec3::new(6194863.12535486, 2944437.90016286, 0.0); - let vel = DVec3::new(-2539.71254827, 9668.69568539, 0.0); - let (pos1, vel1) = - keplerian_to_perifocal(semi_latus, eccentricity, true_anomaly, grav_param); - - assert_float_eq!(pos.x, pos1.x, rel <= 1e-8); - assert_float_eq!(pos.y, pos1.y, rel <= 1e-8); - assert_float_eq!(pos.z, pos1.z, rel <= 1e-8); - assert_float_eq!(vel.x, vel1.x, rel <= 1e-8); - assert_float_eq!(vel.y, vel1.y, rel <= 1e-8); - assert_float_eq!(vel.z, vel1.z, rel <= 1e-8); - } - - #[test] - fn test_elliptic() { - let grav_param = 3.9860047e14; - let semi_major = 24464560.0; - let eccentricity = 0.7311; - let inclination = 0.122138; - let ascending_node = 1.00681; - let periapsis_arg = 3.10686; - let true_anomaly = 0.44369564302687126; - let pos = DVec3::new( - -0.107622532467967e7, - -0.676589636432773e7, - -0.332308783350379e6, - ); - let vel = DVec3::new( - 0.935685775154103e4, - -0.331234775037644e4, - -0.118801577532701e4, - ); - - let ( - semi_major1, - eccentricity1, - inclination1, - ascending_node1, - periapsis_arg1, - true_anomaly1, - ) = cartesian_to_keplerian(grav_param, pos, vel); - let (pos1, vel1) = keplerian_to_cartesian( - grav_param, - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, - ); - - assert_float_eq!(pos.x, pos1.x, rel <= 1e-8); - assert_float_eq!(pos.y, pos1.y, rel <= 1e-8); - assert_float_eq!(pos.z, pos1.z, rel <= 1e-8); - assert_float_eq!(vel.x, vel1.x, rel <= 1e-8); - assert_float_eq!(vel.y, vel1.y, rel <= 1e-8); - assert_float_eq!(vel.z, vel1.z, rel <= 1e-8); - - assert_float_eq!(semi_major, semi_major1, rel <= 1e-8); - assert_float_eq!(eccentricity, eccentricity1, rel <= 1e-8); - assert_float_eq!(inclination, inclination1, rel <= 1e-8); - assert_float_eq!(ascending_node, ascending_node1, rel <= 1e-8); - assert_float_eq!(periapsis_arg, periapsis_arg1, rel <= 1e-8); - assert_float_eq!(true_anomaly, true_anomaly1, rel <= 1e-8); - } - - #[test] - fn test_circular() { - let grav_param = 3.986004418e14; - let semi_major = 6778136.6; - let eccentricity = 0.0; - let inclination = 15f64.to_radians(); - let ascending_node = 20f64.to_radians(); - let periapsis_arg = 0.0; - let true_anomaly = 30f64.to_radians(); - let pos = DVec3::new(4396398.60746266, 5083838.45333733, 877155.42119322); - let vel = DVec3::new(-5797.06004014, 4716.60916063, 1718.86034246); - - let ( - semi_major1, - eccentricity1, - inclination1, - ascending_node1, - periapsis_arg1, - true_anomaly1, - ) = cartesian_to_keplerian(grav_param, pos, vel); - let (pos1, vel1) = keplerian_to_cartesian( - grav_param, - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, - ); - - assert_float_eq!(pos.x, pos1.x, rel <= 1e-8); - assert_float_eq!(pos.y, pos1.y, rel <= 1e-8); - assert_float_eq!(pos.z, pos1.z, rel <= 1e-8); - assert_float_eq!(vel.x, vel1.x, rel <= 1e-4); - assert_float_eq!(vel.y, vel1.y, rel <= 1e-4); - assert_float_eq!(vel.z, vel1.z, rel <= 1e-4); - - assert_float_eq!(semi_major, semi_major1, rel <= 1e-8); - assert_float_eq!(eccentricity, eccentricity1, abs <= 1e-8); - assert_float_eq!(inclination, inclination1, rel <= 1e-8); - assert_float_eq!(ascending_node, ascending_node1, rel <= 1e-8); - assert_float_eq!(periapsis_arg, periapsis_arg1, rel <= 1e-8); - assert_float_eq!(true_anomaly, true_anomaly1, rel <= 1e-8); - } - - #[test] - fn test_circular_orekit() { - let grav_param = 3.9860047e14; - let semi_major = 24464560.0; - let eccentricity = 0.0; - let inclination = 0.122138; - let ascending_node = 1.00681; - let periapsis_arg = 0.0; - let true_anomaly = 0.048363; - - let (pos, vel) = keplerian_to_cartesian( - grav_param, - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, - ); - - let ( - semi_major1, - eccentricity1, - inclination1, - ascending_node1, - periapsis_arg1, - true_anomaly1, - ) = cartesian_to_keplerian(grav_param, pos, vel); - - assert_float_eq!(semi_major, semi_major1, rel <= 1e-8); - assert_float_eq!(eccentricity, eccentricity1, abs <= 1e-8); - assert_float_eq!(inclination, inclination1, rel <= 1e-8); - assert_float_eq!(ascending_node, ascending_node1, rel <= 1e-8); - assert_float_eq!(periapsis_arg, periapsis_arg1, rel <= 1e-8); - assert_float_eq!(true_anomaly, true_anomaly1, rel <= 1e-8); - } - - #[test] - fn test_hyperbolic_orekit() { - let grav_param = 3.9860047e14; - let semi_major = -24464560.0; - let eccentricity = 1.7311; - let inclination = 0.122138; - let ascending_node = 1.00681; - let periapsis_arg = 3.10686; - let true_anomaly = 0.12741601769795755; - - let (pos, vel) = keplerian_to_cartesian( - grav_param, - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, - ); - - let ( - semi_major1, - eccentricity1, - inclination1, - ascending_node1, - periapsis_arg1, - true_anomaly1, - ) = cartesian_to_keplerian(grav_param, pos, vel); - - assert_float_eq!(semi_major, semi_major1, rel <= 1e-8); - assert_float_eq!(eccentricity, eccentricity1, abs <= 1e-8); - assert_float_eq!(inclination, inclination1, rel <= 1e-8); - assert_float_eq!(ascending_node, ascending_node1, rel <= 1e-8); - assert_float_eq!(periapsis_arg, periapsis_arg1, rel <= 1e-8); - assert_float_eq!(true_anomaly, true_anomaly1, rel <= 1e-8); - } - - #[test] - fn test_equatorial() { - let grav_param = 3.9860047e14; - let semi_major = 24464560.0; - let eccentricity = 0.7311; - let inclination = 0.0; - let ascending_node = 0.0; - let periapsis_arg = 3.10686; - let true_anomaly = 0.44369564302687126; - - let (pos, vel) = keplerian_to_cartesian( - grav_param, - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, - ); - - let ( - semi_major1, - eccentricity1, - inclination1, - ascending_node1, - periapsis_arg1, - true_anomaly1, - ) = cartesian_to_keplerian(grav_param, pos, vel); - - assert_float_eq!(semi_major, semi_major1, rel <= 1e-8); - assert_float_eq!(eccentricity, eccentricity1, abs <= 1e-8); - assert_float_eq!(inclination, inclination1, rel <= 1e-8); - assert_float_eq!(ascending_node, ascending_node1, rel <= 1e-8); - assert_float_eq!(periapsis_arg, periapsis_arg1, rel <= 1e-8); - assert_float_eq!(true_anomaly, true_anomaly1, rel <= 1e-8); - } - - #[test] - fn test_circular_equatorial() { - let grav_param = 3.9860047e14; - let semi_major = 24464560.0; - let eccentricity = 0.0; - let inclination = 0.0; - let ascending_node = 0.0; - let periapsis_arg = 0.0; - let true_anomaly = 0.44369564302687126; - - let (pos, vel) = keplerian_to_cartesian( - grav_param, - semi_major, - eccentricity, - inclination, - ascending_node, - periapsis_arg, - true_anomaly, - ); - - let ( - semi_major1, - eccentricity1, - inclination1, - ascending_node1, - periapsis_arg1, - true_anomaly1, - ) = cartesian_to_keplerian(grav_param, pos, vel); - - assert_float_eq!(semi_major, semi_major1, rel <= 1e-8); - assert_float_eq!(eccentricity, eccentricity1, abs <= 1e-8); - assert_float_eq!(inclination, inclination1, rel <= 1e-8); - assert_float_eq!(ascending_node, ascending_node1, rel <= 1e-8); - assert_float_eq!(periapsis_arg, periapsis_arg1, rel <= 1e-8); - assert_float_eq!(true_anomaly, true_anomaly1, rel <= 1e-8); - } -} diff --git a/crates/lox_py/Cargo.toml b/crates/lox_py/Cargo.toml index 87d8ab13..007dbd6f 100644 --- a/crates/lox_py/Cargo.toml +++ b/crates/lox_py/Cargo.toml @@ -12,6 +12,11 @@ name = "lox_space" crate-type = ["cdylib"] [dependencies] -pyo3 = "0.20.0" +pyo3 = "0.20.2" lox_core.workspace = true +float_eq.workspace = true thiserror.workspace = true + +[dev-dependencies] +pyo3 = { version = "0.20.2", features = ["auto-initialize"] } +rstest.workspace = true diff --git a/crates/lox_py/lox_space.pyi b/crates/lox_py/lox_space.pyi index 135b64e2..cc6386a8 100644 --- a/crates/lox_py/lox_space.pyi +++ b/crates/lox_py/lox_space.pyi @@ -1,72 +1,105 @@ +class Epoch: + def __init__( + self, + scale: str, + year: int, + month: int, + day: int, + hour: int = 0, + minute: int = 0, + second: int = 0, + milli: int = 0, + micro: int = 0, + nano: int = 0, + pico: int = 0, + femto: int = 0, + atto: int = 0, + ): ... + class Sun: def __init__(self): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - def mean_radius(self) -> float: ... - def polar_radius(self) -> float: ... - def equatorial_radius(self) -> float: ... class Barycenter: def __init__(self, name: str): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... class Planet: def __init__(self, name: str): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - def mean_radius(self) -> float: ... - def polar_radius(self) -> float: ... - def equatorial_radius(self) -> float: ... class Satellite: def __init__(self, name: str): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - def mean_radius(self) -> float: ... - def polar_radius(self) -> float: ... - def subplanetary_radius(self) -> float: ... - def along_orbit_radius(self) -> float: ... class MinorBody: def __init__(self, name: str): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - def mean_radius(self) -> float: ... - def polar_radius(self) -> float: ... - def subplanetary_radius(self) -> float: ... - def along_orbit_radius(self) -> float: ... + +class Cartesian: + def __init__( + self, + time: Epoch, + body: Sun | Barycenter | Planet | Satellite | MinorBody, + frame: str, + x: float, + y: float, + z: float, + vx: float, + vy: float, + vz: float, + ): ... + def time(self) -> Epoch: ... + def reference_frame(self) -> str: ... + def origin(self) -> Sun | Barycenter | Planet | Satellite | MinorBody: ... + def position(self) -> tuple[float, float, float]: ... + def velocity(self) -> tuple[float, float, float]: ... + def to_keplerian(self) -> Keplerian: ... + +class Keplerian: + def __init__( + self, + time: Epoch, + body: Sun | Barycenter | Planet | Satellite | MinorBody, + frame: str, + semi_major_axis: float, + eccentricity: float, + inclination: float, + ascending_node: float, + periapsis_argument: float, + true_anomaly: float, + ): ... + def time(self) -> Epoch: ... + def reference_frame(self) -> str: ... + def origin(self) -> Sun | Barycenter | Planet | Satellite | MinorBody: ... + def semi_major_axis(self) -> float: ... + def eccentricity(self) -> float: ... + def inclination(self) -> float: ... + def ascending_node(self) -> float: ... + def periapsis_argument(self) -> float: ... + def true_anomaly(self) -> float: ... + def to_cartesian(self) -> Cartesian: ... diff --git a/crates/lox_py/pyproject.toml b/crates/lox_py/pyproject.toml index dfaec801..698a7e54 100644 --- a/crates/lox_py/pyproject.toml +++ b/crates/lox_py/pyproject.toml @@ -11,6 +11,15 @@ classifiers = [ "Programming Language :: Python :: Implementation :: PyPy", ] dynamic = ["version"] +dependencies = [ + "numpy~=1.26" +] + +[project.optional-dependencies] +dev = [ + "black", + "pytest", +] [tool.maturin] features = ["pyo3/extension-module"] diff --git a/crates/lox_py/src/bodies.rs b/crates/lox_py/src/bodies.rs index 01fae1a8..3ae8cee1 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,56 +14,58 @@ use lox_core::bodies::*; use crate::LoxPyError; #[pyclass(name = "Sun")] +#[derive(Clone, Default)] pub struct PySun; #[pymethods] impl PySun { #[new] - fn new() -> Self { + pub fn new() -> Self { Self } - fn __repr__(&self) -> PyResult { - Ok("Sun()".to_string()) + fn __repr__(&self) -> &str { + "Sun()" } - fn __str__(&self) -> PyResult { - Ok("Sun".to_string()) + fn __str__(&self) -> &str { + "Sun" } - fn id(&self) -> i32 { + pub fn id(&self) -> i32 { Sun.id().0 } - fn name(&self) -> &'static str { + pub fn name(&self) -> &'static str { Sun.name() } - fn gravitational_parameter(&self) -> f64 { + pub fn gravitational_parameter(&self) -> f64 { Sun.gravitational_parameter() } - fn mean_radius(&self) -> f64 { + pub fn mean_radius(&self) -> f64 { Sun.mean_radius() } - fn polar_radius(&self) -> f64 { + pub fn polar_radius(&self) -> f64 { Sun.polar_radius() } - fn equatorial_radius(&self) -> f64 { + pub fn equatorial_radius(&self) -> f64 { Sun.equatorial_radius() } } #[pyclass(name = "Barycenter")] -pub struct PyBarycenter(Box); +#[derive(Clone)] +pub struct PyBarycenter(Box); #[pymethods] impl PyBarycenter { #[new] - fn new(name: &str) -> Result { - let barycenter: Option> = match name { + pub fn new(name: &str) -> Result { + let barycenter: Option> = match name { "ssb" | "SSB" | "solar system barycenter" | "Solar System Barycenter" => { Some(Box::new(SolarSystemBarycenter)) } @@ -83,34 +86,35 @@ impl PyBarycenter { } } - fn __repr__(&self) -> PyResult { - Ok(format!("Barycenter(\"{}\")", self.name())) + fn __repr__(&self) -> String { + format!("Barycenter(\"{}\")", self.name()) } - fn __str__(&self) -> PyResult { - Ok(self.name().to_string()) + fn __str__(&self) -> &str { + self.name() } - fn id(&self) -> i32 { + pub fn id(&self) -> i32 { self.0.id().0 } - fn name(&self) -> &'static str { + pub fn name(&self) -> &'static str { self.0.name() } - fn gravitational_parameter(&self) -> f64 { + pub fn gravitational_parameter(&self) -> f64 { self.0.gravitational_parameter() } } #[pyclass(name = "Planet")] +#[derive(Clone)] pub struct PyPlanet(Box); #[pymethods] impl PyPlanet { #[new] - fn new(name: &str) -> Result { + pub fn new(name: &str) -> Result { let planet: Option> = match name { "mercury" | "Mercury" => Some(Box::new(Mercury)), "venus" | "Venus" => Some(Box::new(Venus)), @@ -129,46 +133,47 @@ impl PyPlanet { } } - fn __repr__(&self) -> PyResult { - Ok(format!("Planet(\"{}\")", self.name())) + fn __repr__(&self) -> String { + format!("Planet(\"{}\")", self.name()) } - fn __str__(&self) -> PyResult { - Ok(self.name().to_string()) + fn __str__(&self) -> &str { + self.name() } - fn id(&self) -> i32 { + pub fn id(&self) -> i32 { self.0.id().0 } - fn name(&self) -> &'static str { + pub fn name(&self) -> &'static str { self.0.name() } - fn gravitational_parameter(&self) -> f64 { + pub fn gravitational_parameter(&self) -> f64 { self.0.gravitational_parameter() } - fn mean_radius(&self) -> f64 { + pub fn mean_radius(&self) -> f64 { self.0.mean_radius() } - fn polar_radius(&self) -> f64 { + pub fn polar_radius(&self) -> f64 { self.0.polar_radius() } - fn equatorial_radius(&self) -> f64 { + pub fn equatorial_radius(&self) -> f64 { self.0.equatorial_radius() } } #[pyclass(name = "Satellite")] +#[derive(Clone)] pub struct PySatellite(Box); #[pymethods] impl PySatellite { #[new] - fn new(name: &str) -> Result { + pub fn new(name: &str) -> Result { let satellite: Option> = match name { "moon" | "Moon" | "luna" | "Luna" => Some(Box::new(Moon)), "phobos" | "Phobos" => Some(Box::new(Phobos)), @@ -218,50 +223,51 @@ impl PySatellite { } } - fn __repr__(&self) -> PyResult { - Ok(format!("Satellite(\"{}\")", self.name())) + fn __repr__(&self) -> String { + format!("Satellite(\"{}\")", self.name()) } - fn __str__(&self) -> PyResult { - Ok(self.name().to_string()) + fn __str__(&self) -> &str { + self.name() } - fn id(&self) -> i32 { + pub fn id(&self) -> i32 { self.0.id().0 } - fn name(&self) -> &'static str { + pub fn name(&self) -> &'static str { self.0.name() } - fn gravitational_parameter(&self) -> f64 { + pub fn gravitational_parameter(&self) -> f64 { self.0.gravitational_parameter() } - fn mean_radius(&self) -> f64 { + pub fn mean_radius(&self) -> f64 { self.0.mean_radius() } - fn polar_radius(&self) -> f64 { + pub fn polar_radius(&self) -> f64 { self.0.polar_radius() } - fn subplanetary_radius(&self) -> f64 { + pub fn subplanetary_radius(&self) -> f64 { self.0.subplanetary_radius() } - fn along_orbit_radius(&self) -> f64 { + pub fn along_orbit_radius(&self) -> f64 { self.0.along_orbit_radius() } } #[pyclass(name = "MinorBody")] +#[derive(Clone)] pub struct PyMinorBody(Box); #[pymethods] impl PyMinorBody { #[new] - fn new(name: &str) -> Result { + pub fn new(name: &str) -> Result { let minor: Option> = match name { "ceres" | "Ceres" => Some(Box::new(Ceres)), "vesta" | "Vesta" => Some(Box::new(Vesta)), @@ -276,39 +282,402 @@ impl PyMinorBody { } } - fn __repr__(&self) -> PyResult { - Ok(format!("MinorBody(\"{}\")", self.name())) + fn __repr__(&self) -> String { + format!("MinorBody(\"{}\")", self.name()) } - fn __str__(&self) -> PyResult { - Ok(self.name().to_string()) + fn __str__(&self) -> &str { + self.name() } - fn id(&self) -> i32 { + pub fn id(&self) -> i32 { self.0.id().0 } - fn name(&self) -> &'static str { + pub fn name(&self) -> &'static str { self.0.name() } - fn gravitational_parameter(&self) -> f64 { + pub fn gravitational_parameter(&self) -> f64 { self.0.gravitational_parameter() } - fn mean_radius(&self) -> f64 { + pub fn mean_radius(&self) -> f64 { self.0.mean_radius() } - fn polar_radius(&self) -> f64 { + pub fn polar_radius(&self) -> f64 { self.0.polar_radius() } - fn subplanetary_radius(&self) -> f64 { + pub fn subplanetary_radius(&self) -> f64 { self.0.subplanetary_radius() } - fn along_orbit_radius(&self) -> f64 { + pub fn along_orbit_radius(&self) -> f64 { self.0.along_orbit_radius() } } + +#[derive(Clone)] +pub enum PyBody { + Barycenter(PyBarycenter), + Sun(PySun), + Planet(PyPlanet), + 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(), + } + } +} + +#[cfg(test)] +mod tests { + use rstest::rstest; + + use super::*; + + #[test] + fn test_sun() { + let sun = PySun::new(); + assert_eq!(sun.__repr__(), "Sun()"); + assert_eq!(sun.__str__(), "Sun"); + assert_eq!(sun.id(), Sun.id().0); + assert_eq!(sun.name(), Sun.name()); + assert_eq!(sun.gravitational_parameter(), Sun.gravitational_parameter()); + assert_eq!(sun.mean_radius(), Sun.mean_radius()); + assert_eq!(sun.polar_radius(), Sun.polar_radius()); + assert_eq!(sun.equatorial_radius(), Sun.equatorial_radius()); + } + + #[rstest] + #[case("Solar System Barycenter", SolarSystemBarycenter)] + #[case("SSB", SolarSystemBarycenter)] + #[case("Mercury Barycenter", MercuryBarycenter)] + #[case("Venus Barycenter", VenusBarycenter)] + #[case("Earth Barycenter", EarthBarycenter)] + #[case("Mars Barycenter", MarsBarycenter)] + #[case("Jupiter Barycenter", JupiterBarycenter)] + #[case("Saturn Barycenter", SaturnBarycenter)] + #[case("Uranus Barycenter", UranusBarycenter)] + #[case("Neptune Barycenter", NeptuneBarycenter)] + #[case("Pluto Barycenter", PlutoBarycenter)] + fn test_barycenter(#[case] name: &str, #[case] barycenter: impl Barycenter) { + let py_barycenter = PyBarycenter::new(name).expect("barycenter should be valid"); + assert_eq!( + py_barycenter.__repr__(), + format!("Barycenter(\"{}\")", barycenter.name()) + ); + assert_eq!(py_barycenter.__str__(), barycenter.name()); + assert_eq!(py_barycenter.name(), barycenter.name()); + let py_barycenter = + PyBarycenter::new(&name.to_lowercase()).expect("barycenter should be valid"); + assert_eq!( + py_barycenter.__repr__(), + format!("Barycenter(\"{}\")", barycenter.name()) + ); + assert_eq!(py_barycenter.__str__(), barycenter.name()); + assert_eq!(py_barycenter.name(), barycenter.name()); + assert_eq!(py_barycenter.id(), barycenter.id().0); + assert_eq!( + py_barycenter.gravitational_parameter(), + barycenter.gravitational_parameter() + ); + } + + #[test] + fn test_invalid_barycenter() { + let barycenter = PyBarycenter::new("Rupert Barycenter"); + assert!(barycenter.is_err()); + } + + #[rstest] + #[case("Mercury", Mercury)] + #[case("Venus", Venus)] + #[case("Earth", Earth)] + #[case("Mars", Mars)] + #[case("Jupiter", Jupiter)] + #[case("Saturn", Saturn)] + #[case("Uranus", Uranus)] + #[case("Neptune", Neptune)] + #[case("Pluto", Pluto)] + fn test_planet(#[case] name: &str, #[case] planet: impl Planet) { + let py_planet = PyPlanet::new(name).expect("planet should be valid"); + assert_eq!(py_planet.__repr__(), format!("Planet(\"{}\")", name)); + assert_eq!(py_planet.__str__(), name); + assert_eq!(py_planet.name(), name); + let py_planet = PyPlanet::new(&name.to_lowercase()).expect("planet should be valid"); + assert_eq!(py_planet.__repr__(), format!("Planet(\"{}\")", name)); + assert_eq!(py_planet.__str__(), name); + assert_eq!(py_planet.name(), name); + assert_eq!(py_planet.id(), planet.id().0); + assert_eq!( + py_planet.gravitational_parameter(), + planet.gravitational_parameter() + ); + assert_eq!(py_planet.mean_radius(), planet.mean_radius()); + assert_eq!(py_planet.polar_radius(), planet.polar_radius()); + assert_eq!(py_planet.equatorial_radius(), planet.equatorial_radius()); + } + + #[test] + fn test_invalid_planet() { + let planet = PyPlanet::new("Rupert"); + assert!(planet.is_err()); + } + + #[rstest] + #[case("Moon", Moon)] + #[case("Luna", Moon)] + #[case("Phobos", Phobos)] + #[case("Deimos", Deimos)] + #[case("Io", Io)] + #[case("Europa", Europa)] + #[case("Ganymede", Ganymede)] + #[case("Callisto", Callisto)] + #[case("Amalthea", Amalthea)] + #[case("Himalia", Himalia)] + #[case("Thebe", Thebe)] + #[case("Adrastea", Adrastea)] + #[case("Metis", Metis)] + #[case("Mimas", Mimas)] + #[case("Enceladus", Enceladus)] + #[case("Tethys", Tethys)] + #[case("Dione", Dione)] + #[case("Rhea", Rhea)] + #[case("Titan", Titan)] + #[case("Hyperion", Hyperion)] + #[case("Iapetus", Iapetus)] + #[case("Phoebe", Phoebe)] + #[case("Janus", Janus)] + #[case("Epimetheus", Epimetheus)] + #[case("Helene", Helene)] + #[case("Atlas", Atlas)] + #[case("Prometheus", Prometheus)] + #[case("Pandora", Pandora)] + #[case("Ariel", Ariel)] + #[case("Umbriel", Umbriel)] + #[case("Titania", Titania)] + #[case("Oberon", Oberon)] + #[case("Miranda", Miranda)] + #[case("Triton", Triton)] + #[case("Naiad", Naiad)] + #[case("Thalassa", Thalassa)] + #[case("Despina", Despina)] + #[case("Galatea", Galatea)] + #[case("Larissa", Larissa)] + #[case("Proteus", Proteus)] + #[case("Charon", Charon)] + fn test_satellite(#[case] name: &str, #[case] satellite: impl Satellite) { + let py_satellite = PySatellite::new(name).expect("satellite should be valid"); + assert_eq!( + py_satellite.__repr__(), + format!("Satellite(\"{}\")", satellite.name()) + ); + assert_eq!(py_satellite.__str__(), satellite.name()); + assert_eq!(py_satellite.name(), satellite.name()); + let py_satellite = + PySatellite::new(&name.to_lowercase()).expect("satellite should be valid"); + assert_eq!( + py_satellite.__repr__(), + format!("Satellite(\"{}\")", satellite.name()) + ); + assert_eq!(py_satellite.__str__(), satellite.name()); + assert_eq!(py_satellite.name(), satellite.name()); + assert_eq!(py_satellite.id(), satellite.id().0); + assert_eq!( + py_satellite.gravitational_parameter(), + satellite.gravitational_parameter() + ); + assert_eq!(py_satellite.mean_radius(), satellite.mean_radius()); + assert_eq!(py_satellite.polar_radius(), satellite.polar_radius()); + assert_eq!( + py_satellite.subplanetary_radius(), + satellite.subplanetary_radius() + ); + assert_eq!( + py_satellite.along_orbit_radius(), + satellite.along_orbit_radius() + ); + } + + #[test] + fn test_invalid_satellite() { + let satellite = PySatellite::new("Endor"); + assert!(satellite.is_err()); + } + + #[rstest] + #[case("Ceres", Ceres)] + #[case("Vesta", Vesta)] + #[case("Psyche", Psyche)] + #[case("Eros", Eros)] + #[case("Davida", Davida)] + fn test_minor_body(#[case] name: &str, #[case] minor_body: impl MinorBody) { + let py_minor_body = PyMinorBody::new(name).expect("minor body should be valid"); + assert_eq!( + py_minor_body.__repr__(), + format!("MinorBody(\"{}\")", minor_body.name()) + ); + assert_eq!(py_minor_body.__str__(), minor_body.name()); + assert_eq!(py_minor_body.name(), minor_body.name()); + let py_minor_body = + PyMinorBody::new(&name.to_lowercase()).expect("minor body should be valid"); + assert_eq!( + py_minor_body.__repr__(), + format!("MinorBody(\"{}\")", minor_body.name()) + ); + assert_eq!(py_minor_body.__str__(), minor_body.name()); + assert_eq!(py_minor_body.name(), minor_body.name()); + assert_eq!(py_minor_body.id(), minor_body.id().0); + assert_eq!( + py_minor_body.gravitational_parameter(), + minor_body.gravitational_parameter() + ); + assert_eq!(py_minor_body.mean_radius(), minor_body.mean_radius()); + assert_eq!(py_minor_body.polar_radius(), minor_body.polar_radius()); + assert_eq!( + py_minor_body.subplanetary_radius(), + minor_body.subplanetary_radius() + ); + assert_eq!( + py_minor_body.along_orbit_radius(), + minor_body.along_orbit_radius() + ); + } + + #[test] + fn test_invalid_minor_body() { + let minor_body = PyMinorBody::new("Bielefeld"); + assert!(minor_body.is_err()); + } + + #[test] + fn test_body() { + let sun = PyBody::Sun(PySun::new()); + assert_eq!(sun.id(), Sun.id()); + assert_eq!(sun.name(), Sun.name()); + assert_eq!(sun.gravitational_parameter(), Sun.gravitational_parameter()); + let sun1: PyBody = PyObject::from(sun.clone()) + .try_into() + .expect("sun is valid"); + assert_eq!(sun1.id(), sun.id()); + + let barycenter = PyBody::Barycenter(PyBarycenter::new("ssb").expect("barycenter is valid")); + assert_eq!(barycenter.id(), SolarSystemBarycenter.id()); + assert_eq!(barycenter.name(), SolarSystemBarycenter.name()); + assert_eq!( + barycenter.gravitational_parameter(), + SolarSystemBarycenter.gravitational_parameter() + ); + let barycenter1: PyBody = PyObject::from(barycenter.clone()) + .try_into() + .expect("barycenter is valid"); + assert_eq!(barycenter1.id(), barycenter.id()); + + let planet = PyBody::Planet(PyPlanet::new("earth").expect("planet is valid")); + assert_eq!(planet.id(), Earth.id()); + assert_eq!(planet.name(), Earth.name()); + assert_eq!( + planet.gravitational_parameter(), + Earth.gravitational_parameter() + ); + let planet1: PyBody = PyObject::from(planet.clone()) + .try_into() + .expect("planet is valid"); + assert_eq!(planet1.id(), planet.id()); + + let satellite = PyBody::Satellite(PySatellite::new("moon").expect("satellite is valid")); + assert_eq!(satellite.id(), Moon.id()); + assert_eq!(satellite.name(), Moon.name()); + assert_eq!( + satellite.gravitational_parameter(), + Moon.gravitational_parameter() + ); + let satellite1: PyBody = PyObject::from(satellite.clone()) + .try_into() + .expect("satellite is valid"); + assert_eq!(satellite1.id(), satellite.id()); + + let minor_body = PyBody::MinorBody(PyMinorBody::new("ceres").expect("minor body is valid")); + assert_eq!(minor_body.id(), Ceres.id()); + assert_eq!(minor_body.name(), Ceres.name()); + assert_eq!( + minor_body.gravitational_parameter(), + Ceres.gravitational_parameter() + ); + let minor_body1: PyBody = PyObject::from(minor_body.clone()) + .try_into() + .expect("minor_body is valid"); + assert_eq!(minor_body1.id(), minor_body.id()); + + let obj = Python::with_gil(|py| 1.into_py(py)); + let body = PyBody::try_from(obj); + assert!(body.is_err()); + } +} diff --git a/crates/lox_py/src/coords.rs b/crates/lox_py/src/coords.rs new file mode 100644 index 00000000..a0928c36 --- /dev/null +++ b/crates/lox_py/src/coords.rs @@ -0,0 +1,328 @@ +/* + * 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::str::FromStr; + +use pyo3::prelude::*; + +use lox_core::bodies::PointMass; +use lox_core::coords::states::{CartesianState, KeplerianState, TwoBodyState}; +use lox_core::coords::DVec3; + +use crate::bodies::PyBody; +use crate::frames::PyFrame; +use crate::time::PyEpoch; + +#[pyclass(name = "Cartesian")] +pub struct PyCartesian { + state: CartesianState, + origin: PyBody, + frame: PyFrame, +} + +#[pymethods] +impl PyCartesian { + #[allow(clippy::too_many_arguments)] + #[new] + fn new( + time: &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)?; + let state = CartesianState::new(time.0, DVec3::new(x, y, z), DVec3::new(vx, vy, vz)); + Ok(Self { + state, + origin, + frame, + }) + } + + fn time(&self) -> PyEpoch { + PyEpoch(self.state.time()) + } + + fn reference_frame(&self) -> String { + format!("{}", self.frame) + } + + fn origin(&self) -> PyObject { + self.origin.clone().into() + } + + fn position(&self) -> (f64, f64, f64) { + let position = self.state.position(); + (position.x, position.y, position.z) + } + + fn velocity(&self) -> (f64, f64, f64) { + let velocity = self.state.velocity(); + (velocity.x, velocity.y, velocity.z) + } + + fn to_keplerian(&self) -> PyKeplerian { + let mu = self.origin.gravitational_parameter(); + let state = self.state.to_keplerian_state(mu); + PyKeplerian { + origin: self.origin.clone(), + frame: self.frame.clone(), + state, + } + } +} + +#[pyclass(name = "Keplerian")] +pub struct PyKeplerian { + state: KeplerianState, + origin: PyBody, + frame: PyFrame, +} + +#[pymethods] +impl PyKeplerian { + #[new] + #[allow(clippy::too_many_arguments)] + fn new( + t: &PyEpoch, + body: PyObject, + frame: &str, + semi_major_axis: f64, + eccentricity: f64, + inclination: f64, + ascending_node: f64, + periapsis_argument: f64, + true_anomaly: f64, + ) -> PyResult { + let origin: PyBody = body.try_into()?; + let frame = PyFrame::from_str(frame)?; + let state = KeplerianState::new( + t.0, + semi_major_axis, + eccentricity, + inclination, + ascending_node, + periapsis_argument, + true_anomaly, + ); + Ok(Self { + state, + origin, + frame, + }) + } + + fn time(&self) -> PyEpoch { + PyEpoch(self.state.time()) + } + + fn reference_frame(&self) -> String { + format!("{}", self.frame) + } + + fn origin(&self) -> PyObject { + self.origin.clone().into() + } + + fn semi_major_axis(&self) -> f64 { + self.state.semi_major_axis() + } + + 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_argument(&self) -> f64 { + self.state.periapsis_argument() + } + + fn true_anomaly(&self) -> f64 { + self.state.true_anomaly() + } + + fn to_cartesian(&self) -> PyCartesian { + let mu = self.origin.gravitational_parameter(); + let state = self.state.to_cartesian_state(mu); + PyCartesian { + state, + origin: self.origin.clone(), + frame: self.frame.clone(), + } + } +} + +#[cfg(test)] +mod tests { + use float_eq::assert_float_eq; + + use crate::bodies::PyPlanet; + + use super::*; + + #[test] + fn test_cartesian() { + let epoch = PyEpoch::new( + "TDB", + 2023, + 3, + 25, + Some(21), + Some(8), + Some(0), + Some(0), + Some(0), + Some(0), + Some(0), + Some(0), + Some(0), + ) + .expect("time should be valid"); + let body = Python::with_gil(|py| { + PyPlanet::new("Earth") + .expect("body should be valid") + .into_py(py) + }); + let pos = DVec3::new( + -0.107622532467967e7, + -0.676589636432773e7, + -0.332308783350379e6, + ) * 1e-3; + let vel = DVec3::new( + 0.935685775154103e4, + -0.331234775037644e4, + -0.118801577532701e4, + ) * 1e-3; + + let cartesian = PyCartesian::new( + &epoch, + body.clone(), + "ICRF", + pos.x, + pos.y, + pos.z, + vel.x, + vel.y, + vel.z, + ) + .expect("cartesian state should be valid"); + let cartesian1 = cartesian.to_keplerian().to_cartesian(); + + let origin = + Python::with_gil(|py| body.extract::(py)).expect("origin should be a planet"); + let origin1 = Python::with_gil(|py| cartesian1.origin().extract::(py)) + .expect("origin should be a planet"); + assert_eq!(cartesian1.time(), epoch); + assert_eq!(origin1.name(), origin.name()); + assert_eq!(cartesian1.reference_frame(), "ICRF"); + + assert_float_eq!(cartesian.position().0, cartesian1.position().0, rel <= 1e-8); + assert_float_eq!(cartesian.position().1, cartesian1.position().1, rel <= 1e-8); + assert_float_eq!(cartesian.position().2, cartesian1.position().2, rel <= 1e-8); + assert_float_eq!(cartesian.velocity().0, cartesian1.velocity().0, rel <= 1e-6); + assert_float_eq!(cartesian.velocity().1, cartesian1.velocity().1, rel <= 1e-6); + assert_float_eq!(cartesian.velocity().2, cartesian1.velocity().2, rel <= 1e-6); + } + + #[test] + fn test_keplerian() { + let epoch = PyEpoch::new( + "TDB", + 2023, + 3, + 25, + Some(21), + Some(8), + Some(0), + Some(0), + Some(0), + Some(0), + Some(0), + Some(0), + Some(0), + ) + .expect("time should be valid"); + let body = Python::with_gil(|py| { + PyPlanet::new("Earth") + .expect("body should be valid") + .into_py(py) + }); + let semi_major = 24464560.0e-3; + let eccentricity = 0.7311; + let inclination = 0.122138; + let ascending_node = 1.00681; + let periapsis_arg = 3.10686; + let true_anomaly = 0.44369564302687126; + + let keplerian = PyKeplerian::new( + &epoch, + body.clone(), + "ICRF", + semi_major, + eccentricity, + inclination, + ascending_node, + periapsis_arg, + true_anomaly, + ) + .expect("Keplerian state should be valid"); + let keplerian1 = keplerian.to_cartesian().to_keplerian(); + + let origin = + Python::with_gil(|py| body.extract::(py)).expect("origin should be a planet"); + let origin1 = Python::with_gil(|py| keplerian1.origin().extract::(py)) + .expect("origin should be a planet"); + assert_eq!(keplerian1.time(), epoch); + assert_eq!(origin1.name(), origin.name()); + assert_eq!(keplerian1.reference_frame(), "ICRF"); + + assert_float_eq!( + keplerian.semi_major_axis(), + keplerian1.semi_major_axis(), + rel <= 1e-6 + ); + 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(), + rel <= 1e-6 + ); + assert_float_eq!( + keplerian.periapsis_argument(), + keplerian1.periapsis_argument(), + rel <= 1e-6 + ); + assert_float_eq!( + keplerian.true_anomaly(), + keplerian1.true_anomaly(), + rel <= 1e-6 + ); + } +} diff --git a/crates/lox_py/src/frames.rs b/crates/lox_py/src/frames.rs new file mode 100644 index 00000000..0353399c --- /dev/null +++ b/crates/lox_py/src/frames.rs @@ -0,0 +1,101 @@ +/* + * 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, Eq, PartialEq)] +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())), + } + } +} + +#[cfg(test)] +mod tests { + use rstest::rstest; + + use super::*; + + #[rstest] + #[case("icrf", PyFrame::Icrf(Icrf))] + #[case("iau_mercury", PyFrame::IauMercury(BodyFixed(Mercury)))] + #[case("iau_venus", PyFrame::IauVenus(BodyFixed(Venus)))] + #[case("iau_earth", PyFrame::IauEarth(BodyFixed(Earth)))] + #[case("iau_mars", PyFrame::IauMars(BodyFixed(Mars)))] + #[case("iau_jupiter", PyFrame::IauJupiter(BodyFixed(Jupiter)))] + #[case("iau_saturn", PyFrame::IauSaturn(BodyFixed(Saturn)))] + #[case("iau_uranus", PyFrame::IauUranus(BodyFixed(Uranus)))] + #[case("iau_neptune", PyFrame::IauNeptune(BodyFixed(Neptune)))] + #[case("iau_pluto", PyFrame::IauPluto(BodyFixed(Pluto)))] + fn test_frames(#[case] name: &str, #[case] exp: PyFrame) { + let upper = name.to_uppercase(); + let act = PyFrame::from_str(name).expect("frame should be valid"); + assert_eq!(act, exp); + let act = PyFrame::from_str(&upper).expect("frame should be valid"); + assert_eq!(act, exp); + assert_eq!(format!("{}", act), upper); + } + + #[test] + fn test_invalid_frame() { + let frame = PyFrame::from_str("Flat Earth"); + assert!(frame.is_err()); + } +} diff --git a/crates/lox_py/src/lib.rs b/crates/lox_py/src/lib.rs index 4e6fa920..4fcea7e9 100644 --- a/crates/lox_py/src/lib.rs +++ b/crates/lox_py/src/lib.rs @@ -10,21 +10,24 @@ use pyo3::exceptions::PyValueError; use pyo3::prelude::*; use thiserror::Error; -use lox_core::errors::LoxError; -use lox_core::time::dates::{Date, Time}; -use lox_core::time::epochs::Epoch; -use lox_core::time::epochs::TimeScale; - use crate::bodies::{PyBarycenter, PyMinorBody, PyPlanet, PySatellite, PySun}; +use crate::coords::{PyCartesian, PyKeplerian}; +use crate::time::{PyEpoch, PyTimeScale}; +use lox_core::errors::LoxError; mod bodies; +mod coords; +mod frames; +mod time; #[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,115 +37,15 @@ 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, } } } -#[pyclass(name = "TimeScale")] -struct PyTimeScale(TimeScale); - -#[pymethods] -impl PyTimeScale { - #[new] - fn new(name: &str) -> Result { - match name { - "TAI" => Ok(PyTimeScale(TimeScale::TAI)), - "TCB" => Ok(PyTimeScale(TimeScale::TCB)), - "TCG" => Ok(PyTimeScale(TimeScale::TCG)), - "TDB" => Ok(PyTimeScale(TimeScale::TDB)), - "TT" => Ok(PyTimeScale(TimeScale::TT)), - "UT1" => Ok(PyTimeScale(TimeScale::UT1)), - _ => Err(LoxPyError::InvalidTimeScale(name.to_string())), - } - } - - fn __repr__(&self) -> String { - format!("TimeScale(\"{}\")", self.0) - } - - fn __str__(&self) -> String { - format!("{}", self.0) - } -} - -#[pyclass(name = "Epoch")] -struct PyEpoch(Epoch); - -#[pymethods] -impl PyEpoch { - #[allow(clippy::too_many_arguments)] - #[pyo3(signature = ( - scale, - year, - month, - day, - hour = 0, - minute = 0, - second = 0, - milli = 0, - micro = 0, - nano = 0, - pico = 0, - femto = 0, - atto = 0 - ))] - #[new] - fn new( - scale: &str, - year: i64, - month: i64, - day: i64, - hour: Option, - minute: Option, - second: Option, - milli: Option, - micro: Option, - nano: Option, - pico: Option, - femto: Option, - atto: Option, - ) -> Result { - let time_scale = PyTimeScale::new(scale)?; - let date = Date::new(year, month, day)?; - - let hour = hour.unwrap_or(0); - let minute = minute.unwrap_or(0); - let second = second.unwrap_or(0); - let mut time = Time::new(hour, minute, second)?; - if let Some(milli) = milli { - time = time.milli(milli); - } - if let Some(micro) = micro { - time = time.micro(micro); - } - if let Some(nano) = nano { - time = time.nano(nano); - } - if let Some(pico) = pico { - time = time.pico(pico); - } - if let Some(femto) = femto { - time = time.femto(femto); - } - if let Some(atto) = atto { - time = time.atto(atto); - } - Ok(PyEpoch(Epoch::from_date_and_time(time_scale.0, date, time))) - } - - fn attosecond(&self) -> i64 { - self.0.attosecond() - } - - fn __str__(&self) -> String { - "foo".to_string() - } -} - /// A Python module implemented in Rust. #[pymodule] fn lox_space(_py: Python, m: &PyModule) -> PyResult<()> { @@ -153,5 +56,7 @@ fn lox_space(_py: Python, m: &PyModule) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; + m.add_class::()?; Ok(()) } diff --git a/crates/lox_py/src/time.rs b/crates/lox_py/src/time.rs new file mode 100644 index 00000000..6d75272f --- /dev/null +++ b/crates/lox_py/src/time.rs @@ -0,0 +1,156 @@ +/* + * 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 pyo3::{pyclass, pymethods}; + +use lox_core::time::dates::{Date, Time}; +use lox_core::time::epochs::{Epoch, TimeScale}; + +use crate::LoxPyError; + +#[pyclass(name = "TimeScale")] +pub struct PyTimeScale(pub TimeScale); + +#[pymethods] +impl PyTimeScale { + #[new] + fn new(name: &str) -> Result { + match name { + "TAI" => Ok(PyTimeScale(TimeScale::TAI)), + "TCB" => Ok(PyTimeScale(TimeScale::TCB)), + "TCG" => Ok(PyTimeScale(TimeScale::TCG)), + "TDB" => Ok(PyTimeScale(TimeScale::TDB)), + "TT" => Ok(PyTimeScale(TimeScale::TT)), + "UT1" => Ok(PyTimeScale(TimeScale::UT1)), + _ => Err(LoxPyError::InvalidTimeScale(name.to_string())), + } + } + + fn __repr__(&self) -> String { + format!("TimeScale(\"{}\")", self.0) + } + + fn __str__(&self) -> String { + format!("{}", self.0) + } +} + +#[pyclass(name = "Epoch")] +#[derive(Debug, Copy, Clone, Eq, PartialEq)] +pub struct PyEpoch(pub Epoch); + +#[pymethods] +impl PyEpoch { + #[allow(clippy::too_many_arguments)] + #[pyo3(signature = ( + scale, + year, + month, + day, + hour = 0, + minute = 0, + second = 0, + milli = 0, + micro = 0, + nano = 0, + pico = 0, + femto = 0, + atto = 0 + ))] + #[new] + pub fn new( + scale: &str, + year: i64, + month: i64, + day: i64, + hour: Option, + minute: Option, + second: Option, + milli: Option, + micro: Option, + nano: Option, + pico: Option, + femto: Option, + atto: Option, + ) -> Result { + let time_scale = PyTimeScale::new(scale)?; + let date = Date::new(year, month, day)?; + + let hour = hour.unwrap_or(0); + let minute = minute.unwrap_or(0); + let second = second.unwrap_or(0); + let mut time = Time::new(hour, minute, second)?; + if let Some(milli) = milli { + time = time.milli(milli); + } + if let Some(micro) = micro { + time = time.micro(micro); + } + if let Some(nano) = nano { + time = time.nano(nano); + } + if let Some(pico) = pico { + time = time.pico(pico); + } + if let Some(femto) = femto { + time = time.femto(femto); + } + if let Some(atto) = atto { + time = time.atto(atto); + } + Ok(PyEpoch(Epoch::from_date_and_time(time_scale.0, date, time))) + } +} + +#[cfg(test)] +mod tests { + use rstest::rstest; + + use super::*; + + #[rstest] + #[case("TAI", TimeScale::TAI)] + #[case("TCB", TimeScale::TCB)] + #[case("TCG", TimeScale::TCG)] + #[case("TDB", TimeScale::TDB)] + #[case("TT", TimeScale::TT)] + #[case("UT1", TimeScale::UT1)] + fn test_scale(#[case] name: &str, #[case] scale: TimeScale) { + let py_scale = PyTimeScale::new(name).expect("time scale should be valid"); + assert_eq!(py_scale.0, scale); + assert_eq!(py_scale.__str__(), name); + assert_eq!(py_scale.__repr__(), format!("TimeScale(\"{}\")", name)); + } + + #[test] + fn test_invalid_scale() { + let py_scale = PyTimeScale::new("disco time"); + assert!(py_scale.is_err()) + } + + #[test] + fn test_time() { + let time = PyEpoch::new( + "TDB", + 2024, + 1, + 1, + Some(1), + Some(1), + Some(1), + Some(123), + Some(456), + Some(789), + Some(123), + Some(456), + Some(789), + ) + .expect("time should be valid"); + assert_eq!(time.0.attosecond(), 123456789123456789); + } +} diff --git a/crates/lox_py/tests/test_coords.py b/crates/lox_py/tests/test_coords.py new file mode 100644 index 00000000..33926b0d --- /dev/null +++ b/crates/lox_py/tests/test_coords.py @@ -0,0 +1,22 @@ +# 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/. + +import lox_space as lox + + +def test_coords(): + time = lox.Epoch("TDB", 2016, 5, 30, 12) + x = 6068.27927 + y = -1692.84394 + z = -2516.61918 + vx = -0.660415582 + vy = 5.495938726 + vz = -5.303093233 + iss_cartesian = lox.Cartesian( + time, lox.Planet("Earth"), "ICRF", x, y, z, vx, vy, vz + ) + iss = iss_cartesian.to_keplerian() + assert isinstance(iss, lox.Keplerian) diff --git a/tools/lox_gen/Cargo.lock b/tools/lox_gen/Cargo.lock index c235bf80..9290fd94 100644 --- a/tools/lox_gen/Cargo.lock +++ b/tools/lox_gen/Cargo.lock @@ -8,6 +8,12 @@ version = "1.1.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" +[[package]] +name = "dyn-clone" +version = "1.0.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "545b22097d44f8a9581187cdf93de7a71e4722bf51200cfaba810865b49a495d" + [[package]] name = "fast_polynomial" version = "0.1.0" @@ -25,9 +31,9 @@ checksum = "28a80e3145d8ad11ba0995949bbcf48b9df2be62772b3d351ef017dff6ecb853" [[package]] name = "glam" -version = "0.24.2" +version = "0.25.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b5418c17512bdf42730f9032c74e1ae39afc408745ebb2acf72fbc4691c17945" +checksum = "151665d9be52f9bb40fc7966565d39666f2d1e69233571b71b87791c7e0528b3" [[package]] name = "lazy_static" @@ -39,6 +45,7 @@ checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" name = "lox_core" version = "0.1.0" dependencies = [ + "dyn-clone", "fast_polynomial", "float_eq", "glam",