Integration of the photom crate

Cargo.toml

@@ -1,9 +1,9 @@
 [package]
 name = "outfit"
-version = "2.1.0"
+version = "3.0.0"
 edition = "2021"
 license-file = "LICENSE"
-rust-version = "1.82"
+rust-version = "1.94.0"
 description = "Orbit determination toolkit in Rust. Provides astrometric parsing, observer management, and initial orbit determination (Gauss method) with JPL ephemeris support."
 readme = "README.md"
 repository = "https://github.com/FusRoman/Outfit"
@@ -20,6 +20,9 @@ categories = ["science", "parsing", "data-structures"]
 authors = ["Roman Le Montagner <roman.le-montagner@ijclab.in2p3.fr>"]
 
 [dependencies]
+
+photom = { version = "0.1.0", path = "../photom" }
+
 aberth = { version = "0.4.1", default-features = false }
 ahash = { version = "0.8.11", default-features = false }
 arrow-array = { version = "54.3.1", default-features = false }
@@ -54,18 +57,18 @@ rand = { version = "0.9.2", default-features = false, features = [
 ] }
 rand_distr = { version = "0.5.1", default-features = false }
 
-reqwest = { version = "0.12.15", default-features = false, optional = true, features = [
+reqwest = { version = "0.12.15", default-features = false, features = [
     "http2",
     "rustls-tls",
     "stream",
 ] }
-tokio = { version = "1.44.1", default-features = false, optional = true, features = [
+tokio = { version = "1.44.1", default-features = false, features = [
     "fs",
     "rt",
     "rt-multi-thread",
     "io-util",
 ] }
-tokio-stream = { version = "0.1.17", default-features = false, optional = true }
+tokio-stream = { version = "0.1.17", default-features = false }
 indicatif = { version = "0.18", optional = true, default-features = false }
 rayon = { version = "1.11.0", optional = true, default-features = false }
 comfy-table = { version = "7.1.4", default-features = false }
@@ -76,21 +79,12 @@ criterion = { version = "0.5.1", features = ["html_reports"] }
 husky-rs = "0.1.5"
 proptest = "1.7.0"
 
+photom = { version = "0.1.0", path = "../photom", features = ["mpc_80_col", "ades"] }
+
 [features]
-jpl-download = ["dep:reqwest", "dep:tokio", "dep:tokio-stream"]
 progress = ["dep:indicatif"]
 parallel = ["dep:rayon"]
 
-[[test]]
-name = "reader_80col_test"
-path = "tests/reader_80col_test.rs"
-required-features = ["jpl-download"]
-
-[[test]]
-name = "test_read_ades"
-path = "tests/test_read_ades.rs"
-required-features = ["jpl-download"]
-
 [[test]]
 name = "test_large_parquet"
 path = "tests/trajectories_from_parquet.rs"
@@ -142,7 +136,3 @@ debug = false
 lto = "fat"
 codegen-units = 1
 strip = true
-
-[package.metadata.docs.rs]
-features = ["jpl-download"]
-rustdoc-args = ["--cfg", "docsrs"]

rust-toolchain.toml

@@ -1,3 +1,3 @@
 [toolchain]
-channel = "1.88.0"
+channel = "1.94.0"
 components = ["rustfmt", "clippy"]

src/cache/mod.rs

@@ -0,0 +1,77 @@
+pub mod observer_centric_cache;
+pub mod observer_fixed_cache;
+
+use hifitime::ut1::Ut1Provider;
+use photom::{observation_dataset::ObsDataset, observer::dataset::ObserverId, ObsIndex};
+
+use crate::{
+    cache::{
+        observer_centric_cache::{
+            build_centric_observer_cache, CentricObserverCache, ObserverCentricCache,
+            ObserverGeocentricPosition, ObserverGeocentricVelocity, ObserverHeliocentricPosition,
+        },
+        observer_fixed_cache::{
+            build_fixed_observer_cache, BodyFixedObserverCache, ObserverFixedCache,
+        },
+    },
+    JPLEphem, OutfitError,
+};
+
+/// Precomputed observer positions for all observations in a dataset.
+///
+/// Built once before any trajectory fitting. Each entry is keyed by the
+/// observation's [`ObsIndex`], which is stable for the lifetime of the
+/// [`ObsDataset`].
+pub struct OutfitCache {
+    /// len == number of observations in the dataset. Indexed by [`ObsIndex`].
+    observer_centric: CentricObserverCache,
+    /// len == number of observer in the dataset. Indexed by observer ID.
+    observer_fixed: BodyFixedObserverCache,
+}
+
+impl OutfitCache {
+    pub fn get_centric(&self, idx: ObsIndex) -> &ObserverCentricCache {
+        &self.observer_centric[idx]
+    }
+
+    pub fn get_fixed(&self, observer_id: ObserverId) -> Option<&ObserverFixedCache> {
+        self.observer_fixed.get(&observer_id)
+    }
+
+    /// Build the cache for every observation in `obs_dataset`.
+    pub fn build(
+        obs_dataset: &ObsDataset,
+        jpl: &JPLEphem,
+        ut1_provider: &Ut1Provider,
+    ) -> Result<Self, OutfitError> {
+        let observer_iter = obs_dataset.iter_observer()?;
+
+        let observer_fixed_cache = build_fixed_observer_cache(observer_iter)?;
+
+        let observer_centric_cache =
+            build_centric_observer_cache(jpl, ut1_provider, obs_dataset, &observer_fixed_cache)?;
+
+        Ok(Self {
+            observer_centric: observer_centric_cache,
+            observer_fixed: observer_fixed_cache,
+        })
+    }
+
+    pub fn get_observer_fixed_cache(&self, observer_id: ObserverId) -> Option<&ObserverFixedCache> {
+        self.observer_fixed.get(&observer_id)
+    }
+
+    /// Accessor for the precomputed geocentric position of an observer.
+    pub fn get_observer_geocentric_position(&self, idx: ObsIndex) -> &ObserverGeocentricPosition {
+        &self.get_centric(idx).geo_position
+    }
+
+    pub fn get_observer_geocentric_velocity(&self, idx: ObsIndex) -> &ObserverGeocentricVelocity {
+        &self.get_centric(idx).geo_velocity
+    }
+
+    /// Accessor for the precomputed heliocentric position of an observer.
+    pub fn get_helio_position(&self, idx: ObsIndex) -> &ObserverHeliocentricPosition {
+        &self.get_centric(idx).helio_position
+    }
+}

src/cache/observer_centric_cache.rs

@@ -0,0 +1,195 @@
+use hifitime::{ut1::Ut1Provider, Epoch};
+use nalgebra::Vector3;
+use ordered_float::NotNan;
+use photom::{observation_dataset::ObsDataset, observer::Observer, MJDTT};
+
+use crate::{
+    cache::observer_fixed_cache::{BodyFixedObserverCache, ObserverFixedCache},
+    observer_extension::ResolvedObserver,
+    JPLEphem, OutfitError,
+};
+
+/// Geocentric position of the observer at `time` of an observation (AU, ecliptic mean J2000).
+pub type ObserverGeocentricPosition = Vector3<NotNan<f64>>;
+/// Geocentric velocity of the observer at `time` of an observation (AU/day, ecliptic mean J2000).
+pub type ObserverGeocentricVelocity = Vector3<NotNan<f64>>;
+/// Heliocentric position of the observer at `time` of an observation (AU, ecliptic mean J2000).
+pub type ObserverHeliocentricPosition = Vector3<NotNan<f64>>;
+
+/// Geocentric and heliocentric observer positions for a single observation epoch.
+pub struct ObserverCentricCache {
+    pub geo_position: ObserverGeocentricPosition,
+    pub geo_velocity: ObserverGeocentricVelocity,
+    pub helio_position: ObserverHeliocentricPosition,
+}
+
+impl ObserverCentricCache {
+    pub fn new(
+        jpl: &JPLEphem,
+        ut1_provider: &Ut1Provider,
+        obs_time: MJDTT,
+        observer_fixed_cache: &ObserverFixedCache,
+    ) -> Result<Self, OutfitError> {
+        let obs_mjd = Epoch::from_mjd_in_time_scale(obs_time, hifitime::TimeScale::TT);
+        let (geocentric_pos, geocentric_vel) =
+            Observer::pvobs(&obs_mjd, ut1_provider, observer_fixed_cache)?;
+
+        let heliocentric_pos = Observer::helio_position(jpl, &obs_mjd, &geocentric_pos)?;
+
+        Ok(Self {
+            geo_position: geocentric_pos,
+            geo_velocity: geocentric_vel,
+            helio_position: heliocentric_pos,
+        })
+    }
+}
+
+pub type CentricObserverCache = Vec<ObserverCentricCache>;
+
+pub fn build_centric_observer_cache(
+    jpl: &JPLEphem,
+    ut1_provider: &Ut1Provider,
+    obs_dataset: &ObsDataset,
+    observer_fixed_cache: &BodyFixedObserverCache,
+) -> Result<CentricObserverCache, OutfitError> {
+    obs_dataset
+        .iter_observations()
+        .enumerate()
+        .map(|(idx, obs)| {
+            let observer_id = obs
+                .observer_id()
+                .ok_or_else(|| return OutfitError::ObserverIdIsNone(idx as u64))?;
+
+            let fixed_cache = observer_fixed_cache
+                .get(&observer_id)
+                .ok_or_else(|| return OutfitError::ObserverIdIsNone(idx as u64))?;
+
+            ObserverCentricCache::new(jpl, ut1_provider, obs.mjd_tt(), fixed_cache)
+        })
+        .collect()
+}
+
+#[cfg(test)]
+mod observer_test {
+
+    use photom::{Meters, Radians};
+
+    use crate::{
+        cache::observer_centric_cache::ObserverCentricCache,
+        test_fixture::{JPL_EPHEM_HORIZON, UT1_PROVIDER},
+    };
+
+    use super::*;
+
+    fn to_observer(
+        longitude: Radians,
+        latitude: Radians,
+        height: Meters,
+        name: Option<String>,
+        ra_accuracy: Option<f64>,
+        dec_accuracy: Option<f64>,
+    ) -> Observer {
+        let observer = Observer::new(longitude, latitude, height, name, ra_accuracy, dec_accuracy)
+            .expect("Failed to create Observer");
+        observer
+    }
+
+    #[test]
+    fn body_fixed_coord_test() {
+        // longitude, latitude and height of Pan-STARRS 1, Haleakala
+        let (lon, lat, h) = (
+            203.744090000_f64.to_radians(),
+            20.707233557_f64.to_radians(),
+            3067.694,
+        );
+        let pan_starrs = to_observer(lon, lat, h, None, None, None);
+        assert_eq!(
+            pan_starrs
+                .earth_fixed_position()
+                .unwrap()
+                .map(|x| x.into_inner()),
+            Vector3::new(
+                -0.00003653799439776371,
+                -0.00001607260397528885,
+                0.000014988110430544328
+            )
+        );
+    }
+
+    #[test]
+    fn pvobs_test() {
+        let tmjd = 57028.479297592596;
+        let epoch = Epoch::from_mjd_in_time_scale(tmjd, hifitime::TimeScale::TT);
+        // longitude, latitude and height of Pan-STARRS 1, Haleakala
+        let (lon, lat, h) = (
+            203.744090000_f64.to_radians(),
+            20.707233557_f64.to_radians(),
+            3067.694,
+        );
+
+        let pan_starrs = to_observer(lon, lat, h, Some("Pan-STARRS 1".to_string()), None, None);
+
+        let observer_fixed_cache: ObserverFixedCache = (&pan_starrs).try_into().unwrap();
+
+        let (observer_position, observer_velocity) =
+            Observer::pvobs(&epoch, &UT1_PROVIDER, &observer_fixed_cache).unwrap();
+
+        assert_eq!(
+            observer_position.as_slice(),
+            [
+                -2.086211182493635e-5,
+                3.718476815327979e-5,
+                2.4978996447997476e-7
+            ]
+        );
+        assert_eq!(
+            observer_velocity.as_slice(),
+            [
+                -0.0002143246535691577,
+                -0.00012059801691431748,
+                5.262184624215718e-5
+            ]
+        );
+    }
+
+    #[test]
+    fn test_helio_pos_obs() {
+        let (lon, lat, h) = (203.744090000_f64, 20.707233557_f64, 3067.694_f64);
+        let pan_starrs = to_observer(
+            lon.to_radians(),
+            lat.to_radians(),
+            h,
+            Some("Pan-STARRS 1".to_string()),
+            None,
+            None,
+        );
+        let observer_fixed_cache: ObserverFixedCache = (&pan_starrs).try_into().unwrap();
+
+        let cases = [
+            (
+                57_028.479_297_592_596,
+                [-0.2645666171464416, 0.8689351643701766, 0.3766996211107864],
+            ),
+            (
+                57_049.245_147_592_59,
+                [-0.5891631852137064, 0.7238872516824697, 0.3138186516540669],
+            ),
+            (
+                57_063.977_117_592_59,
+                [-0.7743280306286537, 0.5612532665812755, 0.24333415479994636],
+            ),
+        ];
+
+        for (tmjd, expected) in cases {
+            let obs = ObserverCentricCache::new(
+                &JPL_EPHEM_HORIZON,
+                &UT1_PROVIDER,
+                tmjd,
+                &observer_fixed_cache,
+            )
+            .unwrap();
+
+            assert_eq!(obs.helio_position.as_slice(), expected, "tmjd = {tmjd}");
+        }
+    }
+}