From c75fb8bc5cfe6842413ce2f08051d4964498c3f1 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 15:48:14 +0800 Subject: [PATCH] network: add coverage geometry, footprint, and visibility --- Cargo.lock | 5 + crates/network/apsis-coverage/Cargo.toml | 3 + .../network/apsis-coverage/src/footprint.rs | 176 +++++++++++ crates/network/apsis-coverage/src/geometry.rs | 135 +++++++++ crates/network/apsis-coverage/src/lib.rs | 64 +++- .../network/apsis-coverage/src/visibility.rs | 245 +++++++++++++++ .../network/apsis-coverage/tests/footprint.rs | 180 +++++++++++ .../network/apsis-coverage/tests/geometry.rs | 200 +++++++++++++ .../apsis-coverage/tests/visibility.rs | 282 ++++++++++++++++++ 9 files changed, 1289 insertions(+), 1 deletion(-) create mode 100644 crates/network/apsis-coverage/src/footprint.rs create mode 100644 crates/network/apsis-coverage/src/geometry.rs create mode 100644 crates/network/apsis-coverage/src/visibility.rs create mode 100644 crates/network/apsis-coverage/tests/footprint.rs create mode 100644 crates/network/apsis-coverage/tests/geometry.rs create mode 100644 crates/network/apsis-coverage/tests/visibility.rs diff --git a/Cargo.lock b/Cargo.lock index aec4cf6..27bc1d4 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -53,6 +53,11 @@ version = "0.1.0" [[package]] name = "apsis-coverage" version = "0.1.0" +dependencies = [ + "apsis-core", + "apsis-geodesy", + "apsis-math", +] [[package]] name = "apsis-doppler" diff --git a/crates/network/apsis-coverage/Cargo.toml b/crates/network/apsis-coverage/Cargo.toml index 4a01c96..cd76779 100644 --- a/crates/network/apsis-coverage/Cargo.toml +++ b/crates/network/apsis-coverage/Cargo.toml @@ -11,3 +11,6 @@ publish = false workspace = true [dependencies] +apsis-core.workspace = true +apsis-math.workspace = true +apsis-geodesy.workspace = true diff --git a/crates/network/apsis-coverage/src/footprint.rs b/crates/network/apsis-coverage/src/footprint.rs new file mode 100644 index 0000000..869672c --- /dev/null +++ b/crates/network/apsis-coverage/src/footprint.rs @@ -0,0 +1,176 @@ +//! Footprint geometry and coverage fractions. + +use crate::geometry::{earth_central_angle, elevation_from_central_angle}; +use apsis_math::Vec3; +use core::f64::consts::PI; + +/// A circular ground footprint. +/// +/// The footprint is the part of the Earth's surface from which a satellite +/// is visible above a given elevation mask. On a sphere it is an exact +/// circle, described here by the central angle from the sub-satellite point +/// and the radius of the body. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Footprint { + /// The central half angle of the footprint, in radians. + /// + /// This is the angle at the centre of the Earth from the sub-satellite + /// point to the edge of the footprint. + pub central_angle: f64, + /// The radius of the body, in kilometers. + pub body_radius: f64, +} + +impl Footprint { + /// Creates a footprint from its central angle and body radius. + pub const fn new(central_angle: f64, body_radius: f64) -> Self { + Self { + central_angle, + body_radius, + } + } + + /// Returns the footprint radius measured along the surface, in + /// kilometers. + pub fn surface_radius(&self) -> f64 { + self.central_angle * self.body_radius + } + + /// Returns the area of the footprint, in square kilometers. + /// + /// The area is the part of a sphere cut off by a cone, + /// + /// ```text + /// A = 2 pi R^2 (1 - cos(central)) + /// ``` + /// + /// which is exact for a spherical body. + pub fn area(&self) -> f64 { + 2.0 * PI * self.body_radius * self.body_radius * (1.0 - self.central_angle.cos()) + } + + /// Returns the fraction of the whole sphere covered by the footprint. + /// + /// The value is between zero and one half, because a footprint from a + /// point above the surface can never cover more than a hemisphere. + pub fn fraction_of_surface(&self) -> f64 { + self.area() / (4.0 * PI * self.body_radius * self.body_radius) + } + + /// Returns `true` when a ground point at the given central angle from + /// the sub-satellite point lies inside the footprint. + /// + /// The comparison is inclusive, so a point exactly on the edge is + /// inside. + pub fn contains_central_angle(&self, central_angle: f64) -> bool { + central_angle.is_finite() && central_angle >= 0.0 && central_angle <= self.central_angle + } +} + +/// Returns the central half angle of the footprint, in radians. +/// +/// This is the largest central angle at which the satellite is still above +/// the elevation mask. +/// +/// Returns `None` when the geometry is not usable. +pub fn footprint_central_angle( + satellite_radius: f64, + body_radius: f64, + elevation_mask: f64, +) -> Option { + let central = earth_central_angle(satellite_radius, body_radius, elevation_mask)?; + if central.is_finite() && central >= 0.0 { + Some(central) + } else { + None + } +} + +/// Returns the footprint for a satellite and an elevation mask. +/// +/// Returns `None` when the geometry is not usable. +pub fn footprint( + satellite_radius: f64, + body_radius: f64, + elevation_mask: f64, +) -> Option { + Some(Footprint::new( + footprint_central_angle(satellite_radius, body_radius, elevation_mask)?, + body_radius, + )) +} + +/// Returns the area of a footprint, in square kilometers. +/// +/// This is a free function form of [`Footprint::area`]. +/// +/// Returns `None` when the geometry is not usable. +pub fn footprint_area(satellite_radius: f64, body_radius: f64, elevation_mask: f64) -> Option { + Some(footprint(satellite_radius, body_radius, elevation_mask)?.area()) +} + +/// Returns the fraction of the Earth's surface covered by a footprint. +/// +/// Returns `None` when the geometry is not usable. +pub fn footprint_fraction( + satellite_radius: f64, + body_radius: f64, + elevation_mask: f64, +) -> Option { + Some(footprint(satellite_radius, body_radius, elevation_mask)?.fraction_of_surface()) +} + +/// Returns the angle between two surface points, in radians. +/// +/// The vectors must be positions relative to the centre of the body, in the +/// same frame. The angle is the central angle between them, as seen from the +/// centre. +/// +/// Returns `None` when either vector has zero length or a non-finite +/// component. +pub fn central_angle_between(first: Vec3, second: Vec3) -> Option { + if !first.is_finite() || !second.is_finite() { + return None; + } + let first_norm = first.norm(); + let second_norm = second.norm(); + if first_norm == 0.0 || second_norm == 0.0 { + return None; + } + Some( + (first.dot(second) / (first_norm * second_norm)) + .clamp(-1.0, 1.0) + .acos(), + ) +} + +/// Returns `true` when a ground point lies inside a footprint. +/// +/// The point and the sub-satellite point are positions relative to the +/// centre of the body, in kilometers. The footprint supplies the body radius +/// and the central half angle. +/// +/// Returns `None` when either vector is degenerate or not finite. +pub fn point_in_footprint( + footprint: &Footprint, + point: Vec3, + sub_satellite_point: Vec3, +) -> Option { + let central = central_angle_between(point, sub_satellite_point)?; + Some(footprint.contains_central_angle(central)) +} + +/// Returns the elevation at which a footprint edge is seen from the ground, +/// in radians. +/// +/// This is the inverse of [`footprint_central_angle`], useful when a +/// footprint is known and the mask it corresponds to is wanted. +/// +/// Returns `None` when the geometry is not usable. +pub fn footprint_edge_elevation(satellite_radius: f64, footprint: &Footprint) -> Option { + elevation_from_central_angle( + satellite_radius, + footprint.body_radius, + footprint.central_angle, + ) +} diff --git a/crates/network/apsis-coverage/src/geometry.rs b/crates/network/apsis-coverage/src/geometry.rs new file mode 100644 index 0000000..0e69204 --- /dev/null +++ b/crates/network/apsis-coverage/src/geometry.rs @@ -0,0 +1,135 @@ +//! The spherical geometry of a satellite, a ground point, and the Earth's +//! centre. +//! +//! The three points form a triangle. The angles of that triangle are the +//! elevation at the ground point, the nadir angle at the satellite, and the +//! central angle at the centre of the Earth. Because the ground point is +//! directly below the local horizon, the elevation and the adjacent angle at +//! the ground point are complementary, so the three angles satisfy +//! +//! ```text +//! elevation + nadir + central = 90 degrees +//! ``` +//! +//! The quantities of interest follow from the law of sines applied to that +//! triangle: +//! +//! ```text +//! sin(nadir) / body_radius = cos(elevation) / satellite_radius +//! ``` +//! +//! where the cosine appears because the angle at the ground point, inside +//! the triangle, is the complement of the elevation. +//! +//! All angles are in radians and all lengths in kilometers. The body radius +//! is the distance from the centre of the Earth to the ground point, so a +//! station above sea level should pass its own radius rather than the +//! equatorial radius. + +use core::f64::consts::FRAC_PI_2; + +/// Returns the nadir angle at the satellite, in radians. +/// +/// The nadir angle is measured at the satellite between the direction to the +/// centre of the Earth and the direction to the ground point. It is zero when +/// the ground point is directly below the satellite and grows as the point +/// moves toward the horizon. +/// +/// Returns `None` when the radii are not usable, when the body radius +/// exceeds the satellite radius, or when the elevation is not finite. +pub fn nadir_angle(satellite_radius: f64, body_radius: f64, elevation: f64) -> Option { + if !satellite_radius.is_finite() || !body_radius.is_finite() || !elevation.is_finite() { + return None; + } + if body_radius <= 0.0 || satellite_radius < body_radius { + return None; + } + let sine = body_radius * elevation.cos() / satellite_radius; + // Rounding can push the sine a little above one for a grazing geometry. + Some(sine.clamp(-1.0, 1.0).asin()) +} + +/// Returns the earth central angle between the sub-satellite point and the +/// ground point, in radians. +/// +/// The central angle is the angle at the centre of the Earth between the two +/// radial directions. It is the distance along the surface, divided by the +/// body radius. +/// +/// Returns `None` when the geometry is not usable. +pub fn earth_central_angle(satellite_radius: f64, body_radius: f64, elevation: f64) -> Option { + let nadir = nadir_angle(satellite_radius, body_radius, elevation)?; + // Rounding can leave a small negative value at the zenith, where the + // true angle is zero. The negative part is not physical, so it is + // removed before the result is returned. + Some((FRAC_PI_2 - elevation - nadir).max(0.0)) +} + +/// Returns the slant range between the satellite and the ground point, in +/// kilometers. +/// +/// The range is the straight line distance, not the distance along the +/// surface. It is the quantity a link budget needs, because it sets the free +/// space path loss. +/// +/// The relation is the law of cosines applied to the triangle formed by the +/// satellite, the ground point, and the centre of the Earth, +/// +/// ```text +/// range^2 = satellite_radius^2 + body_radius^2 +/// - 2 satellite_radius body_radius cos(central) +/// ``` +/// +/// The law of cosines is used rather than the equivalent form +/// `Rs sin(central) / cos(elevation)` because the latter is singular at the +/// zenith, where the elevation is a right angle and its cosine is zero, +/// even though the range there is simply the altitude. +/// +/// Returns `None` when the geometry is not usable. +pub fn slant_range(satellite_radius: f64, body_radius: f64, elevation: f64) -> Option { + let central = earth_central_angle(satellite_radius, body_radius, elevation)?; + if !central.is_finite() || central < 0.0 { + return None; + } + let squared = satellite_radius * satellite_radius + body_radius * body_radius + - 2.0 * satellite_radius * body_radius * central.cos(); + if squared < 0.0 { + return None; + } + Some(squared.sqrt()) +} + +/// Returns the elevation angle for a central angle, in radians. +/// +/// This inverts [`earth_central_angle`]. It is used when the position of a +/// ground point relative to the sub-satellite point is known and the +/// elevation at that point is wanted. +/// +/// The relation is +/// +/// ```text +/// tan(elevation) = (cos(central) - body_radius / satellite_radius) +/// / sin(central) +/// ``` +/// +/// Returns `None` when the radii are not usable or when the central angle is +/// not finite. +pub fn elevation_from_central_angle( + satellite_radius: f64, + body_radius: f64, + central_angle: f64, +) -> Option { + if !satellite_radius.is_finite() || !body_radius.is_finite() || !central_angle.is_finite() { + return None; + } + if body_radius <= 0.0 || satellite_radius < body_radius { + return None; + } + let ratio = body_radius / satellite_radius; + let sine = central_angle.sin(); + if sine == 0.0 { + // Directly below the satellite the elevation is a right angle. + return Some(FRAC_PI_2); + } + Some(((central_angle.cos() - ratio) / sine).atan()) +} diff --git a/crates/network/apsis-coverage/src/lib.rs b/crates/network/apsis-coverage/src/lib.rs index 554fb36..bcdf332 100644 --- a/crates/network/apsis-coverage/src/lib.rs +++ b/crates/network/apsis-coverage/src/lib.rs @@ -1 +1,63 @@ -//! Ground and regional coverage analysis. +//! Ground coverage analysis. +//! +//! A satellite covers the part of the Earth's surface from which it can be +//! seen above a minimum elevation angle. That part is a circle, called the +//! footprint, whose size depends on the satellite altitude and on the +//! elevation mask. This crate computes the footprint geometry, tests whether +//! a ground point lies inside a footprint, and measures how much of a region +//! is covered. +//! +//! The spherical model in [`geometry`](self) treats the Earth as a sphere. +//! That is accurate enough for constellation design, where the quantities +//! wanted are the footprint radius to a few parts in ten thousand and the +//! number of satellites in view. The exact line-of-sight test in +//! [`line_of_sight`] uses the actual positions, so it accounts for the +//! flattening and for a ground point above the ellipsoid. +//! +//! The geometry of the satellite, the ground point, and the centre of the +//! Earth is a triangle. Its three angles are the elevation at the ground +//! point, the nadir angle at the satellite, and the central angle at the +//! centre, and they sum to a right angle: +//! +//! ```text +//! elevation + nadir + central = 90 degrees +//! ``` +//! +//! Source: the relations used here are the standard satellite geometry, as +//! given for example in Vallado, D. A., *Fundamentals of Astrodynamics and +//! Applications*, 4th edition, microcosm Press, 2013, section 5.3. +//! +//! ``` +//! use apsis_coverage::{footprint_central_angle, slant_range}; +//! +//! // A satellite at 550 km altitude sees to about fifteen degrees of +//! // central angle at a ten degree elevation mask, which is a footprint +//! // radius of some 1665 km along the surface. +//! let central = footprint_central_angle(6_928.137, 6_378.137, 10.0_f64.to_radians()) +//! .expect("valid geometry"); +//! assert!((central.to_degrees() - 14.96).abs() < 0.01); +//! +//! let range = slant_range(6_928.137, 6_378.137, 10.0_f64.to_radians()) +//! .expect("valid geometry"); +//! assert!((range - 1_815.65).abs() < 0.01); +//! ``` + +#![forbid(unsafe_code)] +#![deny(missing_docs)] + +mod footprint; +mod geometry; +mod visibility; + +pub use footprint::{ + Footprint, central_angle_between, footprint, footprint_area, footprint_central_angle, + footprint_edge_elevation, footprint_fraction, point_in_footprint, +}; +pub use geometry::{earth_central_angle, elevation_from_central_angle, nadir_angle, slant_range}; +pub use visibility::{ + GroundPoint, Sensor, SubSatellitePoint, central_angle_to_station, elevation_angle, + line_of_sight, local_up, sub_satellite_point, visible_from, visible_from_geodetic, +}; + +pub use apsis_geodesy::{Ellipsoid, Geodetic, WGS84}; +pub use apsis_math::Vec3; diff --git a/crates/network/apsis-coverage/src/visibility.rs b/crates/network/apsis-coverage/src/visibility.rs new file mode 100644 index 0000000..5e63bc4 --- /dev/null +++ b/crates/network/apsis-coverage/src/visibility.rs @@ -0,0 +1,245 @@ +//! Line of sight and visibility from a ground point. + +use crate::geometry::earth_central_angle; +use apsis_geodesy::{Ellipsoid, Geodetic, WGS84}; +use apsis_math::Vec3; +use core::f64::consts::FRAC_PI_2; + +/// A ground station, given as a geodetic position on an ellipsoid. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct GroundPoint { + /// The geodetic position of the station. + pub position: Geodetic, + /// The ellipsoid the position is referred to. + pub ellipsoid: Ellipsoid, + /// The sensor elevation mask, in radians. + /// + /// The minimum elevation above the local horizon at which the station + /// can see a satellite. A mask of zero means the station sees down to + /// its own horizon. + pub elevation_mask: f64, +} + +impl GroundPoint { + /// Creates a ground point on the WGS 84 ellipsoid with a given mask. + pub fn new(position: Geodetic, elevation_mask: f64) -> Self { + Self { + position, + ellipsoid: WGS84, + elevation_mask, + } + } + + /// Returns the position of the station in the Earth-fixed Cartesian + /// frame, in kilometers. + /// + /// The geodetic conversion yields meters, so the result is divided by a + /// thousand to match the frame the orbit layer uses. + pub fn ecef(&self) -> Vec3 { + self.position.to_cartesian(self.ellipsoid) / 1_000.0 + } + + /// Returns the distance from the centre of the Earth, in kilometers. + /// + /// This is the local radius of the ellipsoid at the station's latitude + /// plus its height, so it includes the effect of the flattening. + pub fn radius(&self) -> f64 { + self.ecef().norm() + } + + /// Returns `true` when the station and its mask describe a usable + /// sensor. + pub fn is_valid(&self) -> bool { + self.position.is_valid() + && self.elevation_mask.is_finite() + && (-FRAC_PI_2..=FRAC_PI_2).contains(&self.elevation_mask) + && self.ellipsoid.semi_major_axis() > 0.0 + } +} + +/// A sensor at a ground point, with a pointing direction. +/// +/// The direction is a unit vector in the Earth-fixed frame. It is optional, +/// because a station that can point anywhere with sufficient gain is +/// modelled by its elevation mask alone; the direction matters only when the +/// station has a fixed beam. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Sensor { + /// The station. + pub ground: GroundPoint, + /// The pointing direction in the Earth-fixed frame, if the station is + /// directional. + pub pointing: Option, +} + +impl Sensor { + /// Creates a station that can point anywhere above its mask. + pub const fn omnidirectional(ground: GroundPoint) -> Self { + Self { + ground, + pointing: None, + } + } + + /// Creates a station with a fixed pointing direction. + /// + /// The direction is normalized when it is used, so the caller may pass + /// an unnormalized vector. + pub const fn directional(ground: GroundPoint, pointing: Vec3) -> Self { + Self { + ground, + pointing: Some(pointing), + } + } +} + +/// The point on the Earth's surface directly below a satellite. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct SubSatellitePoint { + /// The geodetic position of the sub-satellite point. + pub position: Geodetic, + /// The ellipsoid the position is referred to. + pub ellipsoid: Ellipsoid, +} + +/// Returns the sub-satellite point for a position in the Earth-fixed frame. +/// +/// The satellite position is in kilometers. The sub-satellite point is the +/// intersection of the line from the centre of the Earth through the +/// satellite with the ellipsoid, which is the geodetic position of the +/// satellite with its height set to zero at the same latitude and longitude. +/// +/// This is the geocentric projection, not the point at which the ellipsoid +/// normal through the satellite meets the surface. The difference is at most +/// about twenty kilometers for a low Earth satellite, and the geocentric +/// form is the one the footprint geometry uses. +/// +/// Returns `None` when the position is not finite or lies on the polar axis, +/// where no longitude is defined. +pub fn sub_satellite_point( + position_ecef_km: Vec3, + ellipsoid: Ellipsoid, +) -> Option { + if !position_ecef_km.is_finite() { + return None; + } + let position_m = position_ecef_km * 1_000.0; + let geodetic = apsis_geodesy::geodetic_from_cartesian(position_m, ellipsoid)?; + Some(SubSatellitePoint { + position: Geodetic::new(geodetic.latitude, geodetic.longitude, 0.0), + ellipsoid, + }) +} + +/// Returns the elevation of a satellite as seen from a station, in radians. +/// +/// The satellite position is in the Earth-fixed frame, in kilometers. The +/// elevation is measured from the local horizon, which is the plane +/// perpendicular to the ellipsoid normal at the station. +/// +/// A positive value means the satellite is above the horizon, a value of a +/// right angle means it is at the zenith, and a negative value means it is +/// below the horizon. +/// +/// Returns `None` when either position is not usable. +pub fn elevation_angle(satellite_ecef_km: Vec3, ground: &GroundPoint) -> Option { + if !satellite_ecef_km.is_finite() || !ground.is_valid() { + return None; + } + let station = ground.ecef(); + let line_of_sight = satellite_ecef_km - station; + if line_of_sight.norm() == 0.0 { + return None; + } + // The local up direction is the ellipsoid normal at the station, which + // for a geodetic position is the direction the surface normal points. + let up = local_up(ground); + let sine = line_of_sight.normalize()?.dot(up.normalize()?); + Some(sine.clamp(-1.0, 1.0).asin()) +} + +/// Returns the unit vector along the ellipsoid normal at a station. +/// +/// This is the local vertical, the direction a plumb line points. For a +/// geodetic position it follows from the latitude and longitude, because the +/// normal of a biaxial ellipsoid lies in the meridian plane and makes the +/// geodetic latitude with the equatorial plane. +pub fn local_up(ground: &GroundPoint) -> Vec3 { + let latitude = ground.position.latitude; + let longitude = ground.position.longitude; + let (sin_latitude, cos_latitude) = latitude.sin_cos(); + let (sin_longitude, cos_longitude) = longitude.sin_cos(); + Vec3::new( + cos_latitude * cos_longitude, + cos_latitude * sin_longitude, + sin_latitude, + ) +} + +/// Returns the elevation of a satellite as seen from a station, in radians. +/// +/// This is the same as [`elevation_angle`], named for symmetry with the +/// other visibility functions. +/// +/// Returns `None` when either position is not usable. +pub fn line_of_sight(satellite_ecef_km: Vec3, ground: &GroundPoint) -> Option { + elevation_angle(satellite_ecef_km, ground) +} + +/// Returns `true` when a satellite is visible from a sensor. +/// +/// The satellite is visible when its elevation at the station is at least +/// the sensor's mask, and, for a directional sensor, when the direction to +/// the satellite is within a right angle of the pointing direction, so that +/// the satellite is in the forward hemisphere of the beam. +/// +/// Returns `None` when the satellite position or the sensor is not usable. +pub fn visible_from(satellite_ecef_km: Vec3, sensor: &Sensor) -> Option { + let elevation = elevation_angle(satellite_ecef_km, &sensor.ground)?; + if elevation < sensor.ground.elevation_mask { + return Some(false); + } + match sensor.pointing { + None => Some(true), + Some(pointing) => { + let to_satellite = satellite_ecef_km - sensor.ground.ecef(); + let pointing = pointing.normalize()?; + let to_satellite = to_satellite.normalize()?; + Some(pointing.dot(to_satellite) > 0.0) + } + } +} + +/// Returns the elevation of a satellite as seen from a geodetic station. +/// +/// This is a convenience for a station given directly as a geodetic +/// position, without building a [`GroundPoint`]. +/// +/// Returns `None` when the satellite position or the station is not usable. +pub fn visible_from_geodetic( + satellite_ecef_km: Vec3, + position: Geodetic, + elevation_mask: f64, +) -> Option { + let ground = GroundPoint::new(position, elevation_mask); + visible_from(satellite_ecef_km, &Sensor::omnidirectional(ground)) +} + +/// Returns the central angle between a satellite and a ground station, in +/// radians. +/// +/// This uses the spherical relation between the two radii and the elevation, +/// so it is consistent with the footprint geometry. The satellite position +/// is used only for its radius. +/// +/// Returns `None` when the geometry is not usable. +pub fn central_angle_to_station( + satellite_ecef_km: Vec3, + ground: &GroundPoint, + elevation: f64, +) -> Option { + if !satellite_ecef_km.is_finite() || !ground.is_valid() { + return None; + } + earth_central_angle(satellite_ecef_km.norm(), ground.radius(), elevation) +} diff --git a/crates/network/apsis-coverage/tests/footprint.rs b/crates/network/apsis-coverage/tests/footprint.rs new file mode 100644 index 0000000..d4e8cd9 --- /dev/null +++ b/crates/network/apsis-coverage/tests/footprint.rs @@ -0,0 +1,180 @@ +//! Tests for footprint geometry and coverage fractions. + +use apsis_coverage::{ + Footprint, central_angle_between, footprint, footprint_area, footprint_central_angle, + footprint_edge_elevation, footprint_fraction, point_in_footprint, +}; +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; + +/// A low Earth orbit at 550 km altitude. +const LOW_ORBIT: f64 = BODY_RADIUS + 550.0; + +/// A geostationary radius. +const GEOSTATIONARY: f64 = 42_164.0; + +const TOLERANCE: f64 = 1e-9; + +#[test] +fn the_footprint_at_a_ten_degree_mask_matches_a_hand_computation() { + // A 550 km satellite at a ten degree mask has a central angle of + // 14.9566 degrees, a surface radius of 1664.965 km, and an area of + // 8 659 493.8 square kilometers. + let mask = 10.0_f64.to_radians(); + let footprint = footprint(LOW_ORBIT, BODY_RADIUS, mask).expect("valid"); + assert!((footprint.central_angle.to_degrees() - 14.956_635_124_957_29).abs() < 1e-9); + assert!((footprint.surface_radius() - 1_664.965_006_091_034_9).abs() < 1e-6); + assert!((footprint.area() - 8_659_493.809_336_036).abs() < 1e-3); + assert!((footprint.fraction_of_surface() - 0.016_939_280_322_558_092).abs() < 1e-12); +} + +#[test] +fn the_area_matches_the_spherical_cap_relation() { + let footprint = footprint(LOW_ORBIT, BODY_RADIUS, 0.3).expect("valid"); + let expected = 2.0 * PI * BODY_RADIUS * BODY_RADIUS * (1.0 - footprint.central_angle.cos()); + assert!((footprint.area() - expected).abs() < TOLERANCE); +} + +#[test] +fn the_area_grows_with_the_central_angle() { + let mut previous = 0.0; + for mask_degrees in [30.0_f64, 20.0, 10.0, 5.0, 0.0] { + let area = + footprint_area(LOW_ORBIT, BODY_RADIUS, mask_degrees.to_radians()).expect("valid"); + assert!( + area > previous, + "mask {mask_degrees}: {area} against {previous}" + ); + previous = area; + } +} + +#[test] +fn a_footprint_can_never_cover_more_than_a_hemisphere() { + // The limit is a satellite at the surface with no mask, which sees half + // the sphere. Any real satellite sees less. + for satellite_radius in [LOW_ORBIT, GEOSTATIONARY, 1.0e6] { + for mask_degrees in [0.0_f64, 10.0, 45.0] { + let fraction = + footprint_fraction(satellite_radius, BODY_RADIUS, mask_degrees.to_radians()) + .expect("valid"); + assert!( + fraction > 0.0 && fraction <= 0.5, + "radius {satellite_radius}, mask {mask_degrees}: {fraction}" + ); + } + } +} + +#[test] +fn a_geostationary_satellite_sees_about_a_third_of_the_earth() { + // A geostationary satellite at a zero degree mask sees to about 81.3 + // degrees of central angle, which is a fraction of about 0.42. + let fraction = footprint_fraction(GEOSTATIONARY, BODY_RADIUS, 0.0).expect("valid"); + assert!(fraction > 0.40 && fraction < 0.43, "fraction {fraction}"); +} + +#[test] +fn containment_is_inclusive_at_the_edge() { + let footprint = Footprint::new(0.5, BODY_RADIUS); + assert!(footprint.contains_central_angle(0.0)); + assert!(footprint.contains_central_angle(0.25)); + assert!(footprint.contains_central_angle(0.5)); + assert!(!footprint.contains_central_angle(0.500_001)); + assert!(!footprint.contains_central_angle(-0.1)); + assert!(!footprint.contains_central_angle(f64::NAN)); +} + +#[test] +fn the_central_angle_between_two_surface_points_is_the_included_angle() { + // Two points on the equator ninety degrees apart. + let first = Vec3::new(BODY_RADIUS, 0.0, 0.0); + let second = Vec3::new(0.0, BODY_RADIUS, 0.0); + let angle = central_angle_between(first, second).expect("valid"); + assert!((angle - FRAC_PI_2).abs() < TOLERANCE); + + // A point with itself gives zero. + assert!(central_angle_between(first, first).expect("valid").abs() < TOLERANCE); + + // Opposite points give a straight angle. + let opposite = Vec3::new(-BODY_RADIUS, 0.0, 0.0); + assert!((central_angle_between(first, opposite).expect("valid") - PI).abs() < TOLERANCE); +} + +#[test] +fn the_central_angle_between_points_is_scale_invariant() { + // The angle depends only on the directions, not on the radii. + let near = Vec3::new(1.0, 0.5, 0.25); + let far = near * 1_000.0; + let other = Vec3::new(0.0, 1.0, 0.0); + let first = central_angle_between(near, other).expect("valid"); + let second = central_angle_between(far, other).expect("valid"); + assert!((first - second).abs() < 1e-12); +} + +#[test] +fn a_degenerate_vector_has_no_central_angle() { + assert_eq!(central_angle_between(Vec3::ZERO, Vec3::X), None); + assert_eq!(central_angle_between(Vec3::X, Vec3::ZERO), None); + assert_eq!(central_angle_between(Vec3::splat(f64::NAN), Vec3::X), None); +} + +#[test] +fn point_in_footprint_uses_the_central_angle_from_the_sub_satellite_point() { + let footprint = Footprint::new(0.1, BODY_RADIUS); + // The sub-satellite point on the equator. + let sub_satellite = Vec3::new(BODY_RADIUS, 0.0, 0.0); + + // A point within the footprint, about 0.05 radians away. + let inside = Vec3::new( + BODY_RADIUS * 0.05_f64.cos(), + BODY_RADIUS * 0.05_f64.sin(), + 0.0, + ); + assert_eq!( + point_in_footprint(&footprint, inside, sub_satellite), + Some(true) + ); + + // A point well outside it, about a quarter radian away. + let outside = Vec3::new( + BODY_RADIUS * 0.25_f64.cos(), + BODY_RADIUS * 0.25_f64.sin(), + 0.0, + ); + assert_eq!( + point_in_footprint(&footprint, outside, sub_satellite), + Some(false) + ); +} + +#[test] +fn the_edge_elevation_inverts_the_footprint_central_angle() { + let mask = 12.0_f64.to_radians(); + let footprint = footprint(LOW_ORBIT, BODY_RADIUS, mask).expect("valid"); + let recovered = footprint_edge_elevation(LOW_ORBIT, &footprint).expect("valid"); + assert!( + (recovered - mask).abs() < 1e-9, + "{recovered} against {mask}" + ); +} + +#[test] +fn the_footprint_central_angle_helper_matches_the_constructor() { + let mask = 15.0_f64.to_radians(); + let central = footprint_central_angle(LOW_ORBIT, BODY_RADIUS, mask).expect("valid"); + let footprint = footprint(LOW_ORBIT, BODY_RADIUS, mask).expect("valid"); + assert!((central - footprint.central_angle).abs() < TOLERANCE); +} + +#[test] +fn invalid_geometry_is_rejected() { + assert_eq!(footprint(LOW_ORBIT, 0.0, 0.1), None); + assert_eq!(footprint(BODY_RADIUS - 1.0, BODY_RADIUS, 0.1), None); + assert_eq!(footprint(LOW_ORBIT, BODY_RADIUS, f64::NAN), None); + assert_eq!(footprint_area(LOW_ORBIT, BODY_RADIUS, f64::INFINITY), None); + assert_eq!(footprint_fraction(LOW_ORBIT, -1.0, 0.1), None); +} diff --git a/crates/network/apsis-coverage/tests/geometry.rs b/crates/network/apsis-coverage/tests/geometry.rs new file mode 100644 index 0000000..9dad599 --- /dev/null +++ b/crates/network/apsis-coverage/tests/geometry.rs @@ -0,0 +1,200 @@ +//! Tests for the satellite geometry of a ground coverage calculation. +//! +//! The expected values are computed independently from the law of sines +//! applied to the triangle formed by the satellite, the ground point, and +//! the centre of the Earth. + +use apsis_coverage::{earth_central_angle, elevation_from_central_angle, nadir_angle, slant_range}; +use core::f64::consts::FRAC_PI_2; + +/// The Earth equatorial radius, in kilometers. +const BODY_RADIUS: f64 = 6_378.137; + +/// A low Earth orbit at 550 km altitude. +const LOW_ORBIT: f64 = BODY_RADIUS + 550.0; + +/// A medium Earth orbit at 20 200 km altitude, near the navigation +/// constellations. +const MEDIUM_ORBIT: f64 = BODY_RADIUS + 20_200.0; + +/// A geostationary radius. +const GEOSTATIONARY: f64 = 42_164.0; + +const TOLERANCE: f64 = 1e-9; + +#[test] +fn the_angles_at_the_zenith_are_consistent() { + // Directly below the satellite the elevation is a right angle, the + // nadir angle is zero, and the central angle is zero. + assert!( + nadir_angle(LOW_ORBIT, BODY_RADIUS, FRAC_PI_2) + .expect("valid") + .abs() + < TOLERANCE + ); + assert!( + earth_central_angle(LOW_ORBIT, BODY_RADIUS, FRAC_PI_2) + .expect("valid") + .abs() + < TOLERANCE + ); + // The slant range is the altitude. + assert!( + (slant_range(LOW_ORBIT, BODY_RADIUS, FRAC_PI_2).expect("valid") - 550.0).abs() < TOLERANCE + ); +} + +#[test] +fn the_three_angles_are_complementary() { + // elevation + nadir + central = 90 degrees, for every geometry. + for satellite_radius in [LOW_ORBIT, MEDIUM_ORBIT, GEOSTATIONARY] { + for elevation_degrees in [0.0_f64, 5.0, 10.0, 30.0, 60.0, 89.0] { + let elevation = elevation_degrees.to_radians(); + let nadir = nadir_angle(satellite_radius, BODY_RADIUS, elevation).expect("valid"); + let central = + earth_central_angle(satellite_radius, BODY_RADIUS, elevation).expect("valid"); + let sum = elevation + nadir + central; + assert!( + (sum - FRAC_PI_2).abs() < 1e-9, + "radius {satellite_radius}, elevation {elevation_degrees}: sum {sum}" + ); + } + } +} + +#[test] +fn the_geometry_matches_a_hand_computation() { + // A 550 km satellite at a ten degree mask. The nadir angle is + // asin(Re cos(el) / Rs) = 65.0434 degrees, the central angle is + // 14.9566 degrees, and the slant range is 1815.652 km. + let elevation = 10.0_f64.to_radians(); + assert!( + (nadir_angle(LOW_ORBIT, BODY_RADIUS, elevation) + .expect("valid") + .to_degrees() + - 65.043_364_875_042_71) + .abs() + < 1e-9 + ); + assert!( + (earth_central_angle(LOW_ORBIT, BODY_RADIUS, elevation) + .expect("valid") + .to_degrees() + - 14.956_635_124_957_29) + .abs() + < 1e-9 + ); + assert!( + (slant_range(LOW_ORBIT, BODY_RADIUS, elevation).expect("valid") - 1_815.652_169_368_375_8) + .abs() + < 1e-9 + ); +} + +#[test] +fn a_higher_satellite_sees_further() { + let elevation = 10.0_f64.to_radians(); + let low = earth_central_angle(LOW_ORBIT, BODY_RADIUS, elevation).expect("valid"); + let medium = earth_central_angle(MEDIUM_ORBIT, BODY_RADIUS, elevation).expect("valid"); + let high = earth_central_angle(GEOSTATIONARY, BODY_RADIUS, elevation).expect("valid"); + assert!(low < medium, "{low} against {medium}"); + assert!(medium < high, "{medium} against {high}"); +} + +#[test] +fn a_lower_elevation_mask_sees_further() { + // Relaxing the mask lets the satellite be seen closer to the horizon, + // which increases the central angle. + let mut previous = 0.0; + for elevation_degrees in [30.0_f64, 20.0, 10.0, 5.0, 0.0] { + let central = earth_central_angle(LOW_ORBIT, BODY_RADIUS, elevation_degrees.to_radians()) + .expect("valid"); + assert!( + central > previous, + "elevation {elevation_degrees}: {central} against {previous}" + ); + previous = central; + } +} + +#[test] +fn a_zero_degree_mask_is_the_horizon() { + // At a zero degree mask the central angle is the horizon angle + // acos(Re / Rs), which is where the line of sight is tangent to the + // surface. + let central = earth_central_angle(LOW_ORBIT, BODY_RADIUS, 0.0).expect("valid"); + let horizon = (BODY_RADIUS / LOW_ORBIT).acos(); + assert!((central - horizon).abs() < 1e-9); +} + +#[test] +fn the_slant_range_grows_as_the_elevation_falls() { + // A satellite near the horizon is further away than one overhead. + let overhead = slant_range(LOW_ORBIT, BODY_RADIUS, FRAC_PI_2).expect("valid"); + let low = slant_range(LOW_ORBIT, BODY_RADIUS, 5.0_f64.to_radians()).expect("valid"); + assert!(low > overhead, "{low} against {overhead}"); + assert!(overhead < 600.0); + assert!(low > 1_000.0); +} + +#[test] +fn the_elevation_inverts_the_central_angle() { + for satellite_radius in [LOW_ORBIT, MEDIUM_ORBIT, GEOSTATIONARY] { + for elevation_degrees in [0.0_f64, 10.0, 30.0, 70.0] { + let elevation = elevation_degrees.to_radians(); + let central = + earth_central_angle(satellite_radius, BODY_RADIUS, elevation).expect("valid"); + let recovered = elevation_from_central_angle(satellite_radius, BODY_RADIUS, central) + .expect("valid"); + assert!( + (recovered - elevation).abs() < 1e-9, + "radius {satellite_radius}, elevation {elevation_degrees}: {recovered}" + ); + } + } +} + +#[test] +fn the_elevation_is_a_right_angle_when_the_central_angle_is_zero() { + let elevation = elevation_from_central_angle(LOW_ORBIT, BODY_RADIUS, 0.0).expect("valid"); + assert!((elevation - FRAC_PI_2).abs() < TOLERANCE); +} + +#[test] +fn a_satellite_inside_the_body_is_rejected() { + // A satellite below the surface has no geometry, so the functions + // report nothing rather than a negative angle. + assert_eq!(nadir_angle(BODY_RADIUS - 1.0, BODY_RADIUS, 0.0), None); + assert_eq!( + earth_central_angle(BODY_RADIUS - 1.0, BODY_RADIUS, 0.0), + None + ); + assert_eq!(slant_range(BODY_RADIUS - 1.0, BODY_RADIUS, 0.0), None); +} + +#[test] +fn non_finite_arguments_are_rejected() { + assert_eq!(nadir_angle(f64::NAN, BODY_RADIUS, 0.0), None); + assert_eq!(nadir_angle(LOW_ORBIT, f64::NAN, 0.0), None); + assert_eq!(nadir_angle(LOW_ORBIT, BODY_RADIUS, f64::NAN), None); + assert_eq!( + earth_central_angle(LOW_ORBIT, BODY_RADIUS, f64::INFINITY), + None + ); + assert_eq!(slant_range(LOW_ORBIT, 0.0, 0.0), None); + assert_eq!( + elevation_from_central_angle(LOW_ORBIT, BODY_RADIUS, f64::NAN), + None + ); +} + +#[test] +fn a_grazing_geometry_is_handled_at_the_horizon() { + // At the horizon the nadir angle is the arc sine of the ratio of the + // radii, which is exactly one when the satellite is at the surface and + // less than one above it. The clamp keeps the arc sine defined. + let nadir = nadir_angle(LOW_ORBIT, BODY_RADIUS, 0.0).expect("valid"); + assert!(nadir.is_finite()); + let expected = (BODY_RADIUS / LOW_ORBIT).asin(); + assert!((nadir - expected).abs() < 1e-9); +} diff --git a/crates/network/apsis-coverage/tests/visibility.rs b/crates/network/apsis-coverage/tests/visibility.rs new file mode 100644 index 0000000..095728c --- /dev/null +++ b/crates/network/apsis-coverage/tests/visibility.rs @@ -0,0 +1,282 @@ +//! Tests for line of sight and visibility from a ground station. + +use apsis_coverage::{ + GroundPoint, Sensor, central_angle_to_station, elevation_angle, line_of_sight, local_up, + sub_satellite_point, visible_from, visible_from_geodetic, +}; +use apsis_geodesy::{Geodetic, WGS84}; +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; + +/// A station at the equator on the prime meridian, at sea level. +fn equatorial_station() -> GroundPoint { + GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 0.0) +} + +#[test] +fn a_station_on_the_equator_is_at_the_equatorial_radius() { + let station = equatorial_station(); + let position = station.ecef(); + assert!((position.x - BODY_RADIUS).abs() < 1e-6); + assert!(position.y.abs() < 1e-6); + assert!(position.z.abs() < 1e-6); + assert!((station.radius() - BODY_RADIUS).abs() < 1e-6); +} + +#[test] +fn a_station_at_the_pole_is_at_the_polar_radius() { + // The polar radius is the semi-minor axis of the ellipsoid, which is + // about 21 km less than the equatorial radius. + let station = GroundPoint::new(Geodetic::new(FRAC_PI_2, 0.0, 0.0), 0.0); + let radius = station.radius(); + assert!((radius - 6_356.752).abs() < 0.01, "radius {radius}"); + assert!(radius < BODY_RADIUS); +} + +#[test] +fn the_local_up_equals_the_position_direction_at_the_equator() { + let station = equatorial_station(); + let up = local_up(&station); + assert!((up.x - 1.0).abs() < 1e-12); + assert!(up.y.abs() < 1e-12); + assert!(up.z.abs() < 1e-12); +} + +#[test] +fn the_local_up_at_the_pole_is_the_third_axis() { + let station = GroundPoint::new(Geodetic::new(FRAC_PI_2, 0.0, 0.0), 0.0); + let up = local_up(&station); + assert!(up.x.abs() < 1e-12); + assert!(up.y.abs() < 1e-12); + assert!((up.z - 1.0).abs() < 1e-12); +} + +#[test] +fn a_satellite_overhead_is_at_the_zenith() { + let station = equatorial_station(); + // Directly above the station along the local normal. + let satellite = Vec3::new(BODY_RADIUS + 550.0, 0.0, 0.0); + let elevation = elevation_angle(satellite, &station).expect("valid"); + assert!( + (elevation - FRAC_PI_2).abs() < 1e-9, + "elevation {elevation}" + ); +} + +#[test] +fn a_satellite_at_the_horizon_has_zero_elevation() { + let station = equatorial_station(); + // The local horizon is where the line of sight is tangent to the + // sphere, at a central angle of about 22.98 degrees for a satellite at + // 550 km. At that geometry the elevation is zero and the satellite is + // at its greatest distance along the surface. + let central = (BODY_RADIUS / (BODY_RADIUS + 550.0)).acos(); + let satellite = Vec3::new( + (BODY_RADIUS + 550.0) * central.cos(), + (BODY_RADIUS + 550.0) * central.sin(), + 0.0, + ); + let elevation = elevation_angle(satellite, &station).expect("valid"); + assert!(elevation.abs() < 1e-9, "elevation {elevation}"); +} + +#[test] +fn a_satellite_below_the_horizon_has_negative_elevation() { + let station = equatorial_station(); + // On the far side of the Earth, the elevation is negative. + let satellite = Vec3::new(-(BODY_RADIUS + 550.0), 0.0, 0.0); + let elevation = elevation_angle(satellite, &station).expect("valid"); + assert!(elevation < -1.0, "elevation {elevation}"); +} + +#[test] +fn the_elevation_is_symmetric_about_the_overhead_direction() { + let station = equatorial_station(); + // Two satellites at the same angular distance from the overhead point, + // one to the east and one to the west, have the same elevation. + let angle = 0.3_f64; + let east = Vec3::new( + (BODY_RADIUS + 550.0) * angle.cos(), + (BODY_RADIUS + 550.0) * angle.sin(), + 0.0, + ); + let west = Vec3::new( + (BODY_RADIUS + 550.0) * angle.cos(), + -(BODY_RADIUS + 550.0) * angle.sin(), + 0.0, + ); + let east_elevation = elevation_angle(east, &station).expect("valid"); + let west_elevation = elevation_angle(west, &station).expect("valid"); + assert!((east_elevation - west_elevation).abs() < 1e-12); +} + +#[test] +fn visibility_respects_the_elevation_mask() { + let station = GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 0.2); + let sensor = Sensor::omnidirectional(station); + + // Overhead is well above the mask. + assert_eq!( + visible_from(Vec3::new(BODY_RADIUS + 550.0, 0.0, 0.0), &sensor), + Some(true) + ); + // On the horizon the satellite is below a 0.2 radian mask, because the + // horizon direction is a right angle from the local up vector. + assert_eq!( + visible_from(Vec3::new(0.0, BODY_RADIUS + 550.0, 0.0), &sensor), + Some(false) + ); + // On the far side is well below it. + assert_eq!( + visible_from(Vec3::new(-(BODY_RADIUS + 550.0), 0.0, 0.0), &sensor), + Some(false) + ); +} + +#[test] +fn a_zero_mask_admits_a_satellite_at_the_horizon() { + let sensor = Sensor::omnidirectional(equatorial_station()); + // A satellite just inside the horizon direction has a small positive + // elevation and is visible with a zero mask. + let central = (BODY_RADIUS / (BODY_RADIUS + 550.0)).acos() * 0.999; + let satellite = Vec3::new( + (BODY_RADIUS + 550.0) * central.cos(), + (BODY_RADIUS + 550.0) * central.sin(), + 0.0, + ); + assert_eq!(visible_from(satellite, &sensor), Some(true)); +} + +#[test] +fn a_directional_sensor_only_sees_its_forward_hemisphere() { + // The station points at the zenith, so it sees a satellite overhead but + // not one on the horizon behind it. + let station = equatorial_station(); + let sensor = Sensor::directional(station, Vec3::new(1.0, 0.0, 0.0)); + + assert_eq!( + visible_from(Vec3::new(BODY_RADIUS + 550.0, 0.0, 0.0), &sensor), + Some(true) + ); + // A satellite to the east is at a right angle to the pointing, so it is + // not in the forward hemisphere. + assert_eq!( + visible_from(Vec3::new(0.0, BODY_RADIUS + 550.0, 0.0), &sensor), + Some(false) + ); +} + +#[test] +fn a_directional_sensor_normalizes_its_pointing() { + // An unnormalized pointing vector must behave the same as a unit one. + let station = equatorial_station(); + let unit = Sensor::directional(station, Vec3::new(1.0, 0.0, 0.0)); + let scaled = Sensor::directional(station, Vec3::new(100.0, 0.0, 0.0)); + let satellite = Vec3::new(BODY_RADIUS + 550.0, 0.0, 0.0); + assert_eq!( + visible_from(satellite, &unit), + visible_from(satellite, &scaled) + ); +} + +#[test] +fn a_degenerate_pointing_is_rejected() { + let station = equatorial_station(); + let sensor = Sensor::directional(station, Vec3::ZERO); + let satellite = Vec3::new(BODY_RADIUS + 550.0, 0.0, 0.0); + assert_eq!(visible_from(satellite, &sensor), None); +} + +#[test] +fn the_geodetic_helper_agrees_with_the_sensor_path() { + let position = Geodetic::new(0.9, -1.2, 100.0); + let mask = 0.15; + let satellite = { + // A point above the station. + let station = GroundPoint::new(position, mask); + station.ecef() * 1.09 + }; + let via_sensor = visible_from_geodetic(satellite, position, mask).expect("valid"); + let via_sensor_type = visible_from( + satellite, + &Sensor::omnidirectional(GroundPoint::new(position, mask)), + ) + .expect("valid"); + assert_eq!(via_sensor, via_sensor_type); + assert!(via_sensor); +} + +#[test] +fn the_sub_satellite_point_is_at_zero_height() { + let satellite = Vec3::new(BODY_RADIUS + 550.0, 0.0, 0.0); + let point = sub_satellite_point(satellite, WGS84).expect("valid"); + assert_eq!(point.position.height, 0.0); + assert!(point.position.latitude.abs() < 1e-9); + assert!(point.position.longitude.abs() < 1e-9); +} + +#[test] +fn the_sub_satellite_point_keeps_the_longitude() { + let satellite = Vec3::new(0.0, 4_000.0, 4_500.0); + let point = sub_satellite_point(satellite, WGS84).expect("valid"); + assert!((point.position.longitude - FRAC_PI_2).abs() < 1e-9); + assert!(point.position.latitude > 0.0); +} + +#[test] +fn a_satellite_on_the_polar_axis_has_no_sub_satellite_point() { + let satellite = Vec3::new(0.0, 0.0, BODY_RADIUS + 550.0); + assert_eq!(sub_satellite_point(satellite, WGS84), None); + assert_eq!(sub_satellite_point(Vec3::splat(f64::NAN), WGS84), None); +} + +#[test] +fn the_central_angle_to_a_station_follows_the_spherical_relation() { + // A satellite at the zenith of the station has a central angle of zero. + let station = equatorial_station(); + let satellite = Vec3::new(BODY_RADIUS + 550.0, 0.0, 0.0); + let central = central_angle_to_station(satellite, &station, FRAC_PI_2).expect("valid"); + assert!(central.abs() < 1e-9); + + // At a ten degree elevation the central angle is the footprint angle. + let central = + central_angle_to_station(satellite, &station, 10.0_f64.to_radians()).expect("valid"); + assert!((central.to_degrees() - 14.956_635_124_957_29).abs() < 1e-6); +} + +#[test] +fn invalid_inputs_are_rejected() { + let station = equatorial_station(); + assert_eq!(elevation_angle(Vec3::splat(f64::NAN), &station), None); + // A station with an out of range mask is not usable. + let bad = GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), PI); + assert!(!bad.is_valid()); + assert_eq!( + visible_from( + Vec3::new(BODY_RADIUS + 550.0, 0.0, 0.0), + &Sensor::omnidirectional(bad) + ), + None + ); + // A satellite at the same place as the station has no line of sight. + let coincident = station.ecef(); + assert_eq!(line_of_sight(coincident, &station), None); +} + +#[test] +fn a_station_above_sea_level_has_a_larger_radius() { + let sea_level = GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 0.0); + let high = GroundPoint::new(Geodetic::new(0.0, 0.0, 3_000.0), 0.0); + assert!((high.radius() - sea_level.radius() - 3.0).abs() < 1e-6); +} + +#[test] +fn the_sensor_and_ground_point_validity_is_reported() { + assert!(equatorial_station().is_valid()); + assert!(GroundPoint::new(Geodetic::new(0.9, 0.9, 0.0), 0.1).is_valid()); + assert!(!GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 2.0).is_valid()); + assert!(!GroundPoint::new(Geodetic::new(0.0, 0.0, f64::NAN), 0.1).is_valid()); +}