diff --git a/crates/orbit/apsis-sgp4/src/batch.rs b/crates/orbit/apsis-sgp4/src/batch.rs new file mode 100644 index 0000000..23ecd1f --- /dev/null +++ b/crates/orbit/apsis-sgp4/src/batch.rs @@ -0,0 +1,63 @@ +//! Batch propagation over many instants, optionally across several threads. +//! +//! A single satellite's track is the common case for a coverage run or an +//! access search: the caller wants the state at every minute of a day, for +//! every satellite in a shell. Each state is independent of the others, so +//! the loop is a candidate for several cores, and the work is large enough +//! that starting a thread pays for itself. +//! +//! The parallel entry point produces the same values as the sequential one, +//! in the same order. It differs only in how the work is distributed. + +use crate::{Satellite, State}; +use apsis_math::parallel::map_index; + +/// Propagates a satellite at a list of minutes from its epoch, in order. +/// +/// The sequential form, provided for comparison and for small lists where +/// starting a thread would cost more than the work. +pub fn propagate_minutes_at(satellite: &Satellite, minutes: &[f64]) -> Vec> { + minutes + .iter() + .map(|&value| satellite.propagate_minutes(value).ok()) + .collect() +} + +/// Propagates a satellite at a list of minutes from its epoch, across several +/// threads. +/// +/// The result matches [`propagate_minutes_at`] element for element and in the +/// same order. A minute at which the model reports a decay yields `None`, as +/// it does in the sequential form; the threads do not stop at the first such +/// minute, because a caller usually wants the whole track and will decide for +/// itself what a decay means. +/// +/// Propagation is deterministic and depends only on the minutes, so the +/// result does not depend on how the work was split. +pub fn propagate_minutes_at_parallel(satellite: &Satellite, minutes: &[f64]) -> Vec> { + map_index(minutes.len(), |index| { + satellite.propagate_minutes(minutes[index]).ok() + }) +} + +/// Propagates several satellites at one instant, in order. +/// +/// The per satellite work is independent, so a shell's state at a single +/// instant is a natural batch. The satellites are read only, so they are +/// shared across the threads by reference. +pub fn propagate_at(satellites: &[Satellite], tai_mjd: f64) -> Vec> { + satellites + .iter() + .map(|satellite| satellite.propagate_at(tai_mjd).ok()) + .collect() +} + +/// Propagates several satellites at one instant, across several threads. +/// +/// The result matches [`propagate_at`] element for element and in the same +/// order. +pub fn propagate_at_parallel(satellites: &[Satellite], tai_mjd: f64) -> Vec> { + map_index(satellites.len(), |index| { + satellites[index].propagate_at(tai_mjd).ok() + }) +} diff --git a/crates/orbit/apsis-sgp4/src/lib.rs b/crates/orbit/apsis-sgp4/src/lib.rs index a3bfcb8..3387aee 100644 --- a/crates/orbit/apsis-sgp4/src/lib.rs +++ b/crates/orbit/apsis-sgp4/src/lib.rs @@ -53,6 +53,7 @@ #![forbid(unsafe_code)] #![deny(missing_docs)] +mod batch; mod deep_space; mod error; mod init; @@ -61,6 +62,9 @@ mod propagate; mod satellite; mod state; +pub use batch::{ + propagate_at, propagate_at_parallel, propagate_minutes_at, propagate_minutes_at_parallel, +}; pub use error::{PropagationError, PropagationErrorKind}; pub use init::Elements; pub use model::GravityModel; diff --git a/crates/orbit/apsis-sgp4/tests/batch.rs b/crates/orbit/apsis-sgp4/tests/batch.rs new file mode 100644 index 0000000..9a79e89 --- /dev/null +++ b/crates/orbit/apsis-sgp4/tests/batch.rs @@ -0,0 +1,99 @@ +//! Tests for the batch propagation entry points. +//! +//! The property that matters is that the parallel forms produce exactly the +//! values the sequential forms do, in the same order. Every test compares +//! against the sequential answer rather than against a constant, so a +//! mismatch in either the propagation or the ordering is caught. + +use apsis_sgp4::{ + GravityModel, Satellite, TwoLineElement, propagate_at, propagate_at_parallel, + propagate_minutes_at, propagate_minutes_at_parallel, +}; + +const SAMPLE: &str = concat!( + "1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927\r\n", + "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537\r\n", +); + +fn satellite() -> Satellite { + let element = TwoLineElement::parse(SAMPLE).expect("the sample parses"); + Satellite::from_element(&element, GravityModel::WGS72).expect("the orbit propagates") +} + +#[test] +fn a_short_track_matches_the_sequential_form() { + let satellite = satellite(); + let minutes: Vec = (0..100).map(|index| index as f64 * 1.0).collect(); + let sequential = propagate_minutes_at(&satellite, &minutes); + let parallel = propagate_minutes_at_parallel(&satellite, &minutes); + assert_eq!(parallel, sequential); +} + +#[test] +fn a_long_track_matches_the_sequential_form() { + // A day at a one minute step, which is long enough that the parallel path + // takes the multi thread branch. + let satellite = satellite(); + let minutes: Vec = (0..1_441).map(|index| index as f64).collect(); + let sequential = propagate_minutes_at(&satellite, &minutes); + let parallel = propagate_minutes_at_parallel(&satellite, &minutes); + assert_eq!(parallel, sequential); + assert_eq!(parallel.len(), 1_441); +} + +#[test] +fn a_backwards_track_matches() { + // Negative minutes propagate backwards and must behave the same way in + // both forms. + let satellite = satellite(); + let minutes: Vec = (0..2_000).map(|index| -(index as f64) * 0.5).collect(); + assert_eq!( + propagate_minutes_at_parallel(&satellite, &minutes), + propagate_minutes_at(&satellite, &minutes) + ); +} + +#[test] +fn an_empty_track_is_empty() { + let satellite = satellite(); + assert!(propagate_minutes_at_parallel(&satellite, &[]).is_empty()); +} + +#[test] +fn a_decayed_instant_is_reported_in_both_forms() { + // Far enough ahead that the drag term brings the satellite down, both + // forms must agree that the instant is unusable rather than one of them + // stopping early. + let satellite = satellite(); + let minutes = vec![0.0, 1.0, 60.0 * 24.0 * 365.0 * 20.0]; + let sequential = propagate_minutes_at(&satellite, &minutes); + let parallel = propagate_minutes_at_parallel(&satellite, &minutes); + assert_eq!(parallel, sequential); +} + +#[test] +fn a_shell_at_one_instant_matches() { + let satellite = satellite(); + let shell = vec![satellite; 500]; + let instant = 54_723.5; + assert_eq!( + propagate_at_parallel(&shell, instant), + propagate_at(&shell, instant) + ); +} + +#[test] +fn a_shell_result_is_in_input_order() { + // The states of a repeated satellite are identical, so ordering is checked + // through the minutes form, where each value differs. + let satellite = satellite(); + let minutes: Vec = (0..3_000).map(|index| index as f64 * 0.25).collect(); + let parallel = propagate_minutes_at_parallel(&satellite, &minutes); + for (index, state) in parallel.iter().enumerate() { + let state = state.as_ref().expect("the orbit has not decayed here"); + let direct = satellite + .propagate_minutes(minutes[index]) + .expect("the orbit propagates"); + assert_eq!(state.position, direct.position, "at index {index}"); + } +}