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

network: add coverage geometry, footprint, and visibility

This commit is contained in:
ApsisRS contributors 2026-09-26 15:48:14 +08:00
parent a3741cafae
commit c75fb8bc5c
9 changed files with 1289 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Cargo.lock generated
View File

@ -53,6 +53,11 @@ version = "0.1.0"
[[package]] [[package]]
name = "apsis-coverage" name = "apsis-coverage"
version = "0.1.0" version = "0.1.0"
dependencies = [
"apsis-core",
"apsis-geodesy",
"apsis-math",
]
[[package]] [[package]]
name = "apsis-doppler" name = "apsis-doppler"

3
crates/network/apsis-coverage/Cargo.toml
View File

@ -11,3 +11,6 @@ publish = false
workspace = true workspace = true
[dependencies] [dependencies]
apsis-core.workspace = true
apsis-math.workspace = true
apsis-geodesy.workspace = true

176
crates/network/apsis-coverage/src/footprint.rs Normal file
View File

@ -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<f64> {
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<Footprint> {
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<f64> {
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<f64> {
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<f64> {
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<bool> {
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<f64> {
elevation_from_central_angle(
satellite_radius,
footprint.body_radius,
footprint.central_angle,
)
}

135
crates/network/apsis-coverage/src/geometry.rs Normal file
View File

@ -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<f64> {
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<f64> {
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<f64> {
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<f64> {
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())
}

64
crates/network/apsis-coverage/src/lib.rs
View File

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

245
crates/network/apsis-coverage/src/visibility.rs Normal file
View File

@ -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<Vec3>,
}
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<SubSatellitePoint> {
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<f64> {
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<f64> {
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<bool> {
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<bool> {
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<f64> {
if !satellite_ecef_km.is_finite() || !ground.is_valid() {
return None;
}
earth_central_angle(satellite_ecef_km.norm(), ground.radius(), elevation)
}

180
crates/network/apsis-coverage/tests/footprint.rs Normal file
View File

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

200
crates/network/apsis-coverage/tests/geometry.rs Normal file
View File

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

282
crates/network/apsis-coverage/tests/visibility.rs Normal file
View File

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