ApsisRS
/
apsis-rs
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

network: add access window finder

This commit is contained in:
ApsisRS contributors 2026-09-26 15:51:47 +08:00
parent c75fb8bc5c
commit 2a41cf1560
7 changed files with 1193 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Cargo.lock generated
View File

@ -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"

4
crates/network/apsis-access/Cargo.toml
View File

@ -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

73
crates/network/apsis-access/src/error.rs Normal file
View File

@ -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<String>,
}
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<String>) -> 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 {}

377
crates/network/apsis-access/src/finder.rs Normal file
View File

@ -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<f64> {
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<PassSummary, AccessError> {
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<F>(
&self,
samples: &[OrbitSample],
mut elevation_at: F,
) -> Result<PassSummary, AccessError>
where
F: FnMut(f64) -> Option<f64>,
{
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<f64> = 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<F>(
&self,
samples: &[OrbitSample],
index: usize,
rising: bool,
elevation_at: &mut F,
) -> f64
where
F: FnMut(f64) -> Option<f64>,
{
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<F>(
&self,
start_seconds: f64,
end_seconds: f64,
mut propagator: F,
) -> Result<PassSummary, AccessError>
where
F: FnMut(f64) -> Option<Vec3>,
{
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,
}),
}
}
}

47
crates/network/apsis-access/src/lib.rs
View File

@ -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;

164
crates/network/apsis-access/src/pass.rs Normal file
View File

@ -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<Pass>,
/// 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<Pass>, 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<f64> {
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<f64> {
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<f64> {
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<f64> {
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)
})
}
}

523
crates/network/apsis-access/tests/finder.rs Normal file
View File

@ -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<OrbitSample> {
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
);
}