network: add inter-satellite link geometry

Cargo.lock

@@ -112,6 +112,10 @@ dependencies = [
 [[package]]
 name = "apsis-isl"
 version = "0.1.0"
+dependencies = [
+ "apsis-core",
+ "apsis-math",
+]
 
 [[package]]
 name = "apsis-kepler"

crates/network/apsis-isl/Cargo.toml

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

crates/network/apsis-isl/src/geometry.rs

@@ -0,0 +1,114 @@
+//! Link geometry and the Earth blocking test.
+
+use crate::link::LinkGeometry;
+use apsis_math::Vec3;
+
+/// Returns the closest distance from the origin to the line segment between
+/// two points.
+///
+/// This is the perpendicular distance when the closest point lies between
+/// the two ends, and the distance to the nearer end otherwise. For an
+/// inter-satellite link the first case is the one that matters, because the
+/// blocking body lies between the terminals.
+///
+/// Returns `None` when either point is not finite.
+pub fn closest_approach_to_centre(first: Vec3, second: Vec3) -> Option<f64> {
+    if !first.is_finite() || !second.is_finite() {
+        return None;
+    }
+    let segment = second - first;
+    let length_squared = segment.norm_squared();
+    if length_squared == 0.0 {
+        // The two points coincide, so the distance is the radius of that
+        // point.
+        return Some(first.norm());
+    }
+
+    // The parameter of the closest point on the infinite line, clamped to
+    // the segment.
+    let parameter = (-first.dot(segment) / length_squared).clamp(0.0, 1.0);
+    let closest = first + segment * parameter;
+    Some(closest.norm())
+}
+
+/// Returns `true` when the two satellites can see each other past the body.
+///
+/// The test is whether the line segment between the two positions stays
+/// outside a sphere of the given radius. A link that grazes the surface
+/// exactly is treated as clear, because the limit is the geometric horizon.
+///
+/// Returns `None` when either position is not finite or the body radius is
+/// not positive.
+pub fn mutual_visibility(first: Vec3, second: Vec3, body_radius: f64) -> Option<bool> {
+    if !body_radius.is_finite() || body_radius <= 0.0 {
+        return None;
+    }
+    let closest = closest_approach_to_centre(first, second)?;
+    Some(closest >= body_radius)
+}
+
+/// Returns the full geometry of a link between two satellites.
+///
+/// The directions are unit vectors from each satellite toward the other. The
+/// closest approach to the centre is included so that the caller can judge
+/// how much margin the link has.
+///
+/// Returns `None` when either position is not finite, when the two coincide,
+/// or when the body radius is not usable. A cell that reports its own
+/// position for both ends has no link, so that case is rejected rather than
+/// returned with a zero range.
+pub fn link_geometry(first: Vec3, second: Vec3, body_radius: f64) -> Option<LinkGeometry> {
+    if !body_radius.is_finite() || body_radius < 0.0 {
+        return None;
+    }
+    if !first.is_finite() || !second.is_finite() {
+        return None;
+    }
+    let separation = second - first;
+    let range = separation.norm();
+    if range == 0.0 {
+        return None;
+    }
+    Some(LinkGeometry {
+        range_km: range,
+        direction_first_to_second: separation / range,
+        direction_second_to_first: -separation / range,
+        centre_distance_km: closest_approach_to_centre(first, second)?,
+    })
+}
+
+/// Returns `true` when the Sun is within the field of view of an optical
+/// link.
+///
+/// The test is the angle between the link direction and the direction to the
+/// Sun. If that angle is within the field of view at either terminal, the
+/// Sun is in the beam and the link cannot be used.
+///
+/// The field of view is given as a half angle in radians, which is the
+/// acceptance cone of the telescope. A typical optical terminal has a half
+/// angle of a few microradians.
+///
+/// The direction to the Sun is a unit vector in the same frame as the
+/// positions. This crate does not compute it, because that needs an
+/// ephemeris of the Sun, which the orbit layer does not currently provide; a
+/// caller with such an ephemeris supplies the direction.
+///
+/// Returns `None` when the directions or the field of view are not usable.
+pub fn sun_blocks_link(
+    link: &LinkGeometry,
+    sun_direction: Vec3,
+    field_of_view: f64,
+) -> Option<bool> {
+    if !sun_direction.is_finite() || !field_of_view.is_finite() {
+        return None;
+    }
+    if field_of_view < 0.0 {
+        return None;
+    }
+    let sun = sun_direction.normalize()?;
+    let forward = link.direction_first_to_second.normalize()?;
+    let backward = link.direction_second_to_first.normalize()?;
+    let forward_angle = forward.dot(sun).clamp(-1.0, 1.0).acos();
+    let backward_angle = backward.dot(sun).clamp(-1.0, 1.0).acos();
+    Some(forward_angle <= field_of_view || backward_angle <= field_of_view)
+}

crates/network/apsis-isl/src/lib.rs

@@ -1 +1,45 @@
-//! Inter-satellite link modeling.
+//! Inter-satellite link geometry and visibility.
+//!
+//! Two satellites can form a link when the straight line between them does
+//! not pass through the Earth. This crate decides that condition, gives the
+//! geometry a link budget needs, and reports the range and the angles at
+//! which each satellite must point.
+//!
+//! The blocking test is whether the line segment between the two positions
+//! has any point inside the body. For a spherical body that is a
+//! perpendicular distance test, which is exact and cheap, so it is the one
+//! used here.
+//!
+//! An optical link adds a second condition: the Sun must not be within the
+//! field of view of either terminal, because the background would swamp the
+//! signal. That is a separate test, offered as [`sun_blocks_link`], which
+//! takes the direction to the Sun.
+//!
+//! Lengths are in kilometers and angles in radians. The positions are in a
+//! frame in which the Earth is a sphere of the given radius, which for link
+//! geometry is the same frame the orbit layer produces.
+//!
+//! ```
+//! use apsis_isl::{LinkGeometry, link_geometry, mutual_visibility};
+//! use apsis_math::Vec3;
+//!
+//! // Two satellites on the same side of the Earth see each other.
+//! let first = Vec3::new(6_928.137, 0.0, 0.0);
+//! let second = Vec3::new(6_868.137, 2_500.0, 0.0);
+//! assert_eq!(mutual_visibility(first, second, 6_378.137), Some(true));
+//!
+//! // Two satellites on opposite sides do not.
+//! let opposite = Vec3::new(-6_928.137, 0.0, 0.0);
+//! assert_eq!(mutual_visibility(first, opposite, 6_378.137), Some(false));
+//! ```
+
+#![forbid(unsafe_code)]
+#![deny(missing_docs)]
+
+mod geometry;
+mod link;
+
+pub use geometry::{closest_approach_to_centre, link_geometry, mutual_visibility, sun_blocks_link};
+pub use link::{LinkGeometry, LinkState};
+
+pub use apsis_math::Vec3;

crates/network/apsis-isl/src/link.rs

@@ -0,0 +1,79 @@
+//! The link geometry type and the state of a link.
+
+use apsis_math::Vec3;
+
+/// The geometry of a link between two satellites.
+///
+/// The range is the straight line distance, which a link budget uses for the
+/// free space path loss. The two pointing vectors are the directions each
+/// terminal must aim, expressed in the frame the positions were given in.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct LinkGeometry {
+    /// The distance between the two satellites, in kilometers.
+    pub range_km: f64,
+    /// The unit vector from the first satellite toward the second.
+    pub direction_first_to_second: Vec3,
+    /// The unit vector from the second satellite toward the first.
+    pub direction_second_to_first: Vec3,
+    /// The closest the line between the two passes to the centre of the
+    /// body, in kilometers.
+    ///
+    /// This is the quantity that decides whether the link is blocked. It is
+    /// reported because a link that clears the limb by a small margin is
+    /// fragile: atmospheric refraction and the oblateness of the body both
+    /// eat into the margin.
+    pub centre_distance_km: f64,
+}
+
+impl LinkGeometry {
+    /// Returns the grazing margin of the link, in kilometers.
+    ///
+    /// This is how far the line between the two satellites clears the
+    /// surface, measured at the closest point. A negative value means the
+    /// link is blocked.
+    pub fn surface_clearance(&self, body_radius: f64) -> f64 {
+        self.centre_distance_km - body_radius
+    }
+
+    /// Returns `true` when the link clears the body.
+    pub fn is_clear(&self, body_radius: f64) -> bool {
+        self.centre_distance_km >= body_radius
+    }
+}
+
+/// The state of a link at one instant.
+///
+/// A link is either carrying traffic or not. The reasons a link can be
+/// unavailable are kept separate, because routing across a constellation has
+/// to distinguish a link that is geometrically blocked from one that is
+/// merely not established yet.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
+pub enum LinkState {
+    /// The link is geometrically clear and may be used.
+    Available,
+    /// The Earth blocks the line between the two satellites.
+    Blocked,
+    /// The Sun is within the field of view of one of the terminals, so an
+    /// optical link cannot be used.
+    SunInField,
+    /// The link is clear but is not established, for example because it has
+    /// not yet been acquired or its capacity is exhausted.
+    NotEstablished,
+}
+
+impl LinkState {
+    /// Returns `true` when the link may carry traffic.
+    pub const fn is_usable(self) -> bool {
+        matches!(self, Self::Available)
+    }
+
+    /// Returns a short name for the state.
+    pub const fn name(self) -> &'static str {
+        match self {
+            Self::Available => "available",
+            Self::Blocked => "blocked",
+            Self::SunInField => "sun in field",
+            Self::NotEstablished => "not established",
+        }
+    }
+}

crates/network/apsis-isl/tests/geometry.rs

@@ -0,0 +1,298 @@
+//! Tests for inter-satellite link geometry.
+//!
+//! For two satellites at the same radius separated by a central angle
+//! `theta`, the closest approach of the chord to the centre is
+//! `Rs cos(theta / 2)`. That closed form is used as the reference here.
+
+use apsis_isl::{
+    LinkState, closest_approach_to_centre, link_geometry, mutual_visibility, sun_blocks_link,
+};
+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 550 km circular orbit radius, in kilometers.
+const ORBIT_RADIUS: f64 = BODY_RADIUS + 550.0;
+
+/// Returns a satellite on the equator at the given central angle from the
+/// first axis.
+fn satellite_at(central_angle: f64) -> Vec3 {
+    Vec3::new(
+        ORBIT_RADIUS * central_angle.cos(),
+        ORBIT_RADIUS * central_angle.sin(),
+        0.0,
+    )
+}
+
+#[test]
+fn the_closest_approach_matches_the_chord_geometry() {
+    // Two satellites at the same radius separated by a central angle have a
+    // closest approach of Rs cos(theta / 2).
+    for degrees in [0.0_f64, 10.0, 30.0, 60.0, 120.0, 180.0] {
+        let central = degrees.to_radians();
+        let first = satellite_at(0.0);
+        let second = satellite_at(central);
+        let closest = closest_approach_to_centre(first, second).expect("valid positions");
+        let expected = ORBIT_RADIUS * (central / 2.0).cos();
+        assert!(
+            (closest - expected).abs() < 1e-6,
+            "central {degrees}: {closest} against {expected}"
+        );
+    }
+}
+
+#[test]
+fn the_closest_approach_is_the_distance_for_coincident_points() {
+    let point = Vec3::new(3.0, 4.0, 0.0);
+    let closest = closest_approach_to_centre(point, point).expect("valid");
+    assert!((closest - 5.0).abs() < 1e-12);
+}
+
+#[test]
+fn the_closest_approach_is_the_nearer_end_when_the_perpendicular_falls_outside() {
+    // Two points on the same side of the origin, so the foot of the
+    // perpendicular lies outside the segment.
+    let first = Vec3::new(10.0, 0.0, 0.0);
+    let second = Vec3::new(20.0, 0.0, 0.0);
+    let closest = closest_approach_to_centre(first, second).expect("valid");
+    assert!((closest - 10.0).abs() < 1e-12);
+}
+
+#[test]
+fn satellites_on_the_same_side_see_each_other() {
+    for degrees in [10.0_f64, 20.0, 30.0] {
+        let first = satellite_at(0.0);
+        let second = satellite_at(degrees.to_radians());
+        assert_eq!(
+            mutual_visibility(first, second, BODY_RADIUS),
+            Some(true),
+            "central {degrees}"
+        );
+    }
+}
+
+#[test]
+fn satellites_on_opposite_sides_do_not() {
+    for degrees in [60.0_f64, 120.0, 180.0] {
+        let first = satellite_at(0.0);
+        let second = satellite_at(degrees.to_radians());
+        assert_eq!(
+            mutual_visibility(first, second, BODY_RADIUS),
+            Some(false),
+            "central {degrees}"
+        );
+    }
+}
+
+#[test]
+fn the_visibility_threshold_is_where_the_chord_grazes_the_surface() {
+    // The chord has a closest approach equal to the body radius when
+    // cos(theta / 2) = Re / Rs.
+    let threshold = 2.0 * (BODY_RADIUS / ORBIT_RADIUS).acos();
+    let first = satellite_at(0.0);
+
+    // Just inside the threshold the link is clear.
+    let inside = satellite_at(threshold - 1e-6);
+    assert_eq!(mutual_visibility(first, inside, BODY_RADIUS), Some(true));
+
+    // Just outside it is blocked.
+    let outside = satellite_at(threshold + 1e-6);
+    assert_eq!(mutual_visibility(first, outside, BODY_RADIUS), Some(false));
+
+    // Exactly at the threshold the grazing link counts as clear.
+    let grazing = satellite_at(threshold);
+    assert_eq!(mutual_visibility(first, grazing, BODY_RADIUS), Some(true));
+}
+
+#[test]
+fn a_higher_orbit_allows_a_wider_link() {
+    // Two satellites at the same radius can see each other while the central
+    // angle between them is less than 2 acos(Re / Rs). For a 550 km orbit
+    // that limit is 45.97 degrees, and for a 20 200 km orbit it is 152.23
+    // degrees. A separation of 90 degrees therefore falls between the two.
+    let central = FRAC_PI_2;
+    let low_radius = BODY_RADIUS + 550.0;
+    let high_radius = BODY_RADIUS + 20_200.0;
+
+    let low_first = Vec3::new(low_radius, 0.0, 0.0);
+    let low_second = Vec3::new(low_radius * central.cos(), low_radius * central.sin(), 0.0);
+    assert_eq!(
+        mutual_visibility(low_first, low_second, BODY_RADIUS),
+        Some(false)
+    );
+
+    let high_first = Vec3::new(high_radius, 0.0, 0.0);
+    let high_second = Vec3::new(
+        high_radius * central.cos(),
+        high_radius * central.sin(),
+        0.0,
+    );
+    assert_eq!(
+        mutual_visibility(high_first, high_second, BODY_RADIUS),
+        Some(true)
+    );
+}
+
+#[test]
+fn the_link_geometry_reports_range_and_directions() {
+    let first = Vec3::new(7_000.0, 0.0, 0.0);
+    let second = Vec3::new(7_000.0, 3_000.0, 0.0);
+    let geometry = link_geometry(first, second, BODY_RADIUS).expect("valid");
+
+    // The range is the distance between the positions.
+    assert!((geometry.range_km - first.distance(second)).abs() < 1e-12);
+    // The two directions are opposite unit vectors.
+    assert!((geometry.direction_first_to_second.norm() - 1.0).abs() < 1e-12);
+    assert!((geometry.direction_second_to_first.norm() - 1.0).abs() < 1e-12);
+    assert!(
+        (geometry.direction_first_to_second + geometry.direction_second_to_first).norm() < 1e-12
+    );
+    // The direction points from the first to the second.
+    let expected = (second - first).normalize().expect("non-zero separation");
+    assert!((geometry.direction_first_to_second - expected).norm() < 1e-12);
+}
+
+#[test]
+fn the_link_geometry_reports_the_surface_clearance() {
+    let first = satellite_at(0.0);
+    let second = satellite_at(20.0_f64.to_radians());
+    let geometry = link_geometry(first, second, BODY_RADIUS).expect("valid");
+    let clearance = geometry.surface_clearance(BODY_RADIUS);
+    assert!(clearance > 0.0, "clearance {clearance}");
+    assert!(geometry.is_clear(BODY_RADIUS));
+
+    let blocked = link_geometry(first, satellite_at(PI), BODY_RADIUS).expect("valid");
+    assert!(blocked.surface_clearance(BODY_RADIUS) < 0.0);
+    assert!(!blocked.is_clear(BODY_RADIUS));
+}
+
+#[test]
+fn coincident_satellites_have_no_link() {
+    let point = Vec3::new(7_000.0, 0.0, 0.0);
+    assert_eq!(link_geometry(point, point, BODY_RADIUS), None);
+}
+
+#[test]
+fn invalid_link_arguments_are_rejected() {
+    let first = satellite_at(0.0);
+    let second = satellite_at(0.1);
+    assert_eq!(mutual_visibility(first, second, 0.0), None);
+    assert_eq!(mutual_visibility(first, second, -1.0), None);
+    assert_eq!(
+        mutual_visibility(Vec3::splat(f64::NAN), second, BODY_RADIUS),
+        None
+    );
+    assert_eq!(
+        mutual_visibility(first, Vec3::splat(f64::INFINITY), BODY_RADIUS),
+        None
+    );
+    assert_eq!(link_geometry(first, second, -1.0), None);
+    assert_eq!(
+        closest_approach_to_centre(Vec3::splat(f64::NAN), second),
+        None
+    );
+}
+
+#[test]
+fn a_link_along_a_radius_has_the_altitude_as_its_closest_approach() {
+    // Two satellites on the same radial line, one above the other. The
+    // lower one is the closest point of the segment.
+    let low = Vec3::new(7_000.0, 0.0, 0.0);
+    let high = Vec3::new(8_000.0, 0.0, 0.0);
+    let closest = closest_approach_to_centre(low, high).expect("valid");
+    assert!((closest - 7_000.0).abs() < 1e-12);
+}
+
+#[test]
+fn the_sun_blocks_a_link_when_it_lies_along_it() {
+    let first = Vec3::new(7_000.0, 0.0, 0.0);
+    let second = Vec3::new(7_000.0, 1_000.0, 0.0);
+    let geometry = link_geometry(first, second, BODY_RADIUS).expect("valid");
+    let field_of_view = 1e-5;
+
+    // The Sun in the same direction as the link, seen from the first
+    // terminal.
+    let along = geometry.direction_first_to_second;
+    assert_eq!(sun_blocks_link(&geometry, along, field_of_view), Some(true));
+
+    // The Sun behind the link, seen from the second terminal.
+    let behind = geometry.direction_second_to_first;
+    assert_eq!(
+        sun_blocks_link(&geometry, behind, field_of_view),
+        Some(true)
+    );
+
+    // The Sun at right angles to the link.
+    let transverse = Vec3::new(0.0, 0.0, 1.0);
+    assert_eq!(
+        sun_blocks_link(&geometry, transverse, field_of_view),
+        Some(false)
+    );
+}
+
+#[test]
+fn a_zero_field_of_view_blocks_only_an_exact_alignment() {
+    let geometry = link_geometry(
+        Vec3::new(7_000.0, 0.0, 0.0),
+        Vec3::new(7_000.0, 1_000.0, 0.0),
+        BODY_RADIUS,
+    )
+    .expect("valid");
+    // Exactly along the link, the angle is zero, which is within a zero
+    // field of view.
+    let along = geometry.direction_first_to_second;
+    assert_eq!(sun_blocks_link(&geometry, along, 0.0), Some(true));
+    // A tiny offset is outside a zero field of view.
+    let offset = (along + Vec3::new(0.0, 0.0, 0.01))
+        .normalize()
+        .expect("valid");
+    assert_eq!(sun_blocks_link(&geometry, offset, 0.0), Some(false));
+}
+
+#[test]
+fn the_sun_test_rejects_invalid_arguments() {
+    let geometry = link_geometry(
+        Vec3::new(7_000.0, 0.0, 0.0),
+        Vec3::new(7_000.0, 1_000.0, 0.0),
+        BODY_RADIUS,
+    )
+    .expect("valid");
+    assert_eq!(
+        sun_blocks_link(&geometry, Vec3::splat(f64::NAN), 1e-5),
+        None
+    );
+    assert_eq!(
+        sun_blocks_link(&geometry, Vec3::new(0.0, 0.0, 1.0), -1.0),
+        None
+    );
+    assert_eq!(sun_blocks_link(&geometry, Vec3::ZERO, 1e-5), None);
+}
+
+#[test]
+fn the_link_state_reports_usability() {
+    assert!(LinkState::Available.is_usable());
+    assert!(!LinkState::Blocked.is_usable());
+    assert!(!LinkState::SunInField.is_usable());
+    assert!(!LinkState::NotEstablished.is_usable());
+
+    assert_eq!(LinkState::Available.name(), "available");
+    assert_eq!(LinkState::Blocked.name(), "blocked");
+    assert_eq!(LinkState::SunInField.name(), "sun in field");
+    assert_eq!(LinkState::NotEstablished.name(), "not established");
+}
+
+#[test]
+fn a_link_between_neighbouring_planes_is_clear() {
+    // Two satellites in adjacent planes of a Walker pattern at 550 km are
+    // separated by about 120 degrees of right ascension at the equator,
+    // which at the same radius would be blocked. Satellites that are near
+    // each other in time are much closer, which is what makes cross-plane
+    // links viable only between nearby satellites.
+    let first = satellite_at(0.0);
+    let neighbour = satellite_at(5.0_f64.to_radians());
+    assert_eq!(mutual_visibility(first, neighbour, BODY_RADIUS), Some(true));
+    let far = satellite_at(FRAC_PI_2);
+    assert_eq!(mutual_visibility(first, far, BODY_RADIUS), Some(false));
+}