From 29a10713208961c54d5e765b5796bd2085c46f00 Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Thu, 9 Oct 2025 21:27:47 +0200 Subject: [PATCH 1/5] Consistent `Time` usage throughout rust backend --- CHANGELOG.md | 19 +++++ src/kete/rust/covariance.rs | 9 +- src/kete/rust/elements.rs | 22 ++--- src/kete/rust/fovs/checks.rs | 19 +++-- src/kete/rust/fovs/collection.rs | 4 +- src/kete/rust/fovs/definitions.rs | 76 +++++++++-------- src/kete/rust/frame.rs | 2 +- src/kete/rust/horizons.rs | 6 +- src/kete/rust/kepler.rs | 27 +++--- src/kete/rust/propagation.rs | 17 ++-- src/kete/rust/simult_states.rs | 5 +- src/kete/rust/spice/pck.rs | 4 +- src/kete/rust/spice/spk.rs | 10 +-- src/kete/rust/state.rs | 10 +-- src/kete/rust/state_transition.rs | 4 +- src/kete_core/benches/propagation.rs | 18 ++-- src/kete_core/benches/spice.rs | 8 +- src/kete_core/src/elements.rs | 53 +++++++----- src/kete_core/src/errors.rs | 12 ++- src/kete_core/src/fov/fov_like.rs | 16 ++-- src/kete_core/src/fov/generic.rs | 33 ++++--- src/kete_core/src/fov/neos.rs | 2 +- src/kete_core/src/fov/ptf.rs | 2 +- src/kete_core/src/fov/spherex.rs | 2 +- src/kete_core/src/fov/ztf.rs | 2 +- src/kete_core/src/frames/earth.rs | 12 +-- src/kete_core/src/io/parquet.rs | 4 +- src/kete_core/src/propagation/acceleration.rs | 26 +++--- src/kete_core/src/propagation/kepler.rs | 35 ++++---- src/kete_core/src/propagation/mod.rs | 39 +++++---- src/kete_core/src/propagation/nongrav.rs | 2 +- src/kete_core/src/propagation/picard.rs | 58 +++++++------ src/kete_core/src/propagation/radau.rs | 32 +++---- src/kete_core/src/propagation/runge_kutta.rs | 56 ++++++------ .../src/propagation/state_transition.rs | 9 +- src/kete_core/src/propagation/util.rs | 9 +- src/kete_core/src/simult_states.rs | 11 +-- src/kete_core/src/spice/daf.rs | 5 +- src/kete_core/src/spice/mod.rs | 18 ++-- src/kete_core/src/spice/pck.rs | 3 +- src/kete_core/src/spice/pck_segments.rs | 7 +- src/kete_core/src/spice/sclk.rs | 5 +- src/kete_core/src/spice/spk.rs | 31 +++---- src/kete_core/src/spice/spk_segments.rs | 4 +- src/kete_core/src/state.rs | 53 +++++++----- src/kete_core/src/time/mod.rs | 85 +++++++++++++++++-- 46 files changed, 525 insertions(+), 361 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fe4f93b..75bf43b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,25 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## [Unreleased] + +### Added + +- Added a new numerical integrator based on Picard-Chebyshev integration. This + integrator has only been added to the rust backend at this point, until more + testing can be done and it be made available on the frontend. + +### Changed + +- Throughout the rust code, `Time` is being enforced as inputs for functions instead + of accepting `f64` in a large number of places. + +### Fixed + +- Fixed epoch times in PCK Type 2 frames that were not being converted correctly. + + ## [v2.1.5] ### Added diff --git a/src/kete/rust/covariance.rs b/src/kete/rust/covariance.rs index e26cdff..c99d7ba 100644 --- a/src/kete/rust/covariance.rs +++ b/src/kete/rust/covariance.rs @@ -3,6 +3,7 @@ use std::{collections::HashMap, fmt::Debug}; use crate::state::PyState; +use crate::time::PyTime; use crate::{elements::PyCometElements, vector::VectorLike}; use kete_core::{errors::Error, io::FileIO}; use pyo3::prelude::*; @@ -34,13 +35,13 @@ impl Covariance { #[allow(clippy::too_many_arguments)] pub fn new( desig: String, - epoch: f64, + epoch: PyTime, params: Vec<(String, f64)>, cov_matrix: Vec>, ) -> Self { Self { desig, - epoch, + epoch: epoch.jd(), params, cov_matrix, } @@ -98,12 +99,12 @@ impl Covariance { .ok_or(Error::ValueError("Covariance missing 'inclination'".into()))?; let elem = PyCometElements::new( desig, - epoch, + epoch.into(), eccentricity, inclination, peri_dist, peri_arg, - peri_time, + peri_time.into(), lon_of_ascending, ); elem.state() diff --git a/src/kete/rust/elements.rs b/src/kete/rust/elements.rs index 3d04181..66ac4c7 100644 --- a/src/kete/rust/elements.rs +++ b/src/kete/rust/elements.rs @@ -4,7 +4,7 @@ use kete_core::frames::Ecliptic; use kete_core::prelude; use pyo3::{PyResult, pyclass, pymethods}; -use crate::state::PyState; +use crate::{state::PyState, time::PyTime}; /// Cometary Elements class made accessible to python. /// @@ -60,21 +60,21 @@ impl PyCometElements { #[allow(clippy::too_many_arguments)] pub fn new( desig: String, - epoch: f64, + epoch: PyTime, eccentricity: f64, inclination: f64, peri_dist: f64, peri_arg: f64, - peri_time: f64, + peri_time: PyTime, lon_of_ascending: f64, ) -> Self { Self(elements::CometElements { desig: prelude::Desig::Name(desig), - epoch, + epoch: epoch.into(), eccentricity, inclination: inclination.to_radians(), lon_of_ascending: lon_of_ascending.to_radians(), - peri_time, + peri_time: peri_time.into(), peri_arg: peri_arg.to_radians(), peri_dist, }) @@ -93,8 +93,8 @@ impl PyCometElements { /// Epoch of the elements in JD. #[getter] - pub fn epoch(&self) -> f64 { - self.0.epoch + pub fn epoch(&self) -> PyTime { + self.0.epoch.into() } /// Designation of the object. @@ -128,8 +128,8 @@ impl PyCometElements { /// Perihelion time of the orbit in JD. #[getter] - pub fn peri_time(&self) -> f64 { - self.0.peri_time + pub fn peri_time(&self) -> PyTime { + self.0.peri_time.into() } /// Argument of Perihelion of the orbit in degrees. @@ -196,11 +196,11 @@ impl PyCometElements { format!( "CometElements(desig={:?}, epoch={}, eccentricity={}, inclination={}, lon_of_ascending={}, peri_time={}, peri_arg={}, peri_dist={})", self.desig(), - self.epoch(), + self.epoch().jd(), self.eccentricity(), self.inclination(), self.lon_of_ascending(), - self.peri_time(), + self.peri_time().jd(), self.peri_arg(), self.peri_dist() ) diff --git a/src/kete/rust/fovs/checks.rs b/src/kete/rust/fovs/checks.rs index 8eea1da..1d88c06 100644 --- a/src/kete/rust/fovs/checks.rs +++ b/src/kete/rust/fovs/checks.rs @@ -42,10 +42,10 @@ pub fn fov_checks_py( chunk.push(fov); continue; }; - let jd_start = chunk.first().unwrap().observer().jd; + let jd_start = chunk.first().unwrap().observer().epoch; // chunk is complete - if (fov.observer().jd - jd_start).abs() >= 2.0 * dt_limit { + if (fov.observer().epoch - jd_start).elapsed.abs() >= 2.0 * dt_limit { fov_chunks.push(chunk); chunk = vec![fov]; } else { @@ -55,21 +55,24 @@ pub fn fov_checks_py( if !chunk.is_empty() { fov_chunks.push(chunk); } - let mut jd = pop.jd; + let mut jd = pop.epoch.jd; let mut big_jd = jd; let mut states = pop.states; let mut big_step_states = states.clone(); let mut visible = Vec::new(); for fovs in fov_chunks.into_iter() { - let jd_mean = - (fovs.last().unwrap().observer().jd + fovs.first().unwrap().observer().jd) / 2.0; + let jd_mean = (fovs.last().unwrap().observer().epoch.jd + + fovs.first().unwrap().observer().epoch.jd) + / 2.0; // Take large steps which are 10x the smaller steps, this helps long term numerical stability if (jd_mean - big_jd).abs() >= dt_limit * 50.0 { big_jd = jd_mean; big_step_states = big_step_states .into_par_iter() - .filter_map(|state| propagate_n_body_spk(state, jd, include_asteroids, None).ok()) + .filter_map(|state| { + propagate_n_body_spk(state, jd.into(), include_asteroids, None).ok() + }) .collect(); }; // Take small steps based off of the large steps. @@ -80,7 +83,9 @@ pub fn fov_checks_py( jd = jd_mean; states = states .into_par_iter() - .filter_map(|state| propagate_n_body_spk(state, jd, include_asteroids, None).ok()) + .filter_map(|state| { + propagate_n_body_spk(state, jd.into(), include_asteroids, None).ok() + }) .collect(); }; diff --git a/src/kete/rust/fovs/collection.rs b/src/kete/rust/fovs/collection.rs index c098fa8..840238e 100644 --- a/src/kete/rust/fovs/collection.rs +++ b/src/kete/rust/fovs/collection.rs @@ -25,7 +25,7 @@ impl FOVListLike { FOVListLike::Vec(v) => v, FOVListLike::FOVList(list) => list.0, }; - fovs.sort_by(|a, b| a.jd().total_cmp(&b.jd())); + fovs.sort_by(|a, b| a.jd().jd.total_cmp(&b.jd().jd)); fovs.into_iter().map(|fov| fov.unwrap()).collect() } } @@ -45,7 +45,7 @@ impl FOVList { /// Sort the list of FOVs by their JD. pub fn sort(&mut self) { - self.0.sort_by(|a, b| a.jd().total_cmp(&b.jd())) + self.0.sort_by(|a, b| a.jd().jd.total_cmp(&b.jd().jd)) } fn __len__(&self) -> usize { diff --git a/src/kete/rust/fovs/definitions.rs b/src/kete/rust/fovs/definitions.rs index bd080a3..20b91b9 100644 --- a/src/kete/rust/fovs/definitions.rs +++ b/src/kete/rust/fovs/definitions.rs @@ -1,8 +1,10 @@ use kete_core::fov::{self}; use kete_core::fov::{FovLike, SkyPatch}; use kete_core::frames::Vector; +use kete_core::time::{TDB, Time}; use pyo3::{exceptions, prelude::*}; +use crate::time::PyTime; use crate::vector::VectorLike; use crate::{state::PyState, vector::PyVector}; @@ -256,20 +258,20 @@ pub enum AllowedFOV { impl AllowedFOV { #[allow(missing_docs)] - pub fn jd(&self) -> f64 { + pub fn jd(&self) -> Time { match self { - AllowedFOV::NEOS(fov) => fov.0.observer().jd, - AllowedFOV::WISE(fov) => fov.0.observer().jd, - AllowedFOV::Rectangle(fov) => fov.0.observer().jd, - AllowedFOV::ZTF(fov) => fov.0.observer().jd, - AllowedFOV::ZTFField(fov) => fov.0.observer().jd, - AllowedFOV::NEOSVisit(fov) => fov.0.observer().jd, - AllowedFOV::Cone(fov) => fov.0.observer().jd, - AllowedFOV::OmniDirectional(fov) => fov.0.observer().jd, - AllowedFOV::PTF(fov) => fov.0.observer().jd, - AllowedFOV::PTFField(fov) => fov.0.observer().jd, - AllowedFOV::SPHEREx(fov) => fov.0.observer().jd, - AllowedFOV::SPHERExField(fov) => fov.0.observer().jd, + AllowedFOV::NEOS(fov) => fov.0.observer().epoch, + AllowedFOV::WISE(fov) => fov.0.observer().epoch, + AllowedFOV::Rectangle(fov) => fov.0.observer().epoch, + AllowedFOV::ZTF(fov) => fov.0.observer().epoch, + AllowedFOV::ZTFField(fov) => fov.0.observer().epoch, + AllowedFOV::NEOSVisit(fov) => fov.0.observer().epoch, + AllowedFOV::Cone(fov) => fov.0.observer().epoch, + AllowedFOV::OmniDirectional(fov) => fov.0.observer().epoch, + AllowedFOV::PTF(fov) => fov.0.observer().epoch, + AllowedFOV::PTFField(fov) => fov.0.observer().epoch, + AllowedFOV::SPHEREx(fov) => fov.0.observer().epoch, + AllowedFOV::SPHERExField(fov) => fov.0.observer().epoch, } } @@ -396,8 +398,8 @@ impl PyWiseCmos { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Direction that the observer is looking. @@ -494,8 +496,8 @@ impl PyGenericRectangle { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Direction that the observer is looking. @@ -560,8 +562,8 @@ impl PyGenericCone { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Direction that the observer is looking. @@ -602,8 +604,8 @@ impl PyOmniDirectional { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } fn __repr__(&self) -> String { @@ -653,8 +655,8 @@ impl PyNeosCmos { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Direction that the observer is looking. @@ -798,8 +800,8 @@ impl PyNeosVisit { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Metadata about where this FOV is in the Survey. @@ -946,8 +948,8 @@ impl PyZtfCcdQuad { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Metadata about where this FOV is in the Survey. @@ -1053,8 +1055,8 @@ impl PyZtfField { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Metadata about where this FOV is in the Survey. @@ -1184,8 +1186,8 @@ impl PyPtfCcd { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Metadata about where this FOV is in the Survey. @@ -1279,8 +1281,8 @@ impl PyPtfField { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Metadata about where this FOV is in the Survey. @@ -1376,8 +1378,8 @@ impl PySpherexCmos { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Metadata about where this FOV is in the Survey. @@ -1454,8 +1456,8 @@ impl PySpherexField { /// JD of the observer location. #[getter] - pub fn jd(&self) -> f64 { - self.0.observer().jd + pub fn jd(&self) -> PyTime { + self.0.observer().epoch.into() } /// Metadata about where this FOV is in the Survey. diff --git a/src/kete/rust/frame.rs b/src/kete/rust/frame.rs index 828d7e6..c0c6f9e 100644 --- a/src/kete/rust/frame.rs +++ b/src/kete/rust/frame.rs @@ -152,7 +152,7 @@ pub fn calc_obliquity_py(time: f64) -> f64 { /// Time in TDB scaled Julian Days. #[pyfunction] #[pyo3(name = "earth_precession_rotation")] -pub fn calc_earth_precession(time: f64) -> Vec> { +pub fn calc_earth_precession(time: PyTime) -> Vec> { earth_precession_rotation(time.into()) .rotations_to_equatorial() .unwrap() diff --git a/src/kete/rust/horizons.rs b/src/kete/rust/horizons.rs index badb26f..938cc1b 100644 --- a/src/kete/rust/horizons.rs +++ b/src/kete/rust/horizons.rs @@ -118,7 +118,8 @@ impl HorizonsProperties { desig: prelude::Desig::Name(self.desig.clone()), epoch: self .epoch - .ok_or(prelude::Error::ValueError("No Epoch defined".into()))?, + .ok_or(prelude::Error::ValueError("No Epoch defined".into()))? + .into(), eccentricity: self .eccentricity .ok_or(prelude::Error::ValueError("No Eccentricity defined".into()))?, @@ -135,7 +136,8 @@ impl HorizonsProperties { .ok_or(prelude::Error::ValueError("No peri_dist defined".into()))?, peri_time: self .peri_time - .ok_or(prelude::Error::ValueError("No peri_time defined".into()))?, + .ok_or(prelude::Error::ValueError("No peri_time defined".into()))? + .into(), lon_of_ascending: self .lon_of_ascending .ok_or(prelude::Error::ValueError( diff --git a/src/kete/rust/kepler.rs b/src/kete/rust/kepler.rs index 0b8d6e4..e49c8e9 100644 --- a/src/kete/rust/kepler.rs +++ b/src/kete/rust/kepler.rs @@ -55,7 +55,7 @@ pub fn compute_eccentric_anomaly_py( /// ---------- /// state : /// List of states, which are in units of AU from the Sun and velocity is in AU/Day. -/// jd : +/// epoch : /// Time to integrate to in JD days with TDB scaling. /// observer_pos : /// A vector of length 3 describing the position of an observer. If this is @@ -67,15 +67,15 @@ pub fn compute_eccentric_anomaly_py( /// State /// Final states after propagating to the target time. #[pyfunction] -#[pyo3(name = "propagate_two_body", signature = (states, jd, observer_pos=None))] +#[pyo3(name = "propagate_two_body", signature = (states, epoch, observer_pos=None))] pub fn propagation_kepler_py( py: Python<'_>, states: MaybeVec, - jd: PyTime, + epoch: PyTime, observer_pos: Option, ) -> PyResult { let (states, was_vec): (Vec<_>, bool) = states.into(); - let jd = jd.jd(); + let epoch = epoch.into(); let states = states .par_iter() .with_min_len(10) @@ -85,24 +85,25 @@ pub fn propagation_kepler_py( let Some(state) = state.change_center(crate::desigs::NaifIDLike::Int(10)).ok() else { let nan_state: PyState = - State::::new_nan(state.raw.desig.clone(), jd, center).into(); + State::::new_nan(state.raw.desig.clone(), epoch, center).into(); return nan_state.change_frame(frame); }; - let Some(mut new_state) = propagation::propagate_two_body(&state.raw, jd).ok() else { + let Some(mut new_state) = propagation::propagate_two_body(&state.raw, epoch).ok() + else { let nan_state: PyState = - State::::new_nan(state.raw.desig.clone(), jd, center).into(); + State::::new_nan(state.raw.desig.clone(), epoch, center).into(); return nan_state.change_frame(frame); }; if let Some(observer_pos) = &observer_pos { let observer_pos: Vector = observer_pos.clone().into(); let delay = -(new_state.pos - observer_pos).norm() / constants::C_AU_PER_DAY; - new_state = match propagation::propagate_two_body(&new_state, new_state.jd + delay) - { - Ok(state) => state, - Err(_) => State::new_nan(state.raw.desig.clone(), jd, center), - }; + new_state = + match propagation::propagate_two_body(&new_state, new_state.epoch + delay) { + Ok(state) => state, + Err(_) => State::new_nan(state.raw.desig.clone(), epoch, center), + }; } let new_pystate: PyState = new_state.into(); @@ -110,7 +111,7 @@ pub fn propagation_kepler_py( .change_frame(frame) .change_center(crate::desigs::NaifIDLike::Int(center)) .unwrap_or( - State::::new_nan(state.raw.desig, jd, state.raw.center_id).into(), + State::::new_nan(state.raw.desig, epoch, state.raw.center_id).into(), ) }) .collect(); diff --git a/src/kete/rust/propagation.rs b/src/kete/rust/propagation.rs index b9a26e5..3c4cc17 100644 --- a/src/kete/rust/propagation.rs +++ b/src/kete/rust/propagation.rs @@ -6,7 +6,6 @@ use kete_core::{ propagation::{self, NonGravModel, moid}, spice::{self, LOADED_SPK}, state::State, - time::{TDB, Time}, }; use pyo3::{PyObject, PyResult, Python, pyfunction}; use rayon::prelude::*; @@ -48,7 +47,7 @@ pub fn moid_py( .map(|x| x.raw) .unwrap_or(LOADED_SPK.read().unwrap().try_get_state_with_center( 399, - states[0].raw.jd, + states[0].raw.epoch, 10, )?); @@ -116,7 +115,7 @@ pub fn propagation_n_body_spk_py( } let mut res: Vec = Vec::new(); - let jd = jd.jd(); + let jd = jd.into(); // propagation is broken into chunks, every time a chunk is completed // python is checked for signals. This allows keyboard interrupts to be caught @@ -164,8 +163,8 @@ pub fn propagation_n_body_spk_py( "Impact detected between ({}) <-> {} at time {} ({})", desig, spice::try_name_from_id(id).unwrap_or(id.to_string()), - time, - Time::::new(time).utc().to_iso().unwrap() + time.jd, + time.utc().to_iso().unwrap() ); } // if we get an impact, we return a state with NaNs @@ -256,7 +255,7 @@ pub fn propagation_n_body_py( let non_gravs: Vec> = non_gravs.into_iter().map(|y| y.map(|z| z.0)).collect(); - let jd = jd_final.jd(); + let jd = jd_final.into(); let res = states .into_iter() .zip(non_gravs.into_iter()) @@ -309,7 +308,7 @@ pub fn picard( } let mut res: Vec = Vec::new(); - let jd = jd.jd(); + let jd = jd.into(); // propagation is broken into chunks, every time a chunk is completed // python is checked for signals. This allows keyboard interrupts to be caught @@ -357,8 +356,8 @@ pub fn picard( "Impact detected between ({}) <-> {} at time {} ({})", desig, spice::try_name_from_id(id).unwrap_or(id.to_string()), - time, - Time::::new(time).utc().to_iso().unwrap() + time.jd, + time.utc().to_iso().unwrap() ); } // if we get an impact, we return a state with NaNs diff --git a/src/kete/rust/simult_states.rs b/src/kete/rust/simult_states.rs index 9db4285..ea811db 100644 --- a/src/kete/rust/simult_states.rs +++ b/src/kete/rust/simult_states.rs @@ -7,6 +7,7 @@ use pyo3::exceptions; use pyo3::prelude::*; use pyo3::{PyResult, pyclass, pymethods}; +use crate::time::PyTime; use crate::vector::PyVector; use crate::{fovs::AllowedFOV, state::PyState}; @@ -90,8 +91,8 @@ impl PySimultaneousStates { /// The time of the simultaneous states. #[getter] - pub fn jd(&self) -> f64 { - self.0.jd + pub fn jd(&self) -> PyTime { + self.0.epoch.into() } /// The reference center NAIF ID for this state. diff --git a/src/kete/rust/spice/pck.rs b/src/kete/rust/spice/pck.rs index 7f5bdf6..5d796db 100644 --- a/src/kete/rust/spice/pck.rs +++ b/src/kete/rust/spice/pck.rs @@ -46,7 +46,7 @@ pub fn pck_earth_frame_py( new_center: i32, name: Option, ) -> PyResult { - let jd = jd.jd(); + let jd = jd.into(); let desig = { match name { Some(d) => Desig::Name(d), @@ -89,7 +89,7 @@ pub fn pck_state_to_earth(state: PyState) -> PyResult<(f64, f64, f64)> { let state = state .change_center(crate::desigs::NaifIDLike::Int(399))? .raw; - let frame = pcks.try_get_orientation(3000, state.jd)?; + let frame = pcks.try_get_orientation(3000, state.epoch)?; let (pos, _) = frame.from_equatorial(state.pos.into(), state.vel.into())?; let [x, y, z] = pos.into(); diff --git a/src/kete/rust/spice/spk.rs b/src/kete/rust/spice/spk.rs index ea78686..3db2e81 100644 --- a/src/kete/rust/spice/spk.rs +++ b/src/kete/rust/spice/spk.rs @@ -33,7 +33,7 @@ pub fn spk_load_py(py: Python<'_>, filenames: Vec) -> PyResult<()> { /// (name, JD_start, JD_end, Center Naif ID, Frame ID, SPK Segment type ID) #[pyfunction] #[pyo3(name = "_loaded_object_info")] -pub fn spk_available_info_py(naif_id: NaifIDLike) -> Vec<(String, f64, f64, i32, i32, i32)> { +pub fn spk_available_info_py(naif_id: NaifIDLike) -> Vec<(String, PyTime, PyTime, i32, i32, i32)> { let (name, naif_id) = naif_id.try_into().unwrap(); let singleton = &LOADED_SPK.try_read().unwrap(); singleton @@ -42,8 +42,8 @@ pub fn spk_available_info_py(naif_id: NaifIDLike) -> Vec<(String, f64, f64, i32, .map(|(jd_start, jd_end, center_id, frame_id, segment_id)| { ( name.clone(), - jd_start, - jd_end, + jd_start.into(), + jd_end.into(), center_id, frame_id, segment_id, @@ -130,7 +130,7 @@ pub fn spk_state_py( center: NaifIDLike, frame: PyFrames, ) -> PyResult { - let jd = jd.jd(); + let jd = jd.into(); let (_, center) = center.try_into()?; match id.clone().try_into() { Ok((_, id)) => { @@ -171,7 +171,7 @@ pub fn spk_state_py( #[pyo3(name = "spk_raw_state")] pub fn spk_raw_state_py(id: NaifIDLike, jd: PyTime) -> PyResult { let (_, id) = id.try_into()?; - let jd = jd.jd(); + let jd = jd.into(); let spk = &LOADED_SPK.try_read().unwrap(); Ok(PyState { raw: spk.try_get_state(id, jd)?, diff --git a/src/kete/rust/state.rs b/src/kete/rust/state.rs index bc7bf08..37dbcaf 100644 --- a/src/kete/rust/state.rs +++ b/src/kete/rust/state.rs @@ -85,7 +85,7 @@ impl PyState { let vel = vel.into_vector(frame); let center_id = center_id.unwrap_or(10); - let state = State::new(desig, jd.jd(), pos, vel, center_id); + let state = State::new(desig, jd.into(), pos, vel, center_id); Self { raw: state, frame, @@ -144,8 +144,8 @@ impl PyState { /// JD of the object's state in TDB scaled time. #[getter] - pub fn jd(&self) -> f64 { - self.raw.jd + pub fn jd(&self) -> PyTime { + self.raw.epoch.into() } /// Position of the object in AU with respect to the central object. @@ -171,7 +171,7 @@ impl PyState { pub fn is_finite(&self) -> bool { self.raw.pos.norm().is_finite() && self.raw.vel.norm().is_finite() - && self.raw.jd.is_finite() + && self.raw.epoch.jd.is_finite() } /// Central ID of the object used as reference for the coordinate frame. @@ -209,7 +209,7 @@ impl PyState { /// Perihelion time of the orbit in JD. #[getter] - pub fn peri_time(&mut self) -> f64 { + pub fn peri_time(&mut self) -> PyTime { self.elements().peri_time() } diff --git a/src/kete/rust/state_transition.rs b/src/kete/rust/state_transition.rs index 78fd51a..24beefc 100644 --- a/src/kete/rust/state_transition.rs +++ b/src/kete/rust/state_transition.rs @@ -16,13 +16,13 @@ pub fn compute_stm_py( ) -> ([[f64; 3]; 2], [[f64; 6]; 6]) { let mut state = State::new( Desig::Empty, - jd_start.jd(), + jd_start.into(), [state[0], state[1], state[2]].into(), [state[3], state[4], state[5]].into(), 10, ); - let (final_state, stm) = compute_state_transition(&mut state, jd_end.jd(), central_mass); + let (final_state, stm) = compute_state_transition(&mut state, jd_end.into(), central_mass); (final_state, stm.into()) } diff --git a/src/kete_core/benches/propagation.rs b/src/kete_core/benches/propagation.rs index 24d4aee..159aabc 100644 --- a/src/kete_core/benches/propagation.rs +++ b/src/kete_core/benches/propagation.rs @@ -14,7 +14,7 @@ use rayon::iter::{IntoParallelIterator, ParallelIterator}; static CIRCULAR: std::sync::LazyLock> = std::sync::LazyLock::new(|| { State::new( Desig::Name("Circular".into()), - 2451545.0, + 2451545.0.into(), [0.0, 1., 0.0].into(), [-constants::GMS_SQRT, 0.0, 0.0].into(), 0, @@ -23,7 +23,7 @@ static CIRCULAR: std::sync::LazyLock> = std::sync::LazyLock::new static ELLIPTICAL: std::sync::LazyLock> = std::sync::LazyLock::new(|| { State::new( Desig::Name("Elliptical".into()), - 2451545.0, + 2451545.0.into(), [0.0, 1.5, 0.0].into(), [-constants::GMS_SQRT, 0.0, 0.0].into(), 0, @@ -32,7 +32,7 @@ static ELLIPTICAL: std::sync::LazyLock> = std::sync::LazyLock::n static PARABOLIC: std::sync::LazyLock> = std::sync::LazyLock::new(|| { State::new( Desig::Name("Parabolic".into()), - 2451545.0, + 2451545.0.into(), [0.0, 2., 0.0].into(), [-constants::GMS_SQRT, 0.0, 0.0].into(), 0, @@ -42,7 +42,7 @@ static PARABOLIC: std::sync::LazyLock> = std::sync::LazyLock::ne static HYPERBOLIC: std::sync::LazyLock> = std::sync::LazyLock::new(|| { State::new( Desig::Name("Hyperbolic".into()), - 2451545.0, + 2451545.0.into(), [0.0, 3., 0.0].into(), [-constants::GMS_SQRT, 0.0, 0.0].into(), 0, @@ -50,14 +50,14 @@ static HYPERBOLIC: std::sync::LazyLock> = std::sync::LazyLock::n }); fn prop_n_body_radau(state: State, dt: f64) { - let jd = state.jd + dt; + let jd = state.epoch + dt; let _ = propagate_n_body_spk(state.into_frame(), jd, false, None).unwrap(); } fn prop_n_body_vec_radau(mut state: State, dt: f64) { let spk = &LOADED_SPK.read().unwrap(); spk.try_change_center(&mut state, 10).unwrap(); - let jd = state.jd + dt; + let jd = state.epoch + dt; let states = vec![state.into_frame().clone(); 100]; let non_gravs = vec![None; 100]; let _ = propagation::propagate_n_body_vec(states, jd, None, non_gravs).unwrap(); @@ -68,19 +68,19 @@ fn prop_n_body_radau_par(state: State, dt: f64) { let _tmp: Vec> = states .into_par_iter() .map(|s| { - let jd = s.jd + dt; + let jd = s.epoch + dt; propagate_n_body_spk(s.into_frame(), jd, false, None).unwrap() }) .collect(); } fn prop_2_body_radau(state: State, dt: f64) { - let jd = state.jd + dt; + let jd = state.epoch + dt; let _ = propagation::propagation_central(&state.into_frame(), jd).unwrap(); } fn prop_2_body_kepler(state: State, dt: f64) { - let _ = propagate_two_body(&state, state.jd + dt).unwrap(); + let _ = propagate_two_body(&state, state.epoch + dt).unwrap(); } /// Benchmark functions for the propagation algorithms diff --git a/src/kete_core/benches/spice.rs b/src/kete_core/benches/spice.rs index a51189c..5f0b261 100644 --- a/src/kete_core/benches/spice.rs +++ b/src/kete_core/benches/spice.rs @@ -14,7 +14,7 @@ use std::hint::black_box; fn spice_get_raw_state(jd: f64) { let spice = &LOADED_SPK.try_read().unwrap(); for _ in 0..1000 { - let _: State = spice.try_get_state(5, jd).unwrap(); + let _: State = spice.try_get_state(5, jd.into()).unwrap(); } } @@ -29,13 +29,15 @@ fn spice_change_center(mut state: State) { fn spice_get_state(jd: f64) { let spice = &LOADED_SPK.try_read().unwrap(); for _ in 0..1000 { - let _: State = spice.try_get_state_with_center(5, jd, 10).unwrap(); + let _: State = spice.try_get_state_with_center(5, jd.into(), 10).unwrap(); } } pub fn spice_benchmark(c: &mut Criterion) { let spice = &LOADED_SPK.try_read().unwrap(); - let state = spice.try_get_state_with_center(5, 2451545.0, 10).unwrap(); + let state = spice + .try_get_state_with_center(5, 2451545.0.into(), 10) + .unwrap(); c.bench_function("spice_get_raw_state", |b| { b.iter(|| spice_get_raw_state(black_box(2451545.0))); }); diff --git a/src/kete_core/src/elements.rs b/src/kete_core/src/elements.rs index 3ae1c48..abdb1da 100644 --- a/src/kete_core/src/elements.rs +++ b/src/kete_core/src/elements.rs @@ -35,6 +35,7 @@ use crate::constants::GMS_SQRT; use crate::frames::Ecliptic; use crate::prelude::{Desig, KeteResult, State}; use crate::propagation::{PARABOLIC_ECC_LIMIT, compute_eccentric_anomaly, compute_true_anomaly}; +use crate::time::{TDB, Time}; use nalgebra::Vector3; use std::f64::consts::TAU; @@ -53,7 +54,7 @@ pub struct CometElements { pub desig: Desig, /// Epoch of fit - pub epoch: f64, + pub epoch: Time, /// Eccentricity pub eccentricity: f64, @@ -65,7 +66,7 @@ pub struct CometElements { pub lon_of_ascending: f64, /// Time of perihelion passage in JD TDB scaled time - pub peri_time: f64, + pub peri_time: Time, /// Argument of perihelion in radians pub peri_arg: f64, @@ -79,7 +80,7 @@ impl CometElements { pub fn from_state(state: &State) -> Self { Self::from_pos_vel( state.desig.clone(), - state.jd, + state.epoch, &state.pos.into(), &state.vel.into(), ) @@ -90,7 +91,13 @@ impl CometElements { /// The units of the vectors are AU and AU/Day, with the central mass assumed /// to be the Sun. /// - fn from_pos_vel(desig: Desig, epoch: f64, pos: &Vector3, vel: &Vector3) -> Self { + fn from_pos_vel( + desig: Desig, + epoch: Time, + pos: &Vector3, + vel: &Vector3, + ) -> Self { + let epoch = epoch.jd; let vel_scaled = vel / GMS_SQRT; let v_mag2 = vel_scaled.norm_squared(); let p_mag = pos.norm(); @@ -181,11 +188,11 @@ impl CometElements { Self { desig, - epoch, + epoch: epoch.into(), eccentricity: ecc, inclination: incl, lon_of_ascending: lon_of_asc, - peri_time, + peri_time: peri_time.into(), peri_arg, peri_dist, } @@ -226,7 +233,7 @@ impl CometElements { semi_major.abs().powf(-1.5) * GMS_SQRT }; - let mean_anom = mean_motion * (self.epoch - self.peri_time); + let mean_anom = mean_motion * (self.epoch - self.peri_time).elapsed; let ecc_anom = compute_eccentric_anomaly(self.eccentricity, mean_anom, self.peri_dist)?; let x: f64; @@ -335,7 +342,7 @@ impl CometElements { /// Compute the mean anomaly in radians. pub fn mean_anomaly(&self) -> f64 { let mm = self.mean_motion(); - let mean_anomaly = (self.epoch - self.peri_time) * mm; + let mean_anomaly = (self.epoch - self.peri_time).elapsed * mm; match self.eccentricity { ecc if ecc < 1.0 - PARABOLIC_ECC_LIMIT => mean_anomaly.rem_euclid(TAU), _ => mean_anomaly, @@ -360,11 +367,11 @@ mod tests { // This was previously a failed instance. let elem = CometElements { desig: Desig::Empty, - epoch: 2461722.5, + epoch: 2461722.5.into(), eccentricity: 0.7495474422690582, inclination: 0.1582845445910239, lon_of_ascending: 1.247985615390004, - peri_time: 2459273.227910867, + peri_time: 2459273.227910867.into(), peri_arg: 4.229481513899533, peri_dist: 0.5613867506855604, }; @@ -374,11 +381,11 @@ mod tests { // This was previously a failed instance. let elem = CometElements { desig: Desig::Empty, - epoch: 2455341.243793971, + epoch: 2455341.243793971.into(), eccentricity: 1.001148327267, inclination: 2.433767, lon_of_ascending: -1.24321, - peri_time: 2454482.5825015577, + peri_time: 2454482.5825015577.into(), peri_arg: 0.823935226897, peri_dist: 5.594792535298549, }; @@ -388,11 +395,11 @@ mod tests { { let elem = CometElements { desig: Desig::Empty, - epoch: 2455562.5, + epoch: 2455562.5.into(), eccentricity: 0.99999, inclination: 2.792526803, lon_of_ascending: 0.349065850, - peri_time: 2455369.7, + peri_time: 2455369.7.into(), peri_arg: -0.8726646259, peri_dist: 0.5, }; @@ -409,11 +416,11 @@ mod tests { for peri_dist in [0.1, 0.5, 10.0] { let elem = CometElements { desig: Desig::Empty, - epoch: 10.0, + epoch: 10.0.into(), eccentricity: ecc, inclination: incl, lon_of_ascending: lon_of_asc, - peri_time: 10.0, + peri_time: 10.0.into(), peri_arg, peri_dist, }; @@ -424,7 +431,7 @@ mod tests { ); let new_elem = CometElements::from_pos_vel( Desig::Empty, - 10.0, + 10.0.into(), &pos.into(), &vel.into(), ); @@ -448,11 +455,11 @@ mod tests { for peri_dist in [0.3, 0.5] { let elem = CometElements { desig: Desig::Empty, - epoch, + epoch: epoch.into(), eccentricity: ecc, inclination: incl, lon_of_ascending: lon_of_asc, - peri_time, + peri_time: peri_time.into(), peri_arg, peri_dist, }; @@ -460,7 +467,7 @@ mod tests { elem.to_pos_vel().expect("Failed to convert to state."); let new_elem = CometElements::from_pos_vel( Desig::Empty, - epoch, + epoch.into(), &pos.into(), &vel.into(), ); @@ -509,18 +516,18 @@ mod tests { for peri_dist in [0.3, 0.5] { let elem = CometElements { desig: Desig::Empty, - epoch, + epoch: epoch.into(), eccentricity: ecc, inclination: incl, lon_of_ascending: lon_of_asc, - peri_time, + peri_time: peri_time.into(), peri_arg, peri_dist, }; let [pos, vel] = elem.to_pos_vel().unwrap(); let new_elem = CometElements::from_pos_vel( Desig::Empty, - epoch, + epoch.into(), &pos.into(), &vel.into(), ); diff --git a/src/kete_core/src/errors.rs b/src/kete_core/src/errors.rs index 41d4cb3..d5d7b27 100644 --- a/src/kete_core/src/errors.rs +++ b/src/kete_core/src/errors.rs @@ -58,7 +58,7 @@ pub enum Error { IOError(String), /// Propagator detected an impact. - Impact(i32, f64), + Impact(i32, Time), } impl error::Error for Error {} @@ -76,6 +76,7 @@ impl fmt::Display for Error { write!(f, "This reference frame is not supported.") } Self::Impact(s, t) => { + let t = t.jd; write!(f, "Propagation detected an impact with {s} at time {t}") } } @@ -85,6 +86,8 @@ impl fmt::Display for Error { #[cfg(feature = "pyo3")] use pyo3::{PyErr, exceptions}; +use crate::time::{TDB, Time}; + #[cfg(feature = "pyo3")] impl From for PyErr { fn from(err: Error) -> Self { @@ -98,9 +101,10 @@ impl From for PyErr { Self::new::("This reference frame is not supported.") } - Error::Impact(s, t) => Self::new::(format!( - "Propagation detected an impact with {s} at time {t}" - )), + Error::Impact(s, t) => Self::new::({ + let t = t.jd; + format!("Propagation detected an impact with {s} at time {t}") + }), } } } diff --git a/src/kete_core/src/fov/fov_like.rs b/src/kete_core/src/fov/fov_like.rs index 687c764..f5be6f9 100644 --- a/src/kete_core/src/fov/fov_like.rs +++ b/src/kete_core/src/fov/fov_like.rs @@ -82,13 +82,13 @@ pub trait FovLike: Sync + Sized { let rel_pos = pos - obs_pos; // This also accounts for first order light delay. - let dt = obs.jd - state.jd - rel_pos.norm() * C_AU_PER_DAY_INV; + let dt = obs.epoch.jd - state.epoch.jd - rel_pos.norm() * C_AU_PER_DAY_INV; let new_pos = pos + vel * dt; let new_rel_pos = new_pos - obs_pos; let (idx, contains) = self.contains(&new_rel_pos); let new_state = State::new( state.desig.clone(), - obs.jd + dt, + obs.epoch + dt, new_pos, vel, obs.center_id, @@ -107,11 +107,11 @@ pub trait FovLike: Sync + Sized { let obs_pos = obs.pos; // bring state up to observer time. - let final_state = propagate_two_body(state, obs.jd)?; + let final_state = propagate_two_body(state, obs.epoch)?; // correct for light delay let dt = -(final_state.pos - obs_pos).norm() * C_AU_PER_DAY_INV; - let final_state = propagate_two_body(&final_state, obs.jd + dt)?; + let final_state = propagate_two_body(&final_state, obs.epoch + dt)?; let rel_pos = final_state.pos - obs_pos; let (idx, contains) = self.contains(&rel_pos); @@ -129,11 +129,11 @@ pub trait FovLike: Sync + Sized { let obs = self.observer(); let obs_pos = obs.pos; - let exact_state = propagate_n_body_spk(state.clone(), obs.jd, include_asteroids, None)?; + let exact_state = propagate_n_body_spk(state.clone(), obs.epoch, include_asteroids, None)?; // correct for light delay let dt = -(exact_state.pos - obs_pos).norm() * C_AU_PER_DAY_INV; - let final_state = propagate_two_body(&exact_state, obs.jd + dt)?; + let final_state = propagate_two_body(&exact_state, obs.epoch + dt)?; let rel_pos = final_state.pos - obs_pos; let (idx, contains) = self.contains(&rel_pos); @@ -180,7 +180,7 @@ pub trait FovLike: Sync + Sized { // to the observer? Then add a factor of 2 for safety let max_dist = (state.vel - obs_state.vel).norm() * dt_limit * 2.0; - if (state.jd - obs_state.jd).abs() < dt_limit { + if (state.epoch - obs_state.epoch).elapsed.abs() < dt_limit { let (_, contains, _) = self.check_linear(state); if let Contains::Outside(dist) = contains && dist > max_dist @@ -234,7 +234,7 @@ pub trait FovLike: Sync + Sized { let states: Vec<_> = obj_ids .into_par_iter() .filter_map(|&obj_id| { - match spk.try_get_state_with_center(obj_id, obs.jd, obs.center_id) { + match spk.try_get_state_with_center(obj_id, obs.epoch, obs.center_id) { Ok(state) => match self.check_two_body(&state) { Ok((idx, Contains::Inside, state)) => Some((idx, state)), _ => None, diff --git a/src/kete_core/src/fov/generic.rs b/src/kete_core/src/fov/generic.rs index bd6e448..9f665fa 100644 --- a/src/kete_core/src/fov/generic.rs +++ b/src/kete_core/src/fov/generic.rs @@ -255,23 +255,23 @@ mod tests { fn test_check_rectangle_visible() { let circular = State::new( Desig::Empty, - 2451545.0, + 2451545.0.into(), [0.0, 1., 0.0].into(), [-GMS_SQRT, 0.0, 0.0].into(), 0, ); let circular_back = State::new( Desig::Empty, - 2451545.0, + 2451545.0.into(), [1.0, 0.0, 0.0].into(), [0.0, GMS_SQRT, 0.0].into(), 0, ); - for offset in [-10.0, -5.0, 0.0, 5.0, 10.0] { + for offset in [-10.0_f64, -5.0, 0.0, 5.0, 10.0] { let off_state = propagate_n_body_spk( circular_back.clone(), - circular_back.jd - offset, + circular_back.epoch - offset, false, None, ) @@ -300,7 +300,7 @@ mod tests { let spk = &LOADED_SPK.read().unwrap(); let observer = State::new( Desig::Empty, - 2451545.0, + 2451545.0.into(), [0.0, 1., 0.0].into(), [-GMS_SQRT, 0.0, 0.0].into(), 10, @@ -308,7 +308,7 @@ mod tests { for offset in [-10.0, -5.0, 0.0, 5.0, 10.0] { let ceres = spk - .try_get_state_with_center(20000001, observer.jd + offset, 10) + .try_get_state_with_center(20000001, observer.epoch + offset, 10) .unwrap(); let fov = OmniDirectional::new(observer.clone()); @@ -318,9 +318,12 @@ mod tests { assert!(two_body.is_ok()); let (_, _, two_body) = two_body.unwrap(); let dist = (two_body.pos - observer.pos).norm(); - assert!((observer.jd - two_body.jd - dist * constants::C_AU_PER_DAY_INV).abs() < 1e-6); + assert!( + (observer.epoch.jd - two_body.epoch.jd - dist * constants::C_AU_PER_DAY_INV).abs() + < 1e-6 + ); let ceres_exact = spk - .try_get_state_with_center(20000001, two_body.jd, 10) + .try_get_state_with_center(20000001, two_body.epoch, 10) .unwrap(); // check that we are within about 150km - not bad for 2 body assert!((two_body.pos - ceres_exact.pos).norm() < 1e-6); @@ -329,9 +332,12 @@ mod tests { let n_body = fov.check_n_body(&ceres, false); assert!(n_body.is_ok()); let (_, _, n_body) = n_body.unwrap(); - assert!((observer.jd - n_body.jd - dist * constants::C_AU_PER_DAY_INV).abs() < 1e-6); + assert!( + (observer.epoch.jd - n_body.epoch.jd - dist * constants::C_AU_PER_DAY_INV).abs() + < 1e-6 + ); let ceres_exact = spk - .try_get_state_with_center(20000001, n_body.jd, 10) + .try_get_state_with_center(20000001, n_body.epoch, 10) .unwrap(); // check that we are within about 150m assert!((n_body.pos - ceres_exact.pos).norm() < 1e-9); @@ -340,9 +346,12 @@ mod tests { let spk_check = &fov.check_spks(&[20000001])[0]; assert!(spk_check.is_some()); let spk_check = &spk_check.as_ref().unwrap().states[0]; - assert!((observer.jd - spk_check.jd - dist * constants::C_AU_PER_DAY_INV).abs() < 1e-6); + assert!( + (observer.epoch.jd - spk_check.epoch.jd - dist * constants::C_AU_PER_DAY_INV).abs() + < 1e-6 + ); let ceres_exact = spk - .try_get_state_with_center(20000001, spk_check.jd, 10) + .try_get_state_with_center(20000001, spk_check.epoch, 10) .unwrap(); // check that we are within about 150 micron assert!((spk_check.pos - ceres_exact.pos).norm() < 1e-12); diff --git a/src/kete_core/src/fov/neos.rs b/src/kete_core/src/fov/neos.rs index ed8062b..5af5108 100644 --- a/src/kete_core/src/fov/neos.rs +++ b/src/kete_core/src/fov/neos.rs @@ -205,7 +205,7 @@ impl NeosVisit { || ccd.subloop_id != subloop_id || ccd.exposure_id != exposure_id || ccd.rotation != rotation - || ccd.observer().jd != observer.jd + || ccd.observer().epoch != observer.epoch || ccd.band != band { Err(Error::ValueError( diff --git a/src/kete_core/src/fov/ptf.rs b/src/kete_core/src/fov/ptf.rs index 0b7534e..262c753 100644 --- a/src/kete_core/src/fov/ptf.rs +++ b/src/kete_core/src/fov/ptf.rs @@ -196,7 +196,7 @@ impl PtfField { let filter = first.filter; for ccd in ccds.iter() { - if ccd.field != field || ccd.filter != filter || ccd.observer().jd != observer.jd { + if ccd.field != field || ccd.filter != filter || ccd.observer().epoch != observer.epoch { Err(Error::ValueError( "All PtfCcds must have matching values except CCD ID etc.".into(), ))?; diff --git a/src/kete_core/src/fov/spherex.rs b/src/kete_core/src/fov/spherex.rs index af61f85..fddfdd0 100644 --- a/src/kete_core/src/fov/spherex.rs +++ b/src/kete_core/src/fov/spherex.rs @@ -139,7 +139,7 @@ impl SpherexField { let observer = first.observer().clone(); for ccd in cmos_frames.iter() { - if ccd.observer().jd != observer.jd { + if ccd.observer().epoch != observer.epoch { Err(Error::ValueError( "All SpherexCMOS must have matching values times".into(), ))?; diff --git a/src/kete_core/src/fov/ztf.rs b/src/kete_core/src/fov/ztf.rs index 0159a1c..c979d75 100644 --- a/src/kete_core/src/fov/ztf.rs +++ b/src/kete_core/src/fov/ztf.rs @@ -179,7 +179,7 @@ impl ZtfField { || ccd.fid != fid || ccd.filtercode != filtercode || ccd.imgtypecode != imgtypecode - || ccd.observer().jd != observer.jd + || ccd.observer().epoch != observer.epoch { Err(Error::ValueError( "All ZtfCcdQuads must have matching values except CCD ID etc.".into(), diff --git a/src/kete_core/src/frames/earth.rs b/src/kete_core/src/frames/earth.rs index 07d7dda..5222e81 100644 --- a/src/kete_core/src/frames/earth.rs +++ b/src/kete_core/src/frames/earth.rs @@ -125,7 +125,7 @@ pub fn geodetic_lat_lon_to_ecef(geodetic_lat: f64, geodetic_lon: f64, h: f64) -> /// pub fn earth_obliquity(jd: Time) -> f64 { // centuries from j2000 - let c = (jd - Time::j2000()) / 365.25 / 100.0; + let c = (jd - Time::j2000()).elapsed / 365.25 / 100.0; (23.439279444444444 + c * (-0.013010213611111 + c * (-5.08611111111111e-08 @@ -175,12 +175,12 @@ pub fn earth_rotation_angle(time: Time) -> f64 { /// /// # Arguments /// -/// * `tdb_time` - Time in TDB scaled Julian Days. +/// * `time` - Time in TDB scaled Julian Days. /// #[inline(always)] pub fn earth_precession_rotation(time: Time) -> NonInertialFrame { // centuries since 2000 - let t = (time - Time::j2000()) / 36525.0; + let t = (time - Time::j2000()).elapsed / 36525.0; // angles as defined in the cited paper, equations (21) // Note that equation 45 is an even more developed model, which takes into @@ -232,7 +232,7 @@ pub fn approx_earth_pos_to_ecliptic( let (pos, vel) = rotation.to_equatorial(pos, [0.0; 3].into())?; - Ok(State::::new(desig, time.tdb().jd, pos.into(), vel.into(), 399).into_frame()) + Ok(State::::new(desig, time.tdb(), pos.into(), vel.into(), 399).into_frame()) } /// Compute the next sunset and sunrise times for a given location. @@ -280,7 +280,7 @@ pub fn next_sunset_sunrise( pub fn approx_sun_dec(time: Time) -> f64 { let obliquity = earth_obliquity(time.tdb()); - let time_since_j2000 = time - Time::j2000(); + let time_since_j2000 = (time - Time::j2000()).elapsed; let mean_lon_of_sun = (280.459 + 0.98564736 * time_since_j2000).rem_euclid(360.0); let mean_anom = (357.529 + 0.98560028 * time_since_j2000) .rem_euclid(360.0) @@ -304,7 +304,7 @@ pub fn approx_sun_dec(time: Time) -> f64 { /// /// pub fn equation_of_time(time: Time) -> f64 { - let time_since_j2000 = time - Time::j2000(); + let time_since_j2000 = (time - Time::j2000()).elapsed; let mean_lon_of_sun = (280.459 + 0.98564736 * time_since_j2000).rem_euclid(360.0); let mean_anom = (357.529 + 0.98560028 * time_since_j2000) .rem_euclid(360.0) diff --git a/src/kete_core/src/io/parquet.rs b/src/kete_core/src/io/parquet.rs index a965e92..b557e14 100644 --- a/src/kete_core/src/io/parquet.rs +++ b/src/kete_core/src/io/parquet.rs @@ -50,7 +50,7 @@ pub fn write_states_parquet(states: &[State], filename: &str) -> Ket ); let jd = Column::new( "jd".into(), - states.iter().map(|state| state.jd).collect_vec(), + states.iter().map(|state| state.epoch.jd).collect_vec(), ); let x = Column::new( "x".into(), @@ -145,7 +145,7 @@ pub fn read_states_parquet(filename: &str) -> KeteResult>> State::new( crate::desigs::Desig::Name(desig.to_string()), - jd, + jd.into(), [x, y, z].into(), [vx, vy, vz].into(), center_id, diff --git a/src/kete_core/src/propagation/acceleration.rs b/src/kete_core/src/propagation/acceleration.rs index e8d1176..b0a29f1 100644 --- a/src/kete_core/src/propagation/acceleration.rs +++ b/src/kete_core/src/propagation/acceleration.rs @@ -51,6 +51,7 @@ use crate::frames::Equatorial; use crate::prelude::KeteResult; use crate::spice::LOADED_SPK; +use crate::time::{TDB, Time}; use crate::{constants, errors::Error, propagation::nongrav::NonGravModel}; use nalgebra::allocator::Allocator; use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, SVector, U1, U2, Vector3}; @@ -60,7 +61,7 @@ use std::ops::AddAssign; #[derive(Debug)] pub struct CentralAccelMeta { /// A vector of times where the central accel function was evaluated at. - pub times: Vec, + pub times: Vec>, /// The position where the central accel function was evaluated. pub pos: Vec>, @@ -95,7 +96,7 @@ impl Default for CentralAccelMeta { /// * `vel` - A vector which defines the velocity with respect to the Sun in AU/Day. /// * `meta` - Metadata object which records values at integration steps. pub fn central_accel( - time: f64, + time: Time, pos: &Vector3, vel: &Vector3, meta: &mut CentralAccelMeta, @@ -115,7 +116,7 @@ pub fn central_accel( pub struct AccelSPKMeta<'a> { /// Closest approach to a massive object. /// This records the ID of the object, time, and distance in AU. - pub close_approach: Option<(i32, f64, f64)>, + pub close_approach: Option<(i32, Time, f64)>, /// The non-gravitational forces. /// If this is not provided, only standard gravitational model is applied. @@ -150,7 +151,7 @@ pub struct AccelSPKMeta<'a> { /// * `vel` - A vector which defines the velocity with respect to the Sun in AU/Day multiplied by `SUN_GMS_SQRT`. /// * `meta` - Metadata object [`AccelSPKMeta`] which records values at each integration step. pub fn spk_accel( - time: f64, + time: Time, pos: &Vector3, vel: &Vector3, meta: &mut AccelSPKMeta<'_>, @@ -158,13 +159,6 @@ pub fn spk_accel( ) -> KeteResult> { let mut accel = Vector3::::zeros(); - if exact_eval - && let Some(close_approach) = meta.close_approach.as_mut() - && close_approach.2 == 0.0 - { - *close_approach = (-1, 1000000.0, 0.0); - } - let spk = &LOADED_SPK.try_read().unwrap(); for grav_params in meta.massive_obj { @@ -177,9 +171,9 @@ pub fn spk_accel( if exact_eval { let r = rel_pos.norm(); if let Some(close_approach) = meta.close_approach.as_mut() - && close_approach.2 >= r + && r <= close_approach.2 { - *close_approach = (id, r, time); + *close_approach = (id, time, r); } if r as f32 <= radius { @@ -204,7 +198,7 @@ pub fn spk_accel( /// The `state_vec` is made up of concatenated position and velocity vectors. /// Otherwise this is just a thin wrapper over the [`spk_accel`] function. pub fn spk_accel_first_order( - time: f64, + time: Time, state_vec: &SVector, meta: &mut AccelSPKMeta<'_>, exact_eval: bool, @@ -246,7 +240,7 @@ pub struct AccelVecMeta<'a> { /// * `vel` - A vector which defines the velocity with respect to the Sun in AU/Day. /// * `meta` - Metadata. pub fn vec_accel( - time: f64, + time: Time, pos: &OVector, vel: &OVector, meta: &mut AccelVecMeta<'_>, @@ -351,7 +345,7 @@ mod tests { #[test] fn check_accelerations_equal() { let spk = &LOADED_SPK.try_read().unwrap(); - let jd = 2451545.0; + let jd = 2451545.0.into(); let mut pos: Vec = Vec::new(); let mut vel: Vec = Vec::new(); diff --git a/src/kete_core/src/propagation/kepler.rs b/src/kete_core/src/propagation/kepler.rs index accb714..7876d43 100644 --- a/src/kete_core/src/propagation/kepler.rs +++ b/src/kete_core/src/propagation/kepler.rs @@ -38,6 +38,7 @@ use crate::fitting::newton_raphson; use crate::frames::InertialFrame; use crate::prelude::{CometElements, KeteResult}; use crate::state::State; +use crate::time::{Duration, TDB, Time}; use argmin::core::{CostFunction, Error as ArgminErr, Executor}; use argmin::solver::neldermead::NelderMead; use core::f64; @@ -283,11 +284,12 @@ fn solve_kepler_universal(mut dt: f64, r0: f64, v0: f64, rv0: f64) -> KeteResult /// * `vel` - Starting velocity, in AU/Day. /// * `depth` - Current Recursion depth, this should be called with `None` initially. pub fn analytic_2_body( - time: f64, + time: Duration, pos: &Vector3, vel: &Vector3, depth: Option, ) -> KeteResult<(Vector3, Vector3)> { + let time = time.elapsed; let mut depth = depth.unwrap_or(0); if depth >= 10 { Err(Error::Convergence( @@ -321,8 +323,9 @@ pub fn analytic_2_body( Ok((new_pos, new_vel)) } Err(_) => { - let (inter_pos, inter_vel) = analytic_2_body(0.5 * time, pos, vel, Some(depth))?; - analytic_2_body(time * 0.5, &inter_pos, &inter_vel, Some(depth)) + let (inter_pos, inter_vel) = + analytic_2_body((0.5 * time).into(), pos, vel, Some(depth))?; + analytic_2_body((time * 0.5).into(), &inter_pos, &inter_vel, Some(depth)) } } } @@ -336,10 +339,10 @@ pub fn analytic_2_body( /// This is the fastest method for getting a relatively good estimate of the orbits. pub fn propagate_two_body( state: &State, - jd_final: f64, + time_final: Time, ) -> KeteResult> { let (pos, vel) = analytic_2_body( - jd_final - state.jd, + time_final - state.epoch, &state.pos.into(), &state.vel.into(), None, @@ -347,7 +350,7 @@ pub fn propagate_two_body( Ok(State::new( state.desig.clone(), - jd_final, + time_final, pos.into(), vel.into(), state.center_id, @@ -368,8 +371,8 @@ impl CostFunction for MoidCost { let dt_a = param.first().unwrap(); let dt_b = param.last().unwrap(); - let s0 = propagate_two_body(&self.state_a, self.state_a.jd + dt_a)?; - let s1 = propagate_two_body(&self.state_b, self.state_b.jd + dt_b)?; + let s0 = propagate_two_body(&self.state_a, (self.state_a.epoch.jd + dt_a).into())?; + let s1 = propagate_two_body(&self.state_b, (self.state_b.epoch.jd + dt_b).into())?; Ok((Vector3::from(s0.pos) - Vector3::from(s1.pos)).norm()) } @@ -400,11 +403,11 @@ pub fn moid(mut state_a: State, mut state_b: State) -> K for idx in (-N_STEPS)..N_STEPS { states_a.push(propagate_two_body( &state_a, - state_a.jd + idx as f64 * state_a_step_size, + (state_a.epoch.jd + idx as f64 * state_a_step_size).into(), )?); states_b.push(propagate_two_body( &state_b, - state_b.jd + idx as f64 * state_b_step_size, + (state_b.epoch.jd + idx as f64 * state_b_step_size).into(), )?); } let mut best = (f64::INFINITY, state_a.clone(), state_b.clone()); @@ -447,12 +450,12 @@ mod tests { let pos = Vector3::new(0.0, r, 0.0); let vel = Vector3::new(-GMS_SQRT / r.sqrt(), 0.0, 0.0); let year = TAU / GMS_SQRT * r.powf(3.0 / 2.0); - let res = analytic_2_body(year, &pos, &vel, None).unwrap(); + let res = analytic_2_body(year.into(), &pos, &vel, None).unwrap(); assert!((res.0 - pos).norm() < 1e-8); assert!((res.1 - vel).norm() < 1e-8); // go backwards - let res = analytic_2_body(-year, &pos, &vel, None).unwrap(); + let res = analytic_2_body((-year).into(), &pos, &vel, None).unwrap(); assert!((res.0 - pos).norm() < 1e-8); assert!((res.1 - vel).norm() < 1e-8); } @@ -474,14 +477,14 @@ mod tests { let pos = Vector3::new(0.0, 2.0, 0.0); let vel = Vector3::new(GMS_SQRT, 0.0, 0.0); let year = -TAU / GMS_SQRT; - let res = analytic_2_body(year, &pos, &vel, None).unwrap(); + let res = analytic_2_body(year.into(), &pos, &vel, None).unwrap(); let pos_exp = Vector3::new(-4.448805955479905, -0.4739843046525608, 0.0); let vel_exp = -Vector3::new(-0.00768983428326951, -0.008552645144187791, 0.0); assert!((res.0 - pos_exp).norm() < 1e-8); assert!((res.1 - vel_exp).norm() < 1e-8); // go backwards - let res = analytic_2_body(-year, &pos_exp, &vel_exp, None).unwrap(); + let res = analytic_2_body((-year).into(), &pos_exp, &vel_exp, None).unwrap(); assert!((res.0 - pos).norm() < 1e-8); assert!((res.1 - vel).norm() < 1e-8); } @@ -491,14 +494,14 @@ mod tests { let pos = Vector3::new(0.0, 3.0, 0.0); let vel = Vector3::new(-GMS_SQRT, 0.0, 0.0); let year = TAU / GMS_SQRT; - let res = analytic_2_body(year, &pos, &vel, None).unwrap(); + let res = analytic_2_body(year.into(), &pos, &vel, None).unwrap(); let pos_exp = Vector3::new(-5.556785268950049, 1.6076633958058089, 0.0); let vel_exp = Vector3::new(-0.013061655543084886, -0.005508140023183166, 0.0); assert!((res.0 - pos_exp).norm() < 1e-8); assert!((res.1 - vel_exp).norm() < 1e-8); // go backwards - let res = analytic_2_body(-year, &pos_exp, &vel_exp, None).unwrap(); + let res = analytic_2_body((-year).into(), &pos_exp, &vel_exp, None).unwrap(); assert!((res.0 - pos).norm() < 1e-8); assert!((res.1 - vel).norm() < 1e-8); } diff --git a/src/kete_core/src/propagation/mod.rs b/src/kete_core/src/propagation/mod.rs index ede93fe..d88d7f8 100644 --- a/src/kete_core/src/propagation/mod.rs +++ b/src/kete_core/src/propagation/mod.rs @@ -38,6 +38,7 @@ use crate::frames::Equatorial; use crate::prelude::{Desig, KeteResult}; use crate::spice::LOADED_SPK; use crate::state::State; +use crate::time::{TDB, Time}; use nalgebra::{DVector, SMatrix, SVector, Vector3}; mod acceleration; @@ -75,8 +76,8 @@ pub fn propagate_two_body_radau(dt: f64, pos: &[f64; 3], vel: &[f64; 3]) -> ([f6 ¢ral_accel, Vector3::from(*pos), Vector3::from(*vel), - 0.0, - dt, + 0.0.into(), + dt.into(), CentralAccelMeta::default(), ) .unwrap(); @@ -87,8 +88,11 @@ pub fn propagate_two_body_radau(dt: f64, pos: &[f64; 3], vel: &[f64; 3]) -> ([f6 /// /// This is a very poor approximation over more than a few minutes/hours, however it /// is very fast. -pub fn propagate_linear(state: &State, jd_final: f64) -> KeteResult> { - let dt = jd_final - state.jd; +pub fn propagate_linear( + state: &State, + jd_final: Time, +) -> KeteResult> { + let dt = (jd_final - state.epoch).elapsed; let mut pos: Vector3 = state.pos.into(); pos.iter_mut() .zip(state.vel) @@ -106,7 +110,7 @@ pub fn propagate_linear(state: &State, jd_final: f64) -> KeteResult< /// Propagate an object using full N-Body physics with the Radau 15th order integrator. pub fn propagate_n_body_spk( mut state: State, - jd_final: f64, + jd_final: Time, include_extended: bool, non_grav_model: Option, ) -> KeteResult> { @@ -133,7 +137,7 @@ pub fn propagate_n_body_spk( &spk_accel, state.pos.into(), state.vel.into(), - state.jd, + state.epoch, jd_final, metadata, )? @@ -147,7 +151,7 @@ pub fn propagate_n_body_spk( /// Initialization function for the picard integrator which initializes the state /// using two body mechanics. fn picard_two_body_init( - times: &[f64; N], + times: &[Time; N], init_pos: &SVector, ) -> SMatrix { let pos: Vector3 = init_pos.fixed_rows::<3>(0).into(); @@ -159,7 +163,7 @@ fn picard_two_body_init( res.fixed_rows_mut::<3>(3).set_column(0, &vel); for (idx, t) in times.iter().enumerate().skip(1) { - let dt = t - t0; + let dt = *t - t0; let (p, v) = analytic_2_body(dt, &pos, &vel, None).unwrap(); res.fixed_rows_mut::<3>(0).set_column(idx, &p); @@ -171,7 +175,7 @@ fn picard_two_body_init( /// Propagate an object using full N-Body physics with the Radau 15th order integrator. pub fn propagate_picard_n_body_spk( mut state: State, - jd_final: f64, + jd_final: Time, include_extended: bool, non_grav_model: Option, ) -> KeteResult> { @@ -208,7 +212,7 @@ pub fn propagate_picard_n_body_spk( &spk_accel_first_order, &picard_two_body_init, state_vec, - state.jd, + state.epoch, jd_final, 1.0, &mut metadata, @@ -229,14 +233,17 @@ pub fn propagate_picard_n_body_spk( /// /// It is *strongly recommended* to use the `kepler.rs` code for this, as /// it will be much more computationally efficient. -pub fn propagation_central(state: &State, jd_final: f64) -> KeteResult<[[f64; 3]; 2]> { +pub fn propagation_central( + state: &State, + jd_final: Time, +) -> KeteResult<[[f64; 3]; 2]> { let pos: Vector3 = state.pos.into(); let vel: Vector3 = state.vel.into(); let (pos, vel, _meta) = RadauIntegrator::integrate( ¢ral_accel, pos, vel, - state.jd, + state.epoch, jd_final, CentralAccelMeta::default(), )?; @@ -248,7 +255,7 @@ pub fn propagation_central(state: &State, jd_final: f64) -> KeteResu #[allow(clippy::type_complexity, reason = "Not practical to avoid this")] pub fn propagate_n_body_vec( states: Vec>, - jd_final: f64, + jd_final: Time, planet_states: Option>>, non_gravs: Vec>, ) -> KeteResult<(Vec>, Vec>)> { @@ -264,7 +271,7 @@ pub fn propagate_n_body_vec( ))?; } - let jd_init = states.first().unwrap().jd; + let jd_init = states.first().unwrap().epoch; let mut pos: Vec = Vec::new(); let mut vel: Vec = Vec::new(); @@ -287,7 +294,7 @@ pub fn propagate_n_body_vec( "Input planet states must contain the correct number of states.".into(), ))?; } - if planet_states.first().unwrap().jd != jd_init { + if planet_states.first().unwrap().epoch != jd_init { Err(Error::ValueError( "Planet states JD must match JD of input state.".into(), ))?; @@ -299,7 +306,7 @@ pub fn propagate_n_body_vec( } for state in states { - if jd_init != state.jd { + if jd_init != state.epoch { Err(Error::ValueError( "All input states must have the same JD".into(), ))?; diff --git a/src/kete_core/src/propagation/nongrav.rs b/src/kete_core/src/propagation/nongrav.rs index 5ad549d..0a89b58 100644 --- a/src/kete_core/src/propagation/nongrav.rs +++ b/src/kete_core/src/propagation/nongrav.rs @@ -175,7 +175,7 @@ impl NonGravModel { let n_vec = t_vec.cross(&pos_norm).normalize(); if !dt.is_zero() { - (pos, _) = analytic_2_body(-dt, &pos, vel, None).unwrap(); + (pos, _) = analytic_2_body((-dt).into(), &pos, vel, None).unwrap(); }; let rr0 = pos.norm() / r_0; let scale = alpha * rr0.powf(-m) * (1.0 + rr0.powf(*n)).powf(-k); diff --git a/src/kete_core/src/propagation/picard.rs b/src/kete_core/src/propagation/picard.rs index 46f1e88..1160daf 100644 --- a/src/kete_core/src/propagation/picard.rs +++ b/src/kete_core/src/propagation/picard.rs @@ -61,14 +61,15 @@ use nalgebra::{SMatrix, SVector}; use crate::errors::{Error, KeteResult}; use crate::propagation::util::FirstOrderODE; +use crate::time::{TDB, Time}; /// Initialization function which fills the integrators initial guess. type PicardInitFunc<'a, const DIM: usize, const N: usize> = - &'a dyn Fn(&[f64; N], &SVector) -> SMatrix; + &'a dyn Fn(&[Time; N], &SVector) -> SMatrix; /// Initialization function for the picard integrator where the first value is repeated pub fn dumb_picard_init( - _times: &[f64; N], + _times: &[Time; N], initial_pos: &SVector, ) -> SMatrix { let mut init: SMatrix = SMatrix::zeros(); @@ -139,10 +140,10 @@ pub struct PicardStep { pub y: SMatrix, /// Start time - pub t0: f64, + pub t0: Time, /// End time - pub t1: f64, + pub t1: Time, } impl<'a, const N: usize, const NM1: usize> PicardIntegrator { @@ -175,11 +176,13 @@ impl<'a, const N: usize, const NM1: usize> PicardIntegrator { func: FirstOrderODE<'a, MType, DIM>, init_func: PicardInitFunc<'a, DIM, N>, initial_pos: SVector, - t0: f64, - t1: f64, + t0: Time, + t1: Time, step_size: f64, metadata: &mut MType, ) -> KeteResult> { + let t0 = t0.jd; + let t1 = t1.jd; let mut cur_t0 = t0; let mut cur_stepsize = step_size; @@ -189,18 +192,18 @@ impl<'a, const N: usize, const NM1: usize> PicardIntegrator { cur_t1 = t1; } - let mut times = self.evaluation_times(cur_t0, cur_t1); + let mut times = self.evaluation_times(cur_t0.into(), cur_t1.into()); let mut initial_pos = init_func(×, &initial_pos); let mut has_failed = false; loop { - let step = self.step(func, cur_t0, cur_t1, initial_pos, metadata); + let step = self.step(func, cur_t0.into(), cur_t1.into(), initial_pos, metadata); match step { Ok(res) => { if cur_t1 == t1 { return Ok(res.y.column(N - 1).into()); } - times = self.evaluation_times(cur_t0, cur_t1); + times = self.evaluation_times(cur_t0.into(), cur_t1.into()); initial_pos = init_func(×, &res.y.column(N - 1).into()); cur_t0 = cur_t1; @@ -240,15 +243,15 @@ impl<'a, const N: usize, const NM1: usize> PicardIntegrator { fn step( &self, func: FirstOrderODE<'a, MType, DIM>, - t0: f64, - t1: f64, + t0: Time, + t1: Time, mut cur_pos: SMatrix, metadata: &mut MType, ) -> KeteResult> { let mut b: SMatrix = SMatrix::zeros(); let mut f: SMatrix = SMatrix::zeros(); let times = self.evaluation_times(t0, t1); - let w2 = (t1 - t0) / 2.0; + let w2 = (t1 - t0).elapsed / 2.0; let mut error: f64 = 100.0; let mut last_error = 50.0; @@ -309,11 +312,14 @@ impl<'a, const N: usize, const NM1: usize> PicardIntegrator { /// Given the start and stop times of a desired integration, return the actual /// times where the function will be evaluated starting at t0 and ending at t1. - fn evaluation_times(&self, t0: f64, t1: f64) -> [f64; N] { - let w2 = (t1 - t0) / 2.0; - let w1 = (t1 + t0) / 2.0; - let mut times = self.tau; - times.iter_mut().for_each(|x| *x = w2 * (*x) + w1); + fn evaluation_times(&self, t0: Time, t1: Time) -> [Time; N] { + let w2 = (t1 - t0).elapsed / 2.0; + let w1 = (t1.jd + t0.jd) / 2.0; + let mut times: [Time; N] = [0.0.into(); N]; + times + .iter_mut() + .zip(self.tau) + .for_each(|(x, tau_v)| x.jd = w2 * tau_v + w1); times } } @@ -322,8 +328,8 @@ impl PicardStep { /// Evaluate the fitted integration solution at the specified time. /// This will fail if the requested time is outside of the integration bounds. pub fn evaluate(&self, t: f64) -> KeteResult<[f64; DIM]> { - let w1 = (self.t0 + self.t1) * 0.5; - let w2 = (self.t1 - self.t0) * 0.5; + let w1 = (self.t0.jd + self.t1.jd) * 0.5; + let w2 = (self.t1 - self.t0).elapsed * 0.5; let tau_time = ((t - w1) * w2).acos(); if tau_time.is_nan() { return Err(Error::ExceedsLimits( @@ -445,7 +451,7 @@ mod tests { // define the f' function fn func( - _: f64, + _: Time, vals: &SVector, _: &mut (), _: bool, @@ -456,13 +462,13 @@ mod tests { // integrate with the initial conditions let p = &PC15; - let t1 = 5.0; + let t1 = 5.0.into(); let res = p .integrate( &func, &dumb_picard_init, [1.0, 2.0, 5.0].into(), - 0.0, + 0.0.into(), t1, 0.1, &mut (), @@ -470,12 +476,12 @@ mod tests { .unwrap(); // test against the analytic solution - assert!((res[0] - (1.0 * (-t1).exp())).abs() < 1e-15); - assert!((res[1] - (2.0 * (-0.5 * t1).exp())).abs() < 1e-15); + assert!((res[0] - (1.0 * (-t1.jd).exp())).abs() < 1e-15); + assert!((res[1] - (2.0 * (-0.5 * t1.jd).exp())).abs() < 1e-15); assert!( - (res[2] - (5.0 * (-0.1 * t1).exp())).abs() < 1e-14, + (res[2] - (5.0 * (-0.1 * t1.jd).exp())).abs() < 1e-14, "{}", - (res[2] - (5.0 * (-0.1 * t1).exp())).abs() + (res[2] - (5.0 * (-0.1 * t1.jd).exp())).abs() ); } } diff --git a/src/kete_core/src/propagation/radau.rs b/src/kete_core/src/propagation/radau.rs index 2d40df1..bde6783 100644 --- a/src/kete_core/src/propagation/radau.rs +++ b/src/kete_core/src/propagation/radau.rs @@ -32,6 +32,7 @@ use crate::errors::Error; use crate::prelude::KeteResult; use crate::propagation::util::SecondOrderODE; +use crate::time::{TDB, Time}; use itertools::izip; use nalgebra::allocator::Allocator; use nalgebra::*; @@ -112,9 +113,9 @@ where func: SecondOrderODE<'a, MType, D>, metadata: MType, - final_time: f64, + final_time: Time, - cur_time: f64, + cur_time: Time, cur_state: OVector, cur_state_der: OVector, cur_state_der_der: OVector, @@ -136,8 +137,8 @@ where func: SecondOrderODE<'a, MType, D>, state_init: OVector, state_der_init: OVector, - time_init: f64, - final_time: f64, + time_init: Time, + final_time: Time, metadata: MType, ) -> KeteResult { let (dim, _) = state_init.shape_generic(); @@ -177,8 +178,8 @@ where func: SecondOrderODE<'a, MType, D>, state_init: OVector, state_der_init: OVector, - time_init: f64, - final_time: f64, + time_init: Time, + final_time: Time, metadata: MType, ) -> RadauResult { let mut integrator = Self::new( @@ -189,24 +190,25 @@ where final_time, metadata, )?; - if (final_time - time_init).abs() < 1e-10 { + if (final_time - time_init).elapsed.abs() < 1e-10 { return Ok(( integrator.cur_state, integrator.cur_state_der, integrator.metadata, )); } - let mut next_step_size: f64 = 0.1_f64.copysign(integrator.final_time - integrator.cur_time); + let mut next_step_size: f64 = + 0.1_f64.copysign((integrator.final_time - integrator.cur_time).elapsed); let mut step_failures = 0; loop { - if (integrator.cur_time - integrator.final_time).abs() <= next_step_size.abs() { - next_step_size = integrator.final_time - integrator.cur_time; + if (integrator.cur_time - integrator.final_time).elapsed.abs() <= next_step_size.abs() { + next_step_size = (integrator.final_time - integrator.cur_time).elapsed; } match integrator.step(next_step_size) { Ok(s) => { next_step_size = s; - if (integrator.cur_time - integrator.final_time).abs() < 1e-12 { + if (integrator.cur_time - integrator.final_time).elapsed.abs() < 1e-12 { return Ok(( integrator.cur_state, integrator.cur_state_der, @@ -297,7 +299,7 @@ where self.eval_scratch.set_column( 0, &(self.func)( - self.cur_time + gauss_radau_frac * step_size, + (self.cur_time.jd + gauss_radau_frac * step_size).into(), &self.state_scratch, &self.state_der_scratch, &mut self.metadata, @@ -344,7 +346,7 @@ where + self.cur_b.row(idx).dot(&U_VEC)); } } - self.cur_time += step_size; + self.cur_time.jd += step_size; self.cur_state_der_der = (self.func)( self.cur_time, &self.cur_state, @@ -375,8 +377,8 @@ mod tests { ¢ral_accel, Vector3::new(0.46937657, -0.8829981, 0.), Vector3::new(0.01518942, 0.00807426, 0.), - 0.0, - 1000., + 0.0.into(), + 1000.0.into(), CentralAccelMeta::default(), ) .unwrap(); diff --git a/src/kete_core/src/propagation/runge_kutta.rs b/src/kete_core/src/propagation/runge_kutta.rs index c6e92bd..e5853a5 100644 --- a/src/kete_core/src/propagation/runge_kutta.rs +++ b/src/kete_core/src/propagation/runge_kutta.rs @@ -27,15 +27,14 @@ // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -use crate::{errors::Error, prelude::KeteResult}; +use crate::{ + errors::Error, + prelude::KeteResult, + propagation::util::FirstOrderODE, + time::{TDB, Time}, +}; use nalgebra::SVector; -/// Function will be of the form `y' = F(t, y, metadata, exact_eval)`. -/// Metadata is passed for every evaluation. The `exact_eval` bool indicates to the -/// function that the input parameters are known to be solutions for the IVP. -type FirstOrderFunc<'a, MType, const D: usize> = - &'a dyn Fn(f64, &SVector, &mut MType, bool) -> KeteResult>; - /// Integrator will return a result of this type. type FirstOrderResult = KeteResult<(SVector, MType)>; @@ -43,12 +42,12 @@ type FirstOrderResult = KeteResult<(SVector, MTyp /// #[allow(missing_debug_implementations, reason = "Not practical to avoid this")] pub struct RK45Integrator<'a, MType, const D: usize> { - func: FirstOrderFunc<'a, MType, D>, + func: FirstOrderODE<'a, MType, D>, metadata: MType, - final_time: f64, + final_time: Time, - cur_time: f64, + cur_time: Time, cur_state: SVector, cur_der: SVector, @@ -132,7 +131,7 @@ impl<'a, MType, const D: usize> RK45Integrator<'a, MType, D> { + k4 * (28561.0 / 56430.0) - k5 * (9.0 / 50.0) + k6 * (2.0 / 55.0); - self.cur_time += step_size; + self.cur_time.jd += step_size; self.cur_der = (self.func)(self.cur_time, &self.cur_state, &mut self.metadata, true)?; Ok(new_step) @@ -143,14 +142,14 @@ impl<'a, MType, const D: usize> RK45Integrator<'a, MType, D> { /// /// This is a relatively fast, but not very precise numerical integration method. pub fn integrate( - func: FirstOrderFunc<'a, MType, D>, + func: FirstOrderODE<'a, MType, D>, state_init: SVector, - time_init: f64, - final_time: f64, + time_init: Time, + final_time: Time, mut metadata: MType, tol: f64, ) -> FirstOrderResult { - if (time_init - final_time).abs() < 1e-30 { + if (time_init - final_time).elapsed.abs() < 1e-30 { return Ok((state_init, metadata)); } let cur_der = func(time_init, &state_init, &mut metadata, true)?; @@ -163,11 +162,11 @@ impl<'a, MType, const D: usize> RK45Integrator<'a, MType, D> { cur_der, tol, }; - let mut next_step = final_time - integrator.cur_time; + let mut next_step = (final_time - integrator.cur_time).elapsed; next_step = integrator.step(next_step)?; - while (integrator.cur_time - integrator.final_time).abs() > tol { - if (integrator.cur_time - integrator.final_time).abs() < next_step.abs() { - next_step = integrator.final_time - integrator.cur_time; + while (integrator.cur_time - integrator.final_time).elapsed.abs() > tol { + if (integrator.cur_time - integrator.final_time).elapsed.abs() < next_step.abs() { + next_step = (integrator.final_time - integrator.cur_time).elapsed; } next_step = integrator.step(next_step)?; } @@ -182,11 +181,12 @@ mod tests { use crate::prelude::KeteResult; use crate::propagation::RK45Integrator; + use crate::time::{TDB, Time}; #[test] fn basic_exp() { fn f( - _t: f64, + _t: Time, state: &Vector1, _meta: &mut (), _eval: bool, @@ -198,8 +198,8 @@ mod tests { let res = RK45Integrator::integrate( &f, Vector1::::repeat(1_f64), - 0.0, - step as f64 * 10.0, + 0.0.into(), + (step as f64 * 10.0).into(), (), 1e-14, ); @@ -211,7 +211,7 @@ mod tests { #[test] fn basic_line() { fn f( - _t: f64, + _t: Time, _state: &Vector1, _meta: &mut (), _eval: bool, @@ -219,8 +219,14 @@ mod tests { Ok(Vector1::::repeat(2.0)) } - let res = - RK45Integrator::integrate(&f, Vector1::::repeat(1_f64), 0.0, 10.0, (), 0.001); + let res = RK45Integrator::integrate( + &f, + Vector1::::repeat(1_f64), + 0.0.into(), + 10.0.into(), + (), + 0.001, + ); assert!(res.is_ok()); assert!((res.unwrap().0[0] - 21.0).abs() < 1e-12); diff --git a/src/kete_core/src/propagation/state_transition.rs b/src/kete_core/src/propagation/state_transition.rs index f2b59eb..5083849 100644 --- a/src/kete_core/src/propagation/state_transition.rs +++ b/src/kete_core/src/propagation/state_transition.rs @@ -33,10 +33,11 @@ use crate::prelude::{KeteResult, State}; use crate::propagation::{ CentralAccelMeta, PC15, central_accel, central_accel_grad, dumb_picard_init, }; +use crate::time::{TDB, Time}; use nalgebra::{Const, Matrix6, SVector, U1, U6, Vector3}; fn stm_ivp_eqn( - jd: f64, + jd: Time, state: &SVector, meta: &mut CentralAccelMeta, exact_eval: bool, @@ -83,7 +84,7 @@ fn stm_ivp_eqn( /// This uses the Picard-Chebyshev integrator [`PC15`]. pub fn compute_state_transition( state: &mut State, - jd: f64, + jd: Time, central_mass: f64, ) -> ([[f64; 3]; 2], Matrix6) { let mut meta = CentralAccelMeta { @@ -111,8 +112,8 @@ pub fn compute_state_transition( &stm_ivp_eqn, &dumb_picard_init, initial_state, - state.jd * GMS_SQRT, - jd * GMS_SQRT, + (state.epoch.jd * GMS_SQRT).into(), + (jd.jd * GMS_SQRT).into(), 1.0, &mut meta, ) diff --git a/src/kete_core/src/propagation/util.rs b/src/kete_core/src/propagation/util.rs index a70189d..9386163 100644 --- a/src/kete_core/src/propagation/util.rs +++ b/src/kete_core/src/propagation/util.rs @@ -1,16 +1,19 @@ use nalgebra::{OVector, SVector}; -use crate::errors::KeteResult; +use crate::{ + errors::KeteResult, + time::{TDB, Time}, +}; /// Function will be of the form y' = F(time, y, metadata, bool) /// This is the first-order general IVP solver. pub(crate) type FirstOrderODE<'a, MType, const DIM: usize> = - &'a dyn Fn(f64, &SVector, &mut MType, bool) -> KeteResult>; + &'a dyn Fn(Time, &SVector, &mut MType, bool) -> KeteResult>; /// Function will be of the form y'' = F(time, y, y', metadata, bool) /// This is the input for a second-order general IVP solver. pub(crate) type SecondOrderODE<'a, MType, D> = &'a dyn Fn( - f64, + Time, &OVector, &OVector, &mut MType, diff --git a/src/kete_core/src/simult_states.rs b/src/kete_core/src/simult_states.rs index a37dd65..bcd36b1 100644 --- a/src/kete_core/src/simult_states.rs +++ b/src/kete_core/src/simult_states.rs @@ -35,6 +35,7 @@ use crate::fov::FOV; use crate::frames::{Equatorial, Vector}; use crate::io::FileIO; use crate::prelude::{Error, KeteResult, State}; +use crate::time::{TDB, Time}; use nalgebra::Vector3; use rayon::iter::{IndexedParallelIterator, IntoParallelRefIterator, ParallelIterator}; use serde::{Deserialize, Serialize}; @@ -46,7 +47,7 @@ pub struct SimultaneousStates { pub states: Vec>, /// Common JD time of all states - pub jd: f64, + pub epoch: Time, /// Center ID of all states. pub center_id: i32, @@ -69,18 +70,18 @@ impl SimultaneousStates { } let (mut jd, center_id) = { let first = states.first().unwrap(); - (first.jd, first.center_id) + (first.epoch, first.center_id) }; if let Some(f) = &fov { - jd = f.observer().jd; + jd = f.observer().epoch; } if states.iter().any(|state| state.center_id != center_id) { return Err(Error::ValueError("Center IDs do not match expected".into())); }; - if fov.is_none() && states.iter().any(|state| state.jd != jd) { + if fov.is_none() && states.iter().any(|state| state.epoch != jd) { return Err(Error::ValueError( "Epoch JDs do not match expected, this is only allowed if there is an associated FOV." .into(), @@ -89,7 +90,7 @@ impl SimultaneousStates { Ok(Self { states, - jd, + epoch: jd, center_id, fov, }) diff --git a/src/kete_core/src/spice/daf.rs b/src/kete_core/src/spice/daf.rs index 9f6a876..516f4f3 100644 --- a/src/kete_core/src/spice/daf.rs +++ b/src/kete_core/src/spice/daf.rs @@ -45,6 +45,7 @@ use crate::io::bytes::{ }; use crate::errors::{Error, KeteResult}; +use crate::time::{TDB, Time}; use std::fmt::Debug; use std::io::{Cursor, Read, Seek}; use std::ops::Index; @@ -359,7 +360,7 @@ pub struct SpkArray { impl SpkArray { /// Is the specified JD within the range of this array. - pub fn contains(&self, jd: f64) -> bool { + pub fn contains(&self, jd: Time) -> bool { let jds = jd_to_spice_jd(jd); (jds >= self.jds_start) && (jds <= self.jds_end) } @@ -430,7 +431,7 @@ pub struct PckArray { impl PckArray { /// Is the specified JD within the range of this array. - pub fn contains(&self, jd: f64) -> bool { + pub fn contains(&self, jd: Time) -> bool { let jds = jd_to_spice_jd(jd); (jds >= self.jds_start) && (jds <= self.jds_end) } diff --git a/src/kete_core/src/spice/mod.rs b/src/kete_core/src/spice/mod.rs index f90dc19..352d781 100644 --- a/src/kete_core/src/spice/mod.rs +++ b/src/kete_core/src/spice/mod.rs @@ -52,6 +52,8 @@ pub use pck::{LOADED_PCK, PckCollection}; pub use sclk::{LOADED_SCLK, SclkCollection}; pub use spk::{LOADED_SPK, SpkCollection}; +use crate::time::{TDB, Time}; + /// Convert seconds from J2000 into JD. /// /// # Arguments @@ -60,16 +62,16 @@ pub use spk::{LOADED_SPK, SpkCollection}; /// # Returns /// The Julian Date (TDB). #[inline(always)] -fn spice_jd_to_jd(jds_sec: f64) -> f64 { +fn spice_jd_to_jd(jds_sec: f64) -> Time { // 86400.0 = 60 * 60 * 24 - jds_sec / 86400.0 + 2451545.0 + (jds_sec / 86400.0 + 2451545.0).into() } /// Convert TDB JD to seconds from J2000. #[inline(always)] -fn jd_to_spice_jd(jd: f64) -> f64 { +fn jd_to_spice_jd(epoch: Time) -> f64 { // 86400.0 = 60 * 60 * 24 - (jd - 2451545.0) * 86400.0 + (epoch.jd - 2451545.0) * 86400.0 } #[cfg(test)] @@ -81,24 +83,24 @@ mod tests { { let jd_sec = 0.0; let jd = spice_jd_to_jd(jd_sec); - assert_eq!(jd, 2451545.0); + assert_eq!(jd, 2451545.0.into()); } { let jd_sec = 86400.0; // 1 day in seconds let jd = spice_jd_to_jd(jd_sec); - assert_eq!(jd, 2451546.0); + assert_eq!(jd, 2451546.0.into()); } } #[test] fn test_jd_to_spice_jd() { { - let jd = 2451545.0; + let jd = 2451545.0.into(); let jd_sec = jd_to_spice_jd(jd); assert_eq!(jd_sec, 0.0); } { - let jd = 2451546.0; // 1 day after J2000 + let jd = 2451546.0.into(); // 1 day after J2000 let jd_sec = jd_to_spice_jd(jd); assert_eq!(jd_sec, 86400.0); } diff --git a/src/kete_core/src/spice/pck.rs b/src/kete_core/src/spice/pck.rs index 9b93465..94df098 100644 --- a/src/kete_core/src/spice/pck.rs +++ b/src/kete_core/src/spice/pck.rs @@ -44,6 +44,7 @@ use super::pck_segments::PckSegment; use crate::cache::cache_path; use crate::errors::{Error, KeteResult}; use crate::frames::NonInertialFrame; +use crate::time::{TDB, Time}; use crossbeam::sync::ShardedLock; /// A collection of segments. @@ -74,7 +75,7 @@ impl PckCollection { /// Get the raw orientation from the loaded PCK files. /// This orientation will have the frame of what was originally present in the file. - pub fn try_get_orientation(&self, id: i32, jd: f64) -> KeteResult { + pub fn try_get_orientation(&self, id: i32, jd: Time) -> KeteResult { for segment in self.segments.iter() { let array: &PckArray = segment.into(); if (array.frame_id == id) & array.contains(jd) { diff --git a/src/kete_core/src/spice/pck_segments.rs b/src/kete_core/src/spice/pck_segments.rs index daff0bc..145290a 100644 --- a/src/kete_core/src/spice/pck_segments.rs +++ b/src/kete_core/src/spice/pck_segments.rs @@ -47,6 +47,7 @@ use super::{PckArray, interpolation::*}; use crate::errors::Error; use crate::frames::NonInertialFrame; use crate::prelude::KeteResult; +use crate::spice::spice_jd_to_jd; use crate::time::{TDB, Time}; use std::fmt::Debug; @@ -90,7 +91,7 @@ impl PckSegment { pub(in crate::spice) fn try_get_orientation( &self, center_id: i32, - jd: f64, + epoch: Time, ) -> KeteResult { let arr_ref: &PckArray = self.into(); @@ -100,7 +101,7 @@ impl PckSegment { ))?; } - let jds = jd_to_spice_jd(jd); + let jds = jd_to_spice_jd(epoch); if jds < arr_ref.jds_start || jds > arr_ref.jds_end { Err(Error::ExceedsLimits( @@ -174,7 +175,7 @@ impl PckSegmentType2 { // rem_euclid is equivalent to the modulo operator, so this maps w to [0, 2pi] let w = w.rem_euclid(std::f64::consts::TAU); - let time = Time::::new(jd_to_spice_jd(jds)); + let time = spice_jd_to_jd(jds); let frame = NonInertialFrame::from_euler::<'Z', 'X', 'Z'>( time, [ra, dec, w], diff --git a/src/kete_core/src/spice/sclk.rs b/src/kete_core/src/spice/sclk.rs index 10c4e07..4536c6c 100644 --- a/src/kete_core/src/spice/sclk.rs +++ b/src/kete_core/src/spice/sclk.rs @@ -256,13 +256,12 @@ impl Sclk { * (clock_rate[2] / (*self.tick_rates.first().unwrap() as f64)) + clock_rate[1]; - Ok(Time::::new(spice_jd_to_jd(par_time))) + Ok(spice_jd_to_jd(par_time)) } /// Convert time in TDB to a spacecraft clock tick count. fn time_to_tick(&self, time: Time) -> KeteResult { - let jd = time.jd; - let par_time = jd_to_spice_jd(jd); + let par_time = jd_to_spice_jd(time); let clock_rate = self.find_parallel_time_rate(par_time)?; let tick = (par_time - clock_rate[1]) diff --git a/src/kete_core/src/spice/spk.rs b/src/kete_core/src/spice/spk.rs index 60d0b1e..5d97790 100644 --- a/src/kete_core/src/spice/spk.rs +++ b/src/kete_core/src/spice/spk.rs @@ -56,6 +56,7 @@ use crate::errors::Error; use crate::frames::InertialFrame; use crate::prelude::KeteResult; use crate::state::State; +use crate::time::{TDB, Time}; use pathfinding::prelude::dijkstra; use std::collections::{HashMap, HashSet}; use std::fs; @@ -92,7 +93,7 @@ impl SpkCollection { /// This state will have the center and frame of whatever was originally loaded /// into the file. #[inline(always)] - pub fn try_get_state(&self, id: i32, jd: f64) -> KeteResult> { + pub fn try_get_state(&self, id: i32, jd: Time) -> KeteResult> { for segment in &self.planet_segments { let arr_ref: &SpkArray = segment.into(); if arr_ref.object_id == id && arr_ref.contains(jd) { @@ -118,7 +119,7 @@ impl SpkCollection { pub fn try_get_state_with_center( &self, id: i32, - jd: f64, + jd: Time, center: i32, ) -> KeteResult> { let mut state = self.try_get_state(id, jd)?; @@ -137,24 +138,24 @@ impl SpkCollection { match (state.center_id, new_center) { (a, b) if a == b => (), (i, 0) if i <= 10 => { - state.try_change_center(self.try_get_state(i, state.jd)?)?; + state.try_change_center(self.try_get_state(i, state.epoch)?)?; } (0, 10) => { - let next = self.try_get_state(10, state.jd)?; + let next = self.try_get_state(10, state.epoch)?; state.try_change_center(next)?; } (i, 10) if i < 10 => { - state.try_change_center(self.try_get_state(i, state.jd)?)?; - state.try_change_center(self.try_get_state(10, state.jd)?)?; + state.try_change_center(self.try_get_state(i, state.epoch)?)?; + state.try_change_center(self.try_get_state(10, state.epoch)?)?; } (10, i) if (i > 1) & (i < 10) => { - state.try_change_center(self.try_get_state(10, state.jd)?)?; - state.try_change_center(self.try_get_state(i, state.jd)?)?; + state.try_change_center(self.try_get_state(10, state.epoch)?)?; + state.try_change_center(self.try_get_state(i, state.epoch)?)?; } _ => { let path = self.find_path(state.center_id, new_center)?; for intermediate in path { - let next = self.try_get_state(intermediate, state.jd)?; + let next = self.try_get_state(intermediate, state.epoch)?; state.try_change_center(next)?; } } @@ -163,8 +164,8 @@ impl SpkCollection { } /// For a given NAIF ID, return all increments of time which are currently loaded. - pub fn available_info(&self, id: i32) -> Vec<(f64, f64, i32, i32, i32)> { - let mut segment_info = Vec::<(f64, f64, i32, i32, i32)>::new(); + pub fn available_info(&self, id: i32) -> Vec<(Time, Time, i32, i32, i32)> { + let mut segment_info = Vec::<(Time, Time, i32, i32, i32)>::new(); if let Some(segments) = self.segments.get(&id) { for segment in segments { let spk_array_ref: &SpkArray = segment.into(); @@ -198,19 +199,19 @@ impl SpkCollection { return segment_info; } - segment_info.sort_by(|a, b| (a.0).total_cmp(&b.0)); + segment_info.sort_by(|a, b| (a.0.jd).total_cmp(&b.0.jd)); - let mut avail_times = Vec::<(f64, f64, i32, i32, i32)>::new(); + let mut avail_times = Vec::<(Time, Time, i32, i32, i32)>::new(); let mut cur_segment = segment_info[0]; for segment in segment_info.iter().skip(1) { // if the segments are overlapped or nearly overlapped, join them together // 1e-8 is approximately a millisecond - if cur_segment.1 <= (segment.0 - 1e-8) { + if cur_segment.1.jd <= (segment.0.jd - 1e-8) { avail_times.push(cur_segment); cur_segment = *segment; } else { - cur_segment.1 = segment.1.max(cur_segment.1); + cur_segment.1.jd = segment.1.jd.max(cur_segment.1.jd); } } avail_times.push(cur_segment); diff --git a/src/kete_core/src/spice/spk_segments.rs b/src/kete_core/src/spice/spk_segments.rs index 8ba5e18..1a4839e 100644 --- a/src/kete_core/src/spice/spk_segments.rs +++ b/src/kete_core/src/spice/spk_segments.rs @@ -109,7 +109,7 @@ impl SpkSegment { #[inline(always)] pub(in crate::spice) fn try_get_state( &self, - jd: f64, + jd: Time, ) -> KeteResult> { let arr_ref: &SpkArray = self.into(); @@ -640,7 +640,7 @@ impl SpkSegmentType10 { ] = *rec; let epoch = julian_years_since_j2000( - &Time::::new(spice_jd_to_jd(epoch)) + &spice_jd_to_jd(epoch) .utc() .to_datetime() .unwrap() diff --git a/src/kete_core/src/state.rs b/src/kete_core/src/state.rs index 3f28bf3..bdaa8a9 100644 --- a/src/kete_core/src/state.rs +++ b/src/kete_core/src/state.rs @@ -52,6 +52,7 @@ use std::fmt::Debug; use crate::desigs::Desig; use crate::errors::{Error, KeteResult}; use crate::frames::{InertialFrame, Vector}; +use crate::time::{TDB, Time}; /// Exact State of an object. /// @@ -67,8 +68,8 @@ where /// Designation number which corresponds to the object. pub desig: Desig, - /// JD of the object's state in TDB scaled time. - pub jd: f64, + /// Epoch of the object's state in TDB scaled time. + pub epoch: Time, /// Position of the object with respect to the `center_id` object, units of AU. pub pos: Vector, @@ -84,10 +85,16 @@ where impl State { /// Construct a new State object. #[inline(always)] - pub fn new(desig: Desig, jd: f64, pos: Vector, vel: Vector, center_id: i32) -> Self { + pub fn new( + desig: Desig, + epoch: Time, + pos: Vector, + vel: Vector, + center_id: i32, + ) -> Self { Self { desig, - jd, + epoch, pos, vel, center_id, @@ -98,13 +105,19 @@ impl State { /// remaining data. This is primarily useful as a place holder when propagation /// has failed and the object needs to be recorded still. #[inline(always)] - pub fn new_nan(desig: Desig, jd: f64, center_id: i32) -> Self { - Self::new(desig, jd, Vector::new_nan(), Vector::new_nan(), center_id) + pub fn new_nan(desig: Desig, epoch: Time, center_id: i32) -> Self { + Self::new( + desig, + epoch, + Vector::new_nan(), + Vector::new_nan(), + center_id, + ) } /// Are all values finite. pub fn is_finite(&self) -> bool { - self.pos.is_finite() & self.vel.is_finite() & self.jd.is_finite() + self.pos.is_finite() & self.vel.is_finite() & self.epoch.jd.is_finite() } /// Trade the center ID and ID values, and flip the direction of the position and @@ -136,8 +149,10 @@ impl State { /// * `state` - [`State`] object which defines the new center point. #[inline(always)] pub fn try_change_center(&mut self, mut state: Self) -> KeteResult<()> { - if self.jd != state.jd { - return Err(Error::ValueError("States don't have matching jds.".into())); + if self.epoch != state.epoch { + return Err(Error::ValueError( + "States don't have matching epochs.".into(), + )); } let Desig::Naif(state_id) = state.desig else { @@ -178,7 +193,7 @@ impl State { State { desig: self.desig, - jd: self.jd, + epoch: self.epoch, pos, vel, center_id: self.center_id, @@ -197,7 +212,7 @@ mod tests { fn flip_center() { let mut a = State::::new( Desig::Naif(1), - 0.0, + 0.0.into(), [1.0, 0.0, 0.0].into(), [0.0, 1.0, 0.0].into(), 0, @@ -215,14 +230,14 @@ mod tests { fn nan_finite() { let a = State::::new( Desig::Naif(1), - 0.0, + 0.0.into(), [1.0, 0.0, 0.0].into(), [0.0, 1.0, 0.0].into(), 0, ); assert!(a.is_finite()); - let b = State::::new_nan(Desig::Empty, 0.0, 1000); + let b = State::::new_nan(Desig::Empty, 0.0.into(), 1000); assert!(!b.is_finite()); } @@ -230,7 +245,7 @@ mod tests { fn naif_name_resolution() { let mut a = State::::new( Desig::Naif(1), - 0.0, + 0.0.into(), [1.0, 0.0, 0.0].into(), [0.0, 1.0, 0.0].into(), 0, @@ -245,14 +260,14 @@ mod tests { fn change_center() { let mut a = State::::new( Desig::Naif(1), - 0.0, + 0.0.into(), [1.0, 0.0, 0.0].into(), [1.0, 0.0, 0.0].into(), 0, ); let b = State::::new( Desig::Naif(3), - 0.0, + 0.0.into(), [0.0, 1.0, 0.0].into(), [0.0, 1.0, 0.0].into(), 0, @@ -268,7 +283,7 @@ mod tests { // try cases which cause errors let diff_jd = State::::new( Desig::Naif(3), - 1.0, + 1.0.into(), [0.0, 1.0, 0.0].into(), [0.0, 1.0, 0.0].into(), 0, @@ -277,7 +292,7 @@ mod tests { let not_naif_id = State::::new( Desig::Empty, - 0.0, + 0.0.into(), [0.0, 1.0, 0.0].into(), [0.0, 1.0, 0.0].into(), 0, @@ -286,7 +301,7 @@ mod tests { let no_matching_id = State::::new( Desig::Naif(2), - 0.0, + 0.0.into(), [0.0, 1.0, 0.0].into(), [0.0, 1.0, 0.0].into(), 1000000000, diff --git a/src/kete_core/src/time/mod.rs b/src/kete_core/src/time/mod.rs index 27dc927..7d68a9c 100644 --- a/src/kete_core/src/time/mod.rs +++ b/src/kete_core/src/time/mod.rs @@ -41,6 +41,7 @@ mod leap_second; mod scales; use chrono::{DateTime, Datelike, NaiveDate, Timelike, Utc}; +use serde::{Deserialize, Serialize}; use crate::prelude::{Error, KeteResult}; @@ -67,7 +68,7 @@ pub use self::scales::{JD_TO_MJD, TAI, TDB, TimeScale, UTC}; /// Any conversions to a single float will by necessity result in some small accuracy /// loss due to the nature of the representation of numbers on computers. /// -#[derive(Debug, Clone, Copy, PartialEq)] +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] pub struct Time { /// Julian Date @@ -77,6 +78,41 @@ pub struct Time { scale_type: PhantomData, } +/// Duration of time. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Duration { + /// Elapsed time in days. + pub elapsed: f64, + + /// [`PhantomData`] is used here as the scale is only a record keeping convenience. + scale_type: PhantomData, +} + +impl Duration { + /// Construct a new [`Duration`] object. + pub fn new(elapsed: f64) -> Self { + Self { + elapsed, + scale_type: PhantomData, + } + } + + /// Cast to UTC scaled time. + pub fn utc(&self) -> Duration { + Duration::::new(UTC::from_tdb(T::to_tdb(self.elapsed))) + } + + /// Cast to TAI scaled time. + pub fn tai(&self) -> Duration { + Duration::::new(TAI::from_tdb(T::to_tdb(self.elapsed))) + } + + /// Cast to TDB scaled time. + pub fn tdb(&self) -> Duration { + Duration::::new(T::to_tdb(self.elapsed)) + } +} + /// Convert Hours/Minutes/Seconds/millisecond to fraction of a day pub fn hour_min_sec_to_day(h: u32, m: u32, s: u32, ms: u32) -> f64 { h as f64 / 24.0 @@ -302,21 +338,54 @@ impl From for Time { } } +impl From for Duration { + fn from(value: f64) -> Self { + Self::new(value) + } +} + impl Sub> for Time { - type Output = f64; + type Output = Duration; - /// Subtract two times, returning the difference in days. + /// Subtract two times, returning the difference in days of TDB. fn sub(self, other: Time) -> Self::Output { - self.tdb().jd - other.tdb().jd + (self.tdb().jd - other.tdb().jd).into() + } +} + +impl Add> for Time { + type Output = Self; + + /// Add a time and duration together + fn add(self, other: Duration) -> Self::Output { + A::from_tdb(self.tdb().jd + other.tdb().elapsed).into() + } +} + +impl Add for Time { + type Output = Self; + + /// Add a time and duration together + fn add(self, other: f64) -> Self::Output { + A::from_tdb(self.tdb().jd + other).into() + } +} + +impl Sub> for Time { + type Output = Self; + + /// Subtract two times, returning the difference in days. + fn sub(self, other: Duration) -> Self::Output { + A::from_tdb(self.tdb().jd - other.tdb().elapsed).into() } } -impl Add> for Time { +impl Sub for Time { type Output = Self; - /// Add two times together, returning a new time in the same scale as `A`. - fn add(self, other: Time) -> Self { - Self::new(A::from_tdb(self.tdb().jd + other.tdb().jd)) + /// Add a time and duration together + fn sub(self, other: f64) -> Self::Output { + A::from_tdb(self.tdb().jd - other).into() } } From 8530a6ca44db09ef7443265122728dec05da93a7 Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Thu, 9 Oct 2025 21:48:47 +0200 Subject: [PATCH 2/5] Fix tests --- src/kete/rust/time.rs | 15 ++++++++++----- src/tests/test_spice.py | 6 +++--- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/src/kete/rust/time.rs b/src/kete/rust/time.rs index 5677033..88adf47 100644 --- a/src/kete/rust/time.rs +++ b/src/kete/rust/time.rs @@ -68,7 +68,8 @@ impl PyTime { /// Construct a new time object, TDB default. #[new] #[pyo3(signature = (jd, scaling="tdb"))] - pub fn new(jd: f64, scaling: &str) -> PyResult { + pub fn new(jd: PyTime, scaling: &str) -> PyResult { + let jd = jd.jd(); Ok(match scaling.to_ascii_lowercase().as_str() { "tt" => PyTime(Time::::new(jd)), "tdb" => PyTime(Time::::new(jd)), @@ -235,15 +236,19 @@ impl PyTime { (self.0.jd + other.0.jd).into() } - fn __sub__(&self, other: PyTime) -> Self { - (self.0.jd - other.0.jd).into() + fn __sub__(&self, other: PyTime) -> f64 { + self.0.jd - other.0.jd } - fn __rsub__(&self, other: PyTime) -> Self { - (other.0.jd - self.0.jd).into() + fn __rsub__(&self, other: PyTime) -> f64 { + other.0.jd - self.0.jd } fn __repr__(&self) -> String { format!("Time({})", self.0.jd) } + + fn __eq__(&self, other: PyTime) -> bool { + self.0 == other.0 + } } diff --git a/src/tests/test_spice.py b/src/tests/test_spice.py index 7684097..3dd7c83 100644 --- a/src/tests/test_spice.py +++ b/src/tests/test_spice.py @@ -70,10 +70,10 @@ def test_loaded_info(self): info = spice.loaded_object_info("ceres") expected = SpkInfo( name="ceres", - jd_start=2415020.5, - jd_end=2488069.5, + jd_start=Time(2415020.5), + jd_end=Time(2488069.5), center=10, - frame=Frames.Equatorial, + frame=1, spk_type=21, ) assert info == [expected] From 7be4413e00fd803514ad4e05222cf05b437ff21a Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Thu, 9 Oct 2025 21:53:21 +0200 Subject: [PATCH 3/5] doc typo --- src/kete_core/src/time/mod.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/kete_core/src/time/mod.rs b/src/kete_core/src/time/mod.rs index 7d68a9c..d5c30e2 100644 --- a/src/kete_core/src/time/mod.rs +++ b/src/kete_core/src/time/mod.rs @@ -347,7 +347,7 @@ impl From for Duration { impl Sub> for Time { type Output = Duration; - /// Subtract two times, returning the difference in days of TDB. + /// Subtract two times, returning the duration in days of TDB. fn sub(self, other: Time) -> Self::Output { (self.tdb().jd - other.tdb().jd).into() } @@ -365,7 +365,7 @@ impl Add> for Time { impl Add for Time { type Output = Self; - /// Add a time and duration together + /// Add a time and duration together, assuming the duration is the same scaling fn add(self, other: f64) -> Self::Output { A::from_tdb(self.tdb().jd + other).into() } @@ -374,7 +374,7 @@ impl Add for Time { impl Sub> for Time { type Output = Self; - /// Subtract two times, returning the difference in days. + /// Subtract a duration from a time fn sub(self, other: Duration) -> Self::Output { A::from_tdb(self.tdb().jd - other.tdb().elapsed).into() } @@ -383,7 +383,7 @@ impl Sub> for Time { impl Sub for Time { type Output = Self; - /// Add a time and duration together + /// Subtract a duration from a time, assuming the duration is the same scaling fn sub(self, other: f64) -> Self::Output { A::from_tdb(self.tdb().jd - other).into() } From eb1ff92ffefef6db2398bf9f10a3817b13799be0 Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Fri, 10 Oct 2025 09:34:22 +0200 Subject: [PATCH 4/5] typo in doc string --- src/kete_core/src/spice/spk.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/kete_core/src/spice/spk.rs b/src/kete_core/src/spice/spk.rs index 5d97790..aba740a 100644 --- a/src/kete_core/src/spice/spk.rs +++ b/src/kete_core/src/spice/spk.rs @@ -14,7 +14,7 @@ //! let singleton = LOADED_SPK.try_read().unwrap(); //! //! // get the state of 399 (Earth) -//! let state = singleton.try_get_state::(399, 2451545.0); +//! let state = singleton.try_get_state::(399, 2451545.0.into()); //! ``` //! //! From 7c20c4fe8005c61c6dfd6ffd25fd55f5c14e2dc9 Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Fri, 10 Oct 2025 22:20:58 +0200 Subject: [PATCH 5/5] Time __lt__ --- src/kete/rust/time.rs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/kete/rust/time.rs b/src/kete/rust/time.rs index 88adf47..46b4ac2 100644 --- a/src/kete/rust/time.rs +++ b/src/kete/rust/time.rs @@ -251,4 +251,8 @@ impl PyTime { fn __eq__(&self, other: PyTime) -> bool { self.0 == other.0 } + + fn __lt__(&self, other: PyTime) -> bool { + self.0.jd < other.0.jd + } }