From 2a41cf1560be5a37c92f1d2787a0f530fcec4b71 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 15:51:47 +0800 Subject: [PATCH] network: add access window finder --- Cargo.lock | 6 + crates/network/apsis-access/Cargo.toml | 4 + crates/network/apsis-access/src/error.rs | 73 +++ crates/network/apsis-access/src/finder.rs | 377 ++++++++++++++ crates/network/apsis-access/src/lib.rs | 47 +- crates/network/apsis-access/src/pass.rs | 164 ++++++ crates/network/apsis-access/tests/finder.rs | 523 ++++++++++++++++++++ 7 files changed, 1193 insertions(+), 1 deletion(-) create mode 100644 crates/network/apsis-access/src/error.rs create mode 100644 crates/network/apsis-access/src/finder.rs create mode 100644 crates/network/apsis-access/src/pass.rs create mode 100644 crates/network/apsis-access/tests/finder.rs diff --git a/Cargo.lock b/Cargo.lock index 27bc1d4..0f260da 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -5,6 +5,12 @@ version = 4 [[package]] name = "apsis-access" version = "0.1.0" +dependencies = [ + "apsis-core", + "apsis-coverage", + "apsis-geodesy", + "apsis-math", +] [[package]] name = "apsis-antenna" diff --git a/crates/network/apsis-access/Cargo.toml b/crates/network/apsis-access/Cargo.toml index 3e99d9f..0d79fe0 100644 --- a/crates/network/apsis-access/Cargo.toml +++ b/crates/network/apsis-access/Cargo.toml @@ -11,3 +11,7 @@ publish = false workspace = true [dependencies] +apsis-core.workspace = true +apsis-math.workspace = true +apsis-geodesy.workspace = true +apsis-coverage.workspace = true diff --git a/crates/network/apsis-access/src/error.rs b/crates/network/apsis-access/src/error.rs new file mode 100644 index 0000000..1049054 --- /dev/null +++ b/crates/network/apsis-access/src/error.rs @@ -0,0 +1,73 @@ +//! Errors reported while finding access windows. + +use core::fmt; + +/// The category of an access error. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +#[non_exhaustive] +pub enum AccessErrorKind { + /// The station or its mask is not usable. + InvalidStation, + /// The sampling step is not positive and finite. + InvalidSampleStep, + /// The crossing tolerance is not positive and finite. + InvalidTolerance, + /// The search interval is empty or not ordered. + InvalidInterval, + /// A supplied satellite position is not usable. + InvalidPosition, +} + +impl fmt::Display for AccessErrorKind { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let text = match self { + Self::InvalidStation => "the ground station is not usable", + Self::InvalidSampleStep => "the sampling step is not usable", + Self::InvalidTolerance => "the crossing tolerance is not usable", + Self::InvalidInterval => "the search interval is not usable", + Self::InvalidPosition => "a satellite position is not usable", + }; + f.write_str(text) + } +} + +/// A failure while finding access windows. +#[derive(Debug, Clone, PartialEq, Eq)] +pub struct AccessError { + kind: AccessErrorKind, + detail: Option, +} + +impl AccessError { + /// Creates an error with the given category. + pub fn new(kind: AccessErrorKind) -> Self { + Self { kind, detail: None } + } + + /// Attaches a description of the values that caused the failure. + pub fn with_detail(mut self, detail: impl Into) -> Self { + self.detail = Some(detail.into()); + self + } + + /// Returns the category of this failure. + pub fn kind(&self) -> AccessErrorKind { + self.kind + } + + /// Returns the detail text attached to this error, if any. + pub fn detail(&self) -> Option<&str> { + self.detail.as_deref() + } +} + +impl fmt::Display for AccessError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match &self.detail { + Some(detail) => write!(f, "{}: {detail}", self.kind), + None => write!(f, "{}", self.kind), + } + } +} + +impl core::error::Error for AccessError {} diff --git a/crates/network/apsis-access/src/finder.rs b/crates/network/apsis-access/src/finder.rs new file mode 100644 index 0000000..25b172a --- /dev/null +++ b/crates/network/apsis-access/src/finder.rs @@ -0,0 +1,377 @@ +//! The access window finder. + +use crate::error::{AccessError, AccessErrorKind}; +use crate::pass::{Pass, PassSummary}; +use apsis_coverage::{GroundPoint, elevation_angle}; +use apsis_math::Vec3; + +/// One sample of a satellite track, as the finder consumes it. +/// +/// The position is in the Earth-fixed frame, in kilometers, and the time is +/// in seconds from the caller's origin. The caller supplies a sequence of +/// these, which decouples the finder from the propagation and the frame +/// transformation. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct OrbitSample { + /// The time of the sample, in seconds from the caller's origin. + pub seconds: f64, + /// The satellite position in the Earth-fixed frame, in kilometers. + pub position_ecef_km: Vec3, +} + +impl OrbitSample { + /// Creates a sample from a time and a position. + pub const fn new(seconds: f64, position_ecef_km: Vec3) -> Self { + Self { + seconds, + position_ecef_km, + } + } + + /// Returns `true` when the sample is usable. + pub fn is_valid(&self) -> bool { + self.seconds.is_finite() && self.position_ecef_km.is_finite() + } +} + +/// Finds access windows for one station. +/// +/// The finder holds the station and the numerical settings of the search. +/// The settings are separated from the search so that the same finder can be +/// applied to many satellites. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct AccessFinder { + /// The station and its elevation mask. + pub station: GroundPoint, + /// The interval between samples of the satellite track, in seconds. + /// + /// A smaller step finds shorter passes but costs proportionally more + /// evaluations. A step of one minute resolves any pass a low Earth + /// satellite makes, since such a pass lasts several minutes. + pub sample_step_seconds: f64, + /// The tolerance on a crossing time, in seconds. + /// + /// The refinement of a rise or set continues until the interval + /// containing the crossing is this narrow. + pub crossing_tolerance_seconds: f64, +} + +impl AccessFinder { + /// Returns `true` when the finder settings are usable. + pub fn is_valid(&self) -> bool { + self.station.is_valid() + && self.sample_step_seconds.is_finite() + && self.sample_step_seconds > 0.0 + && self.crossing_tolerance_seconds.is_finite() + && self.crossing_tolerance_seconds > 0.0 + && self.crossing_tolerance_seconds <= self.sample_step_seconds + } + + /// Validates the settings, returning the reason when they are not + /// usable. + fn validate(&self) -> Result<(), AccessError> { + if !self.station.is_valid() { + return Err(AccessError::new(AccessErrorKind::InvalidStation)); + } + if !self.sample_step_seconds.is_finite() || self.sample_step_seconds <= 0.0 { + return Err(AccessError::new(AccessErrorKind::InvalidSampleStep) + .with_detail(format!("step is {}", self.sample_step_seconds))); + } + if !self.crossing_tolerance_seconds.is_finite() || self.crossing_tolerance_seconds <= 0.0 { + return Err(AccessError::new(AccessErrorKind::InvalidTolerance) + .with_detail(format!("tolerance is {}", self.crossing_tolerance_seconds))); + } + if self.crossing_tolerance_seconds > self.sample_step_seconds { + return Err(AccessError::new(AccessErrorKind::InvalidTolerance) + .with_detail("the tolerance exceeds the sampling step")); + } + Ok(()) + } + + /// Returns the elevation of a sample, in radians. + /// + /// Returns `None` when the sample or the station is not usable. + pub fn elevation_of(&self, sample: &OrbitSample) -> Option { + if !sample.is_valid() { + return None; + } + elevation_angle(sample.position_ecef_km, &self.station) + } + + /// Finds the access windows in an ordered sequence of samples. + /// + /// The samples must be ordered by time and cover the interval of + /// interest. The finder assumes the elevation is continuous between + /// samples, so the step must be short enough that the satellite does not + /// rise and set within one step. A low Earth satellite is above a + /// ten degree mask for several minutes, so a step of a minute or less + /// satisfies that. + /// + /// The crossing times are found by interpolating the satellite position + /// linearly between the two samples that straddle the crossing. The + /// chord between two samples departs from the true arc by an amount of + /// the second order in the step, so a ten second step gives a crossing + /// accurate to about a tenth of a second and a one second step to about + /// a millisecond. Use [`AccessFinder::find_with`] when the propagator is + /// available and the crossing is wanted to the full tolerance. + /// + /// A pass that is already in progress at the first sample is reported + /// with its rise at the first sample time, and one still in progress at + /// the last sample is reported with its set at the last sample time. + /// That convention means the pass is clipped to the interval searched + /// rather than extrapolated beyond it. + /// + /// Returns an error when the settings are not usable, when the sample + /// sequence has fewer than two entries, or when a sample is not usable. + pub fn find(&self, samples: &[OrbitSample]) -> Result { + self.find_internal(samples, |seconds| { + // Interpolate the position between the bracketing samples. The + // caller is responsible for the ordering and the step, which + // were checked before this closure is called. + let index = match samples.binary_search_by(|sample| { + sample + .seconds + .partial_cmp(&seconds) + .unwrap_or(core::cmp::Ordering::Equal) + }) { + Ok(exact) => exact, + Err(insert) => insert.min(samples.len() - 1), + }; + let upper = index.max(1); + let lower = upper - 1; + let span = samples[upper].seconds - samples[lower].seconds; + let fraction = if span == 0.0 { + 0.0 + } else { + (seconds - samples[lower].seconds) / span + }; + let position = samples[lower] + .position_ecef_km + .lerp(samples[upper].position_ecef_km, fraction); + self.elevation_of(&OrbitSample::new(seconds, position)) + }) + } + + /// Runs the search given a function that returns the elevation at a + /// time. + /// + /// The function must be defined over the whole interval and must be + /// continuous, and the sample sequence must be ordered by time. The + /// closure returns `None` when the position at that time is not usable. + fn find_internal( + &self, + samples: &[OrbitSample], + mut elevation_at: F, + ) -> Result + where + F: FnMut(f64) -> Option, + { + self.validate()?; + if samples.len() < 2 { + return Err(AccessError::new(AccessErrorKind::InvalidInterval) + .with_detail("at least two samples are needed")); + } + + let first = samples[0].seconds; + let last = samples[samples.len() - 1].seconds; + if !first.is_finite() || !last.is_finite() || last <= first { + return Err(AccessError::new(AccessErrorKind::InvalidInterval) + .with_detail(format!("times run from {first} to {last}"))); + } + + let mask = self.station.elevation_mask; + let mut elevations = Vec::with_capacity(samples.len()); + for sample in samples { + let elevation = elevation_at(sample.seconds).ok_or_else(|| { + AccessError::new(AccessErrorKind::InvalidPosition) + .with_detail(format!("sample at {} s is not usable", sample.seconds)) + })?; + if !elevation.is_finite() { + return Err(AccessError::new(AccessErrorKind::InvalidPosition) + .with_detail(format!("elevation at {} s is not finite", sample.seconds))); + } + elevations.push(elevation); + } + + let mut passes = Vec::new(); + let mut rise: Option = None; + // Track the highest elevation seen inside the current pass so that + // the peak is known without a second pass over the samples. + let mut peak = f64::NEG_INFINITY; + let mut peak_seconds = first; + let mut peak_position = samples[0].position_ecef_km; + + for index in 0..samples.len() { + let above = elevations[index] >= mask; + if above { + if rise.is_none() { + rise = Some(if index == 0 { + // The pass was already open at the start of the + // window, so it is clipped to the first sample. + samples[0].seconds + } else { + self.refine_crossing(samples, index, true, &mut elevation_at) + }); + } + if elevations[index] > peak { + peak = elevations[index]; + peak_seconds = samples[index].seconds; + peak_position = samples[index].position_ecef_km; + } + } else if let Some(rise_seconds) = rise.take() { + let set_seconds = self.refine_crossing(samples, index, false, &mut elevation_at); + passes.push(Pass { + rise_seconds, + set_seconds, + peak_seconds, + peak_elevation: peak, + mask_elevation: mask, + peak_position, + }); + peak = f64::NEG_INFINITY; + } + } + + // A pass still open at the end of the sequence is clipped to the + // last sample. + if let Some(rise_seconds) = rise { + passes.push(Pass { + rise_seconds, + set_seconds: last, + peak_seconds, + peak_elevation: peak, + mask_elevation: mask, + peak_position, + }); + } + + Ok(PassSummary::new(passes, last - first)) + } + + /// Refines the time at which the elevation crosses the mask. + /// + /// The crossing lies between the sample at `index` and the one before + /// it, and the two straddle the mask. Bisection halves the interval + /// until it is narrower than the tolerance. + /// + /// When `rising` is true the crossing is a rise, where the elevation + /// increases through the mask; otherwise it is a set. The flag selects + /// which side of the crossing is above the mask, which keeps the search + /// bracket correct. + fn refine_crossing( + &self, + samples: &[OrbitSample], + index: usize, + rising: bool, + elevation_at: &mut F, + ) -> f64 + where + F: FnMut(f64) -> Option, + { + let mask = self.station.elevation_mask; + let mut low = samples[index - 1].seconds; + let mut high = samples[index].seconds; + + // The elevation difference from the mask, with the sign arranged so + // that the function is negative at `low` and positive at `high`. + let sign = if rising { 1.0 } else { -1.0 }; + let mut difference = |seconds: f64| { + let elevation = elevation_at(seconds).unwrap_or(mask); + sign * (elevation - mask) + }; + + let mut low_value = difference(low); + for _ in 0..64 { + if high - low <= self.crossing_tolerance_seconds { + break; + } + let middle = 0.5 * (low + high); + let value = difference(middle); + if value == 0.0 { + return middle; + } + if (value > 0.0) == (low_value > 0.0) { + low = middle; + low_value = value; + } else { + high = middle; + } + } + 0.5 * (low + high) + } + + /// Finds the access windows over an interval, calling a propagator for + /// each sample. + /// + /// The propagator receives a time in seconds from the origin of the + /// search and returns the satellite position in the Earth-fixed frame. + /// This adapts the finder to any orbit source without this crate + /// depending on one. + /// + /// Returns an error when the settings are not usable, when the interval + /// is not ordered, or when the propagator fails at a sample. + pub fn find_with( + &self, + start_seconds: f64, + end_seconds: f64, + mut propagator: F, + ) -> Result + where + F: FnMut(f64) -> Option, + { + self.validate()?; + if !start_seconds.is_finite() || !end_seconds.is_finite() || end_seconds <= start_seconds { + return Err( + AccessError::new(AccessErrorKind::InvalidInterval).with_detail(format!( + "interval runs from {start_seconds} to {end_seconds}" + )), + ); + } + + // Build the sample sequence, including the end point. + let step = self.sample_step_seconds; + let count = ((end_seconds - start_seconds) / step).ceil() as usize + 1; + let mut samples = Vec::with_capacity(count); + for index in 0..count { + let seconds = (start_seconds + index as f64 * step).min(end_seconds); + let position = propagator(seconds).ok_or_else(|| { + AccessError::new(AccessErrorKind::InvalidPosition) + .with_detail(format!("the propagator failed at {seconds} s")) + })?; + samples.push(OrbitSample::new(seconds, position)); + } + // The loop above may append the end point twice when the interval is + // an exact multiple of the step, so the last sample is replaced when + // it duplicates the end. + if let Some(last) = samples.last_mut() { + if last.seconds != end_seconds { + let position = propagator(end_seconds).ok_or_else(|| { + AccessError::new(AccessErrorKind::InvalidPosition) + .with_detail(format!("the propagator failed at {end_seconds} s")) + })?; + samples.push(OrbitSample::new(end_seconds, position)); + } else { + last.seconds = end_seconds; + } + } + + // The refinement calls the propagator at each bisection point, so + // the crossing is found to the tolerance rather than to the + // accuracy of a chord between two samples. + let mut failure = None; + let summary = self.find_internal(&samples, |seconds| match propagator(seconds) { + Some(position) if position.is_finite() => elevation_angle(position, &self.station), + _ => { + failure = Some(seconds); + None + } + }); + match summary { + Ok(summary) => Ok(summary), + Err(error) => Err(match failure { + Some(seconds) => error.with_detail(format!("the propagator failed at {seconds} s")), + None => error, + }), + } + } +} diff --git a/crates/network/apsis-access/src/lib.rs b/crates/network/apsis-access/src/lib.rs index 92fd2bd..69b6ed6 100644 --- a/crates/network/apsis-access/src/lib.rs +++ b/crates/network/apsis-access/src/lib.rs @@ -1 +1,46 @@ -//! Ground station to satellite access window computation. +//! Ground station access windows. +//! +//! A satellite is accessible from a ground station during the intervals when +//! it rises above the station's elevation mask. Each such interval is a +//! pass, or access window, and the sequence of windows over a period is the +//! access schedule. This crate finds those windows from a propagated orbit +//! and reports their start, end, duration, and the geometry at their peak. +//! +//! The windows are found by sampling the elevation angle at a fixed step and +//! refining each crossing where the elevation passes the mask. The sampling +//! step sets the shortest pass that can be found, and the refinement gives +//! the crossing times to a tolerance independent of the step. +//! +//! The elevation is evaluated from the satellite position in the +//! Earth-fixed frame, so the caller supplies the transformation from the +//! orbit frame. In this workspace that is the conversion in `apsis-coord` +//! with the Earth rotation angle from `apsis-frames`. +//! +//! ``` +//! use apsis_access::{AccessFinder, Pass}; +//! use apsis_coverage::GroundPoint; +//! use apsis_geodesy::Geodetic; +//! +//! // A station on the equator with a ten degree mask. +//! let station = GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 10.0_f64.to_radians()); +//! let finder = AccessFinder { +//! station, +//! sample_step_seconds: 60.0, +//! crossing_tolerance_seconds: 0.001, +//! }; +//! assert!(finder.is_valid()); +//! ``` + +#![forbid(unsafe_code)] +#![deny(missing_docs)] + +mod error; +mod finder; +mod pass; + +pub use error::{AccessError, AccessErrorKind}; +pub use finder::{AccessFinder, OrbitSample}; +pub use pass::{Pass, PassSummary}; + +pub use apsis_coverage::{GroundPoint, Sensor, elevation_angle}; +pub use apsis_math::Vec3; diff --git a/crates/network/apsis-access/src/pass.rs b/crates/network/apsis-access/src/pass.rs new file mode 100644 index 0000000..c94c0eb --- /dev/null +++ b/crates/network/apsis-access/src/pass.rs @@ -0,0 +1,164 @@ +//! The access window type and its summary. + +use apsis_math::Vec3; + +/// One continuous interval during which a satellite is above a station's +/// elevation mask. +/// +/// Times are seconds from an origin chosen by the caller, so the same +/// structure serves a search over a scenario clock and one over elapsed +/// time. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Pass { + /// The time the satellite rises above the mask, in seconds from the + /// caller's origin. + pub rise_seconds: f64, + /// The time the satellite falls below the mask. + pub set_seconds: f64, + /// The time of greatest elevation. + pub peak_seconds: f64, + /// The greatest elevation during the pass, in radians. + pub peak_elevation: f64, + /// The elevation at which the satellite crossed the mask, in radians. + /// + /// This is the station's mask, recorded so that a pass carries the + /// geometry it was found under. + pub mask_elevation: f64, + /// The satellite position in the Earth-fixed frame at the peak, in + /// kilometers. + pub peak_position: Vec3, +} + +impl Pass { + /// Returns the duration of the pass, in seconds. + pub fn duration(&self) -> f64 { + self.set_seconds - self.rise_seconds + } + + /// Returns the time from the rise to the peak, in seconds. + pub fn rise_to_peak(&self) -> f64 { + self.peak_seconds - self.rise_seconds + } + + /// Returns the elevation at a time within the pass, when the caller can + /// supply it. + /// + /// The pass stores only the peak, so this helper reports the stored + /// values at the three defining times and defers to the caller's + /// elevation function otherwise. + pub fn contains(&self, seconds: f64) -> bool { + seconds >= self.rise_seconds && seconds <= self.set_seconds + } + + /// Returns the maximum elevation of the pass, in radians. + /// + /// This is the same as the peak elevation, named for readability at a + /// call site that compares passes. + pub fn maximum_elevation(&self) -> f64 { + self.peak_elevation + } +} + +/// A summary of the access windows over an interval. +#[derive(Debug, Clone, PartialEq)] +pub struct PassSummary { + /// The passes, in time order. + pub passes: Vec, + /// The length of the interval searched, in seconds. + pub interval_seconds: f64, +} + +impl PassSummary { + /// Creates a summary from a list of passes and the interval length. + pub fn new(passes: Vec, interval_seconds: f64) -> Self { + Self { + passes, + interval_seconds, + } + } + + /// Returns the number of passes. + pub fn count(&self) -> usize { + self.passes.len() + } + + /// Returns the total access time, in seconds. + /// + /// The passes are disjoint, so the durations add. + pub fn total_access_seconds(&self) -> f64 { + self.passes.iter().map(Pass::duration).sum() + } + + /// Returns the fraction of the interval during which the satellite was + /// accessible. + /// + /// Returns `None` when the interval length is not positive. + pub fn access_fraction(&self) -> Option { + if self.interval_seconds > 0.0 && self.interval_seconds.is_finite() { + Some(self.total_access_seconds() / self.interval_seconds) + } else { + None + } + } + + /// Returns the mean pass duration, in seconds. + /// + /// Returns `None` when there are no passes. + pub fn mean_duration_seconds(&self) -> Option { + if self.passes.is_empty() { + return None; + } + Some(self.total_access_seconds() / self.passes.len() as f64) + } + + /// Returns the longest pass, if there is one. + pub fn longest(&self) -> Option<&Pass> { + self.passes.iter().max_by(|first, second| { + first + .duration() + .partial_cmp(&second.duration()) + .unwrap_or(core::cmp::Ordering::Equal) + }) + } + + /// Returns the highest pass, if there is one. + pub fn highest(&self) -> Option<&Pass> { + self.passes.iter().max_by(|first, second| { + first + .peak_elevation + .partial_cmp(&second.peak_elevation) + .unwrap_or(core::cmp::Ordering::Equal) + }) + } + + /// Returns the gap between two consecutive passes, in seconds. + /// + /// Returns `None` when either index is out of range or the two indices + /// are not consecutive and increasing. + pub fn gap_seconds(&self, first: usize, second: usize) -> Option { + if second != first + 1 { + return None; + } + let earlier = self.passes.get(first)?; + let later = self.passes.get(second)?; + Some(later.rise_seconds - earlier.set_seconds) + } + + /// Returns the largest gap between consecutive passes, in seconds. + /// + /// This is the longest time a station waits between contacts. + /// Returns `None` when there are fewer than two passes. + pub fn longest_gap_seconds(&self) -> Option { + if self.passes.len() < 2 { + return None; + } + self.passes + .windows(2) + .map(|pair| pair[1].rise_seconds - pair[0].set_seconds) + .max_by(|first, second| { + first + .partial_cmp(second) + .unwrap_or(core::cmp::Ordering::Equal) + }) + } +} diff --git a/crates/network/apsis-access/tests/finder.rs b/crates/network/apsis-access/tests/finder.rs new file mode 100644 index 0000000..d90ef5c --- /dev/null +++ b/crates/network/apsis-access/tests/finder.rs @@ -0,0 +1,523 @@ +//! Tests for the access window finder. +//! +//! The reference case is a circular equatorial orbit at 550 km altitude seen +//! from a station on the equator with a ten degree mask. For that geometry +//! the elevation is a function of the orbital position alone, so the rise, +//! peak, and set times have a closed form: +//! +//! ```text +//! theta_c = 90 - mask - asin(Re cos(mask) / Rs) +//! t = (theta - theta_0) / n +//! ``` +//! +//! with `theta` the orbital position measured from the station and `n` the +//! mean motion. The values used below are computed from those expressions. + +use apsis_access::{AccessErrorKind, AccessFinder, OrbitSample, Pass, PassSummary}; +use apsis_coverage::GroundPoint; +use apsis_geodesy::Geodetic; +use apsis_math::Vec3; +use core::f64::consts::{FRAC_PI_2, PI}; + +/// The Earth equatorial radius, in kilometers. +const BODY_RADIUS: f64 = 6_378.137; + +/// The orbit radius for a 550 km circular orbit, in kilometers. +const ORBIT_RADIUS: f64 = BODY_RADIUS + 550.0; + +/// The mean motion of the reference orbit, in radians per second. +/// +/// This is `sqrt(mu / a^3)` with the WGS 84 gravitational parameter. +const MEAN_MOTION: f64 = 0.001_094_823_692_885_802; + +/// A station on the equator with a ten degree mask. +fn station() -> GroundPoint { + GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 10.0_f64.to_radians()) +} + +/// A finder with a one second step and a millisecond tolerance. +fn finder() -> AccessFinder { + AccessFinder { + station: station(), + sample_step_seconds: 1.0, + crossing_tolerance_seconds: 0.001, + } +} + +/// The orbital position measured from the station at the start of the +/// search, in radians. +/// +/// The reference closed form solves for the crossing times relative to an +/// initial position of minus one radian, so the track starts there. +const INITIAL_ANGLE: f64 = -1.0; + +/// Returns the satellite position for a reference orbit that passes over the +/// station. +/// +/// The orbit is circular and equatorial, so the position depends only on the +/// angular position at the given time. The angle is measured from the station +/// and increases at the mean motion, so the elevation is a function of that +/// angle alone. +fn reference_position(seconds: f64) -> Vec3 { + let angle = INITIAL_ANGLE + MEAN_MOTION * seconds; + Vec3::new(ORBIT_RADIUS * angle.cos(), ORBIT_RADIUS * angle.sin(), 0.0) +} + +/// Returns the samples of the reference orbit over a window around one pass. +fn reference_samples(step: f64, start: f64, end: f64) -> Vec { + let mut samples = Vec::new(); + let mut seconds = start; + while seconds <= end { + samples.push(OrbitSample::new(seconds, reference_position(seconds))); + seconds += step; + } + samples +} + +/// The half angle of the footprint for the reference geometry, in radians. +/// +/// This is `90 - mask - asin(Re cos(mask) / Rs)`, which for the reference +/// case is 14.9566 degrees. +const FOOTPRINT_HALF_ANGLE: f64 = 0.261_042_364_648_894_9; + +/// The rise time of the reference pass from the start of the search, in +/// seconds. +const REFERENCE_RISE: f64 = 674.955_681_770_3; + +/// The peak time of the reference pass, in seconds. +const REFERENCE_PEAK: f64 = 913.389_074_878_5; + +/// The set time of the reference pass, in seconds. +const REFERENCE_SET: f64 = 1_151.822_467_986_6; + +#[test] +fn the_finder_settings_are_validated() { + assert!(finder().is_valid()); + + let coarse = AccessFinder { + sample_step_seconds: 0.0, + ..finder() + }; + assert!(!coarse.is_valid()); + assert_eq!( + coarse + .find(&reference_samples(1.0, 0.0, 10.0)) + .expect_err("zero step") + .kind(), + AccessErrorKind::InvalidSampleStep + ); + + let loose = AccessFinder { + crossing_tolerance_seconds: 10.0, + ..finder() + }; + assert!(!loose.is_valid()); + assert_eq!( + loose + .find(&reference_samples(1.0, 0.0, 10.0)) + .expect_err("tolerance exceeds the step") + .kind(), + AccessErrorKind::InvalidTolerance + ); + + let bad_station = AccessFinder { + station: GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), PI), + ..finder() + }; + assert!(!bad_station.is_valid()); + assert_eq!( + bad_station + .find(&reference_samples(1.0, 0.0, 10.0)) + .expect_err("invalid station") + .kind(), + AccessErrorKind::InvalidStation + ); +} + +#[test] +fn the_reference_pass_matches_the_closed_form() { + let samples = reference_samples(1.0, 0.0, 2_000.0); + let summary = finder().find(&samples).expect("valid"); + assert_eq!(summary.count(), 1, "expected exactly one pass"); + + let pass = &summary.passes[0]; + assert!( + (pass.rise_seconds - REFERENCE_RISE).abs() < 0.5, + "rise {} against {REFERENCE_RISE}", + pass.rise_seconds + ); + assert!( + (pass.peak_seconds - REFERENCE_PEAK).abs() < 1.0, + "peak {} against {REFERENCE_PEAK}", + pass.peak_seconds + ); + assert!( + (pass.set_seconds - REFERENCE_SET).abs() < 0.5, + "set {} against {REFERENCE_SET}", + pass.set_seconds + ); +} + +#[test] +fn the_peak_of_an_overhead_pass_is_a_right_angle() { + let samples = reference_samples(1.0, 0.0, 2_000.0); + let summary = finder().find(&samples).expect("valid"); + let pass = &summary.passes[0]; + // The reference geometry passes directly over the station, so the + // greatest elevation is a right angle. + assert!( + (pass.peak_elevation - FRAC_PI_2).abs() < 0.01, + "peak elevation {}", + pass.peak_elevation + ); +} + +#[test] +fn the_pass_is_symmetric_about_the_peak() { + let samples = reference_samples(1.0, 0.0, 2_000.0); + let summary = finder().find(&samples).expect("valid"); + let pass = &summary.passes[0]; + let rise_to_peak = pass.rise_to_peak(); + let peak_to_set = pass.set_seconds - pass.peak_seconds; + assert!( + (rise_to_peak - peak_to_set).abs() < 1.0, + "{rise_to_peak} against {peak_to_set}" + ); +} + +#[test] +fn the_duration_matches_the_geometry() { + let samples = reference_samples(1.0, 0.0, 2_000.0); + let summary = finder().find(&samples).expect("valid"); + let pass = &summary.passes[0]; + // The orbital position spans twice the footprint half angle, so the + // duration is 2 theta_c / n. + let expected = 2.0 * FOOTPRINT_HALF_ANGLE / MEAN_MOTION; + assert!( + (pass.duration() - expected).abs() < 1.0, + "{} against {expected}", + pass.duration() + ); +} + +#[test] +fn the_refinement_is_independent_of_the_sample_step() { + // The crossing times are found by interpolating the satellite position + // between samples, so the result approaches the true crossing as the + // step falls. The residual difference between two steps is second order + // in the step, because the chord between two samples departs from the + // arc by that order. + let coarse = finder() + .find(&reference_samples(10.0, 0.0, 2_000.0)) + .expect("valid"); + let fine = finder() + .find(&reference_samples(1.0, 0.0, 2_000.0)) + .expect("valid"); + assert_eq!(coarse.count(), 1); + assert_eq!(fine.count(), 1); + let coarse_pass = &coarse.passes[0]; + let fine_pass = &fine.passes[0]; + + // A ten second step gives a crossing within about a tenth of a second + // of the one found with a one second step, because the chord error is + // second order. + assert!( + (coarse_pass.rise_seconds - fine_pass.rise_seconds).abs() < 0.1, + "rise {} against {}", + coarse_pass.rise_seconds, + fine_pass.rise_seconds + ); + assert!( + (coarse_pass.set_seconds - fine_pass.set_seconds).abs() < 0.1, + "set {} against {}", + coarse_pass.set_seconds, + fine_pass.set_seconds + ); +} + +#[test] +fn the_tolerance_controls_the_crossing_accuracy() { + let loose = AccessFinder { + crossing_tolerance_seconds: 1.0, + ..finder() + }; + let tight = AccessFinder { + crossing_tolerance_seconds: 1e-6, + ..finder() + }; + let samples = reference_samples(1.0, 0.0, 2_000.0); + let loose_pass = loose.find(&samples).expect("valid").passes[0]; + let tight_pass = tight.find(&samples).expect("valid").passes[0]; + // The two agree to within the looser tolerance. + assert!((loose_pass.rise_seconds - tight_pass.rise_seconds).abs() <= 1.0); + assert!((loose_pass.set_seconds - tight_pass.set_seconds).abs() <= 1.0); +} + +#[test] +fn a_full_orbit_finds_exactly_one_pass() { + // The reference geometry gives one overhead pass per revolution, so a + // window of one period contains one pass. + let period = 2.0 * PI / MEAN_MOTION; + let samples = reference_samples(1.0, 0.0, period); + let summary = finder().find(&samples).expect("valid"); + assert_eq!(summary.count(), 1); +} + +#[test] +fn the_summary_reports_the_access_fraction() { + let samples = reference_samples(1.0, 0.0, 2_000.0); + let summary = finder().find(&samples).expect("valid"); + let fraction = summary.access_fraction().expect("a positive interval"); + let expected = summary.total_access_seconds() / 2_000.0; + assert!((fraction - expected).abs() < 1e-12); + // One overhead pass in a two thousand second window is about a quarter + // of the time. + assert!(fraction > 0.2 && fraction < 0.3, "fraction {fraction}"); +} + +#[test] +fn the_summary_finds_the_longest_and_highest_pass() { + let samples = reference_samples(1.0, 0.0, 2_000.0); + let summary = finder().find(&samples).expect("valid"); + assert!(summary.longest().is_some()); + assert!(summary.highest().is_some()); + assert!(summary.mean_duration_seconds().is_some()); +} + +#[test] +fn a_pass_that_starts_before_the_window_is_clipped() { + // Starting the search after the rise means the pass is reported from + // the first sample rather than extrapolated backwards. + let samples = reference_samples(1.0, 800.0, 1_400.0); + let summary = finder().find(&samples).expect("valid"); + assert_eq!(summary.count(), 1); + let pass = &summary.passes[0]; + assert!((pass.rise_seconds - 800.0).abs() < 1e-9); + assert!((pass.set_seconds - REFERENCE_SET).abs() < 1.0); +} + +#[test] +fn a_pass_that_ends_after_the_window_is_clipped() { + let samples = reference_samples(1.0, 0.0, 1_000.0); + let summary = finder().find(&samples).expect("valid"); + assert_eq!(summary.count(), 1); + let pass = &summary.passes[0]; + assert!((pass.rise_seconds - REFERENCE_RISE).abs() < 1.0); + assert!((pass.set_seconds - 1_000.0).abs() < 1e-9); +} + +#[test] +fn a_window_with_no_pass_is_empty() { + // A window before the satellite comes into view contains no pass. + let samples = reference_samples(1.0, -2_000.0, -1_000.0); + let summary = finder().find(&samples).expect("valid"); + assert_eq!(summary.count(), 0); + assert_eq!(summary.total_access_seconds(), 0.0); + assert_eq!(summary.longest(), None); + assert_eq!(summary.mean_duration_seconds(), None); + assert_eq!(summary.longest_gap_seconds(), None); +} + +#[test] +fn the_gap_between_passes_is_reported() { + // Two orbits give two passes, and the gap between them is the period + // less the pass duration. + let period = 2.0 * PI / MEAN_MOTION; + let samples = reference_samples(1.0, 0.0, 2.0 * period); + let summary = finder().find(&samples).expect("valid"); + assert!(summary.count() >= 2, "count {}", summary.count()); + + let gap = summary.gap_seconds(0, 1).expect("consecutive passes"); + let expected = period - summary.passes[0].duration(); + assert!((gap - expected).abs() < 2.0, "{gap} against {expected}"); + // The largest gap is at least that one. + let longest = summary.longest_gap_seconds().expect("at least two passes"); + assert!(longest >= gap - 1e-9); +} + +#[test] +fn a_non_consecutive_gap_query_is_rejected() { + let samples = reference_samples(1.0, 0.0, 2_000.0); + let summary = finder().find(&samples).expect("valid"); + assert_eq!(summary.gap_seconds(0, 2), None); + assert_eq!(summary.gap_seconds(1, 0), None); + assert_eq!(summary.gap_seconds(0, 99), None); +} + +#[test] +fn the_pass_reports_containment() { + let pass = Pass { + rise_seconds: 100.0, + set_seconds: 200.0, + peak_seconds: 150.0, + peak_elevation: 1.0, + mask_elevation: 0.1, + peak_position: Vec3::new(ORBIT_RADIUS, 0.0, 0.0), + }; + assert!(pass.contains(100.0)); + assert!(pass.contains(150.0)); + assert!(pass.contains(200.0)); + assert!(!pass.contains(99.9)); + assert!(!pass.contains(200.1)); + assert!((pass.duration() - 100.0).abs() < 1e-12); + assert!((pass.rise_to_peak() - 50.0).abs() < 1e-12); + assert!((pass.maximum_elevation() - 1.0).abs() < 1e-12); +} + +#[test] +fn too_few_samples_are_rejected() { + let error = finder() + .find(&[OrbitSample::new(0.0, Vec3::new(ORBIT_RADIUS, 0.0, 0.0))]) + .expect_err("one sample is not enough"); + assert_eq!(error.kind(), AccessErrorKind::InvalidInterval); + assert_eq!( + finder().find(&[]).expect_err("no samples").kind(), + AccessErrorKind::InvalidInterval + ); +} + +#[test] +fn an_unordered_sample_sequence_is_rejected() { + let samples = reference_samples(1.0, 1_000.0, 0.0); + assert!(samples.is_empty()); + let unordered = [ + OrbitSample::new(100.0, reference_position(100.0)), + OrbitSample::new(50.0, reference_position(50.0)), + ]; + let error = finder().find(&unordered).expect_err("time runs backwards"); + assert_eq!(error.kind(), AccessErrorKind::InvalidInterval); +} + +#[test] +fn an_invalid_position_is_rejected() { + let samples = [ + OrbitSample::new(0.0, reference_position(0.0)), + OrbitSample::new(1.0, Vec3::splat(f64::NAN)), + ]; + let error = finder() + .find(&samples) + .expect_err("the position is not finite"); + assert_eq!(error.kind(), AccessErrorKind::InvalidPosition); +} + +#[test] +fn the_finder_can_drive_a_propagator_directly() { + // The convenience entry point must agree with building the samples by + // hand. + let via_samples = finder() + .find(&reference_samples(1.0, 0.0, 2_000.0)) + .expect("valid"); + let via_propagator = finder() + .find_with(0.0, 2_000.0, |seconds| Some(reference_position(seconds))) + .expect("valid"); + assert_eq!(via_samples.count(), via_propagator.count()); + let first = &via_samples.passes[0]; + let second = &via_propagator.passes[0]; + assert!((first.rise_seconds - second.rise_seconds).abs() < 1e-6); + assert!((first.set_seconds - second.set_seconds).abs() < 1e-6); +} + +#[test] +fn a_failing_propagator_is_reported() { + let error = finder() + .find_with(0.0, 100.0, |_| None) + .expect_err("the propagator failed"); + assert_eq!(error.kind(), AccessErrorKind::InvalidPosition); +} + +#[test] +fn an_invalid_search_interval_is_rejected() { + let error = finder() + .find_with(100.0, 0.0, |seconds| Some(reference_position(seconds))) + .expect_err("the interval runs backwards"); + assert_eq!(error.kind(), AccessErrorKind::InvalidInterval); + assert_eq!( + finder() + .find_with(0.0, f64::NAN, |seconds| Some(reference_position(seconds))) + .expect_err("the interval is not finite") + .kind(), + AccessErrorKind::InvalidInterval + ); +} + +#[test] +fn a_higher_mask_shortens_the_pass() { + let low = AccessFinder { + station: GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 5.0_f64.to_radians()), + ..finder() + }; + let high = AccessFinder { + station: GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 30.0_f64.to_radians()), + ..finder() + }; + let samples = reference_samples(1.0, 0.0, 2_000.0); + let low_duration = low.find(&samples).expect("valid").passes[0].duration(); + let high_duration = high.find(&samples).expect("valid").passes[0].duration(); + assert!( + low_duration > high_duration, + "{low_duration} against {high_duration}" + ); +} + +#[test] +fn the_summary_is_constructible_directly() { + let pass = Pass { + rise_seconds: 0.0, + set_seconds: 10.0, + peak_seconds: 5.0, + peak_elevation: 1.0, + mask_elevation: 0.1, + peak_position: Vec3::ZERO, + }; + let summary = PassSummary::new(vec![pass], 10.0); + assert_eq!(summary.count(), 1); + assert!((summary.total_access_seconds() - 10.0).abs() < 1e-12); + assert!((summary.access_fraction().expect("valid") - 1.0).abs() < 1e-12); + assert_eq!(PassSummary::new(Vec::new(), 0.0).access_fraction(), None); +} + +#[test] +fn the_propagator_entry_point_reaches_the_tolerance_regardless_of_step() { + // When the propagator is available the refinement evaluates the true + // position at each bisection point rather than interpolating, so the + // crossing is accurate to the tolerance however coarse the sampling. + let coarse = AccessFinder { + sample_step_seconds: 30.0, + crossing_tolerance_seconds: 1e-6, + ..finder() + }; + let fine = AccessFinder { + sample_step_seconds: 1.0, + crossing_tolerance_seconds: 1e-6, + ..finder() + }; + let via_coarse = coarse + .find_with(0.0, 2_000.0, |seconds| Some(reference_position(seconds))) + .expect("valid"); + let via_fine = fine + .find_with(0.0, 2_000.0, |seconds| Some(reference_position(seconds))) + .expect("valid"); + + assert_eq!(via_coarse.count(), 1); + assert_eq!(via_fine.count(), 1); + let first = &via_coarse.passes[0]; + let second = &via_fine.passes[0]; + assert!( + (first.rise_seconds - second.rise_seconds).abs() < 1e-5, + "rise {} against {}", + first.rise_seconds, + second.rise_seconds + ); + assert!( + (first.set_seconds - second.set_seconds).abs() < 1e-5, + "set {} against {}", + first.set_seconds, + second.set_seconds + ); + // The crossing matches the closed form to the tolerance. + assert!( + (first.rise_seconds - REFERENCE_RISE).abs() < 1e-5, + "rise {} against {REFERENCE_RISE}", + first.rise_seconds + ); +}