rf: add Doppler shift and range rate models

Cargo.lock

@@ -47,6 +47,10 @@ version = "0.1.0"
 [[package]]
 name = "apsis-doppler"
 version = "0.1.0"
+dependencies = [
+ "apsis-core",
+ "apsis-math",
+]
 
 [[package]]
 name = "apsis-ephemeris"

crates/rf/apsis-doppler/Cargo.toml

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

crates/rf/apsis-doppler/src/lib.rs

@@ -1 +1,46 @@
 //! Doppler shift computation.
+//!
+//! A satellite and a ground station move relative to each other, so the
+//! frequency that arrives at the receiver differs from the frequency that
+//! left the transmitter. This crate computes that shift from the relative
+//! velocity.
+//!
+//! Two effects are present and they are separate. The classical shift is
+//! proportional to the range rate divided by the speed of light, and it is
+//! what a Doppler tracker compensates. The relativistic time dilation
+//! between the two clocks is a second order correction that is far smaller
+//! for orbital velocities, and it is offered separately so that a caller can
+//! include it when the application needs it.
+//!
+//! The sign convention used here is that a positive range rate, meaning the
+//! two terminals are moving apart, produces a negative shift and therefore a
+//! lower received frequency. This matches the everyday observation that a
+//! receding source sounds lower.
+//!
+//! ```
+//! use apsis_doppler::{DopplerShift, range_rate_from_velocity};
+//! use apsis_math::Vec3;
+//!
+//! // A satellite at the first axis moving away along it.
+//! let range_rate = range_rate_from_velocity(
+//!     Vec3::new(7_000_000.0, 0.0, 0.0),
+//!     Vec3::new(1_000.0, 0.0, 0.0),
+//! );
+//! assert_eq!(range_rate, Some(1_000.0));
+//!
+//! let shift = DopplerShift::from_range_rate(12.0e9, 1_000.0).expect("valid");
+//! assert!(shift.shifted_frequency() < 12.0e9);
+//! ```
+
+#![forbid(unsafe_code)]
+#![deny(missing_docs)]
+
+mod shift;
+
+pub use shift::{
+    DopplerShift, acquisition_span, angle_from_range_rate, classical_shift, maximum_shift,
+    range_rate_from_velocity, relativistic_doppler_factor,
+};
+
+pub use apsis_core::constant::SPEED_OF_LIGHT;
+pub use apsis_math::Vec3;

crates/rf/apsis-doppler/src/shift.rs

@@ -0,0 +1,251 @@
+//! The Doppler shift and the range rate that causes it.
+
+use apsis_core::constant::SPEED_OF_LIGHT;
+use apsis_math::Vec3;
+use core::fmt;
+
+/// Returns the range rate between two terminals, in meters per second.
+///
+/// The range rate is the component of the relative velocity along the line
+/// between the terminals. A positive value means the separation is
+/// increasing.
+///
+/// The first argument is the position of the remote terminal relative to the
+/// local one, so that its magnitude is the separation. The second is the
+/// velocity of the remote terminal relative to the local one. Both are in
+/// SI units.
+///
+/// Returns `None` when the separation is zero, because no line of sight is
+/// defined there, or when either argument is not finite.
+pub fn range_rate_from_velocity(relative_position: Vec3, relative_velocity: Vec3) -> Option<f64> {
+    if !relative_position.is_finite() || !relative_velocity.is_finite() {
+        return None;
+    }
+    let range = relative_position.norm();
+    if range == 0.0 {
+        return None;
+    }
+    Some(relative_position.dot(relative_velocity) / range)
+}
+
+/// Returns the classical Doppler shift of a frequency for a range rate, in
+/// hertz.
+///
+/// The relation is
+///
+/// ```text
+/// delta f = -f_0 (v_r / c)
+/// ```
+///
+/// where `v_r` is the range rate, positive when receding. The shift is
+/// negative for a receding terminal, so the received frequency is lower than
+/// the transmitted one.
+///
+/// The magnitude of the shift is small compared to the carrier: a range rate
+/// of 7 km/s at 12 GHz gives about 280 kHz, which is two parts in a hundred
+/// thousand. The caller must therefore treat the frequency as a double, and
+/// the shifted frequency is returned separately rather than by adding a
+/// small shift to a large number.
+///
+/// Returns `None` when the frequency is not positive and finite, or when the
+/// range rate is not finite.
+pub fn classical_shift(carrier_frequency: f64, range_rate: f64) -> Option<f64> {
+    if !carrier_frequency.is_finite() || carrier_frequency <= 0.0 {
+        return None;
+    }
+    if !range_rate.is_finite() {
+        return None;
+    }
+    Some(-carrier_frequency * range_rate / SPEED_OF_LIGHT)
+}
+
+/// Returns the relativistic Doppler factor for a range rate.
+///
+/// The factor is `sqrt((1 - beta) / (1 + beta))` with `beta = v_r / c`.
+/// Multiplying the transmitted frequency by it gives the received frequency
+/// to all orders in `beta`, so it is not a small correction to the classical
+/// shift but the quantity the classical shift approximates.
+///
+/// Its expansion is `1 - beta + beta^2 / 2 - ...`. The first two terms are
+/// the classical result, which is `1 - beta` to within a factor that is
+/// itself second order, and the remaining terms are the difference between
+/// the two treatments. For a satellite at orbital velocity, `beta` is about
+/// `2.3e-5` and the second order term is about `2.7e-10` of the carrier, so
+/// the two agree to a few parts in ten to the tenth.
+///
+/// Returns `None` when the range rate is not finite, or when its magnitude
+/// reaches the speed of light, where the expression has no real value.
+pub fn relativistic_doppler_factor(range_rate: f64) -> Option<f64> {
+    if !range_rate.is_finite() {
+        return None;
+    }
+    let beta = range_rate / SPEED_OF_LIGHT;
+    if beta.abs() >= 1.0 {
+        return None;
+    }
+    Some(((1.0 - beta) / (1.0 + beta)).sqrt())
+}
+
+/// A Doppler shift applied to a carrier.
+///
+/// The type carries the carrier frequency and the range rate, so the
+/// transmitted and received frequencies and the shift itself are all
+/// available without ambiguity about which sign convention was used.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct DopplerShift {
+    /// The transmitted frequency, in hertz.
+    carrier: f64,
+    /// The range rate, in meters per second, positive when receding.
+    range_rate: f64,
+}
+
+impl DopplerShift {
+    /// Creates a shift from a carrier frequency and a range rate.
+    ///
+    /// Returns `None` when the frequency is not positive and finite, or when
+    /// the range rate is not finite.
+    pub fn from_range_rate(carrier_frequency: f64, range_rate: f64) -> Option<Self> {
+        classical_shift(carrier_frequency, range_rate)?;
+        Some(Self {
+            carrier: carrier_frequency,
+            range_rate,
+        })
+    }
+
+    /// Returns the transmitted frequency, in hertz.
+    pub const fn carrier_frequency(self) -> f64 {
+        self.carrier
+    }
+
+    /// Returns the range rate, in meters per second, positive when receding.
+    pub const fn range_rate(self) -> f64 {
+        self.range_rate
+    }
+
+    /// Returns the shift, in hertz.
+    ///
+    /// The value is negative when the terminals are receding.
+    pub fn shift(self) -> f64 {
+        classical_shift(self.carrier, self.range_rate).unwrap_or(0.0)
+    }
+
+    /// Returns the received frequency, in hertz.
+    pub fn shifted_frequency(self) -> f64 {
+        self.carrier + self.shift()
+    }
+
+    /// Returns the shift including the relativistic clock correction, in
+    /// hertz.
+    ///
+    /// The received frequency is the carrier multiplied by the relativistic
+    /// factor, so the shift is `f_0 (factor - 1)`.
+    ///
+    /// This is not the classical shift plus a small correction: the factor
+    /// already contains the first order term, and its expansion is
+    /// `1 - beta + beta^2 / 2 - ...`, whose leading term is exactly the
+    /// classical shift. The classical value is the first order
+    /// approximation to this one, and the two differ only by the second and
+    /// higher order terms, of order `beta^2`, which for a satellite is a few
+    /// parts in ten to the tenth of the carrier.
+    ///
+    /// The factor is evaluated first and the carrier multiplied afterwards,
+    /// because the difference from one is tiny and would lose precision if
+    /// the shift were added to the carrier before the correction.
+    ///
+    /// Returns `None` when the range rate reaches the speed of light.
+    pub fn shift_including_relativity(self) -> Option<f64> {
+        let factor = relativistic_doppler_factor(self.range_rate)?;
+        Some(self.carrier * (factor - 1.0))
+    }
+
+    /// Returns the shift in parts per million of the carrier.
+    ///
+    /// This is a convenient way to compare the shift of links on different
+    /// carriers.
+    pub fn shift_in_ppm(self) -> f64 {
+        self.shift() / self.carrier * 1.0e6
+    }
+
+    /// Returns the rate of change of the shift, in hertz per second.
+    ///
+    /// This is what a Doppler tracker must follow, and it follows from the
+    /// range acceleration.
+    ///
+    /// Returns `None` when the acceleration is not finite.
+    pub fn shift_rate(self, range_acceleration: f64) -> Option<f64> {
+        if !range_acceleration.is_finite() {
+            return None;
+        }
+        Some(-self.carrier * range_acceleration / SPEED_OF_LIGHT)
+    }
+
+    /// Returns `true` when the shift is a redshift, meaning the terminals
+    /// are receding and the received frequency is lower.
+    pub fn is_receding(self) -> bool {
+        self.range_rate > 0.0
+    }
+}
+
+impl fmt::Display for DopplerShift {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(
+            f,
+            "{} Hz on {} Hz ({:.3} ppm)",
+            self.shift(),
+            self.carrier,
+            self.shift_in_ppm()
+        )
+    }
+}
+
+/// Returns the angle between the line of sight and the velocity, in radians.
+///
+/// This is the complement of the usual elevation angle within the plane
+/// containing both, and it is useful when the range rate is known and the
+/// velocity from the orbit is known: the angle follows from the ratio of the
+/// two.
+///
+/// The result lies in `[0, pi]`.
+///
+/// Returns `None` when the velocity has zero length or either argument is
+/// not finite.
+pub fn angle_from_range_rate(velocity: Vec3, range_rate: f64) -> Option<f64> {
+    if !velocity.is_finite() || !range_rate.is_finite() {
+        return None;
+    }
+    let speed = velocity.norm();
+    if speed == 0.0 {
+        return None;
+    }
+    Some((range_rate / speed).clamp(-1.0, 1.0).acos())
+}
+
+/// Returns the largest possible shift magnitude for a carrier and a speed,
+/// in hertz.
+///
+/// This is the shift when the motion is entirely along the line of sight,
+/// which is the width of the acquisition range a receiver must cover. It is
+/// half the total span between the approach and the recession extremes.
+///
+/// The result is returned as a positive magnitude. Returns `None` when the
+/// carrier or speed is not positive and finite.
+pub fn maximum_shift(carrier_frequency: f64, speed: f64) -> Option<f64> {
+    if !carrier_frequency.is_finite() || carrier_frequency <= 0.0 {
+        return None;
+    }
+    if !speed.is_finite() || speed < 0.0 {
+        return None;
+    }
+    Some(carrier_frequency * speed / SPEED_OF_LIGHT)
+}
+
+/// Returns the total span of frequencies a receiver must accept to cover
+/// every range rate up to the given speed.
+///
+/// Approach and recession give shifts of opposite sign, so the span is twice
+/// the maximum shift.
+///
+/// Returns `None` when the carrier or speed is not usable.
+pub fn acquisition_span(carrier_frequency: f64, speed: f64) -> Option<f64> {
+    maximum_shift(carrier_frequency, speed).map(|shift| 2.0 * shift)
+}

crates/rf/apsis-doppler/tests/shift.rs

@@ -0,0 +1,271 @@
+//! Tests for the Doppler shift.
+
+use apsis_doppler::{
+    DopplerShift, SPEED_OF_LIGHT, acquisition_span, angle_from_range_rate, classical_shift,
+    maximum_shift, range_rate_from_velocity, relativistic_doppler_factor,
+};
+use apsis_math::Vec3;
+use core::f64::consts::PI;
+
+const TOLERANCE: f64 = 1e-9;
+
+#[test]
+fn range_rate_is_the_projection_of_velocity_on_the_line_of_sight() {
+    // Moving away along the line of sight.
+    assert_eq!(
+        range_rate_from_velocity(Vec3::new(1_000.0, 0.0, 0.0), Vec3::new(10.0, 0.0, 0.0)),
+        Some(10.0)
+    );
+    // Moving across the line of sight gives no range rate.
+    assert_eq!(
+        range_rate_from_velocity(Vec3::new(1_000.0, 0.0, 0.0), Vec3::new(0.0, 10.0, 0.0)),
+        Some(0.0)
+    );
+    // Approaching along the line of sight.
+    assert_eq!(
+        range_rate_from_velocity(Vec3::new(1_000.0, 0.0, 0.0), Vec3::new(-10.0, 0.0, 0.0)),
+        Some(-10.0)
+    );
+}
+
+#[test]
+fn range_rate_uses_the_magnitude_of_the_separation() {
+    // The velocity is doubled but the separation is also doubled, so the
+    // projection along the diagonal works out to the same fraction.
+    let rate = range_rate_from_velocity(Vec3::new(3.0, 4.0, 0.0), Vec3::new(6.0, 8.0, 0.0))
+        .expect("non-zero separation");
+    // The unit line of sight is (0.6, 0.8, 0), so the rate is
+    // 0.6 * 6 + 0.8 * 8 = 10.
+    assert!((rate - 10.0).abs() < TOLERANCE);
+}
+
+#[test]
+fn a_zero_separation_has_no_range_rate() {
+    assert_eq!(
+        range_rate_from_velocity(Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0)),
+        None
+    );
+}
+
+#[test]
+fn non_finite_arguments_are_rejected() {
+    assert_eq!(
+        range_rate_from_velocity(Vec3::splat(f64::NAN), Vec3::ZERO),
+        None
+    );
+    assert_eq!(
+        range_rate_from_velocity(Vec3::X, Vec3::splat(f64::INFINITY)),
+        None
+    );
+}
+
+#[test]
+fn the_classical_shift_is_negative_for_a_receding_terminal() {
+    let shift = classical_shift(12.0e9, 1_000.0).expect("valid arguments");
+    assert!(shift < 0.0);
+    assert!((shift + 12.0e9 * 1_000.0 / SPEED_OF_LIGHT).abs() < 1e-3);
+}
+
+#[test]
+fn a_known_range_rate_gives_a_known_shift() {
+    // A range rate of 7 km/s at 12 GHz gives about 280.2 kHz.
+    let shift = classical_shift(12.0e9, 7_000.0).expect("valid");
+    assert!((shift + 280_194.0).abs() < 1.0, "shift {shift}");
+    assert!((shift.abs() / 12.0e9 * 1.0e6 - 23.35).abs() < 0.01);
+}
+
+#[test]
+fn approaching_and_receding_shifts_are_mirror_images() {
+    let receding = classical_shift(1.0e9, 5_000.0).expect("valid");
+    let approaching = classical_shift(1.0e9, -5_000.0).expect("valid");
+    assert!((receding + approaching).abs() < 1e-9);
+}
+
+#[test]
+fn a_zero_range_rate_gives_no_shift() {
+    assert_eq!(classical_shift(12.0e9, 0.0), Some(0.0));
+}
+
+#[test]
+fn invalid_arguments_are_rejected() {
+    assert_eq!(classical_shift(0.0, 1_000.0), None);
+    assert_eq!(classical_shift(-1.0, 1_000.0), None);
+    assert_eq!(classical_shift(12.0e9, f64::NAN), None);
+    assert_eq!(classical_shift(f64::NAN, 1_000.0), None);
+    assert_eq!(classical_shift(12.0e9, f64::INFINITY), None);
+}
+
+#[test]
+fn the_relativistic_factor_is_one_at_rest() {
+    let factor = relativistic_doppler_factor(0.0).expect("valid");
+    assert!((factor - 1.0).abs() < TOLERANCE);
+}
+
+#[test]
+fn the_relativistic_factor_opposes_the_classical_shift() {
+    // For a receding terminal the factor is below one, which reinforces the
+    // classical redshift.
+    let receding = relativistic_doppler_factor(7_000.0).expect("valid");
+    assert!(receding < 1.0);
+    // For an approaching terminal it is above one.
+    let approaching = relativistic_doppler_factor(-7_000.0).expect("valid");
+    assert!(approaching > 1.0);
+    // The two are reciprocals, because swapping the roles of the terminals
+    // inverts the ratio.
+    assert!((receding * approaching - 1.0).abs() < 1e-15);
+}
+
+#[test]
+fn the_relativistic_factor_contains_the_classical_shift() {
+    // The factor's expansion begins with 1 - beta, which is the classical
+    // result, so the two treatments must agree to first order. They differ
+    // only by terms of order beta squared.
+    let beta = 7_000.0 / SPEED_OF_LIGHT;
+    let factor = relativistic_doppler_factor(7_000.0).expect("valid");
+    let difference = (factor - 1.0 + beta).abs();
+    // The leading difference is beta^2 / 2, about 2.7e-10.
+    assert!(difference < beta * beta, "difference {difference}");
+    assert!((difference - beta * beta / 2.0).abs() < 1e-12);
+}
+
+#[test]
+fn the_two_treatments_agree_to_second_order() {
+    // The classical shift and the relativistic one differ by a fraction of
+    // the carrier that is of order beta squared.
+    let carrier = 12.0e9;
+    let shift = DopplerShift::from_range_rate(carrier, 7_000.0).expect("valid");
+    let classical = shift.shift();
+    let relativistic = shift.shift_including_relativity().expect("valid");
+    let relative_difference = (relativistic - classical).abs() / carrier;
+    let beta = 7_000.0 / SPEED_OF_LIGHT;
+    assert!(
+        relative_difference < beta * beta,
+        "difference {relative_difference}"
+    );
+    // The absolute difference is the second order term, about 3.3 Hz at
+    // 12 GHz: tiny compared with the 280 kHz classical shift, but real.
+    assert!(
+        (relativistic - classical).abs() > 1.0,
+        "difference {}",
+        (relativistic - classical).abs()
+    );
+    assert!((relativistic - classical).abs() < 10.0);
+}
+
+#[test]
+fn a_speed_of_light_range_rate_is_rejected() {
+    assert_eq!(relativistic_doppler_factor(SPEED_OF_LIGHT), None);
+    assert_eq!(relativistic_doppler_factor(-SPEED_OF_LIGHT), None);
+    assert_eq!(relativistic_doppler_factor(SPEED_OF_LIGHT * 2.0), None);
+    assert_eq!(relativistic_doppler_factor(f64::NAN), None);
+}
+
+#[test]
+fn the_shift_type_reports_the_carrier_and_rate() {
+    let shift = DopplerShift::from_range_rate(12.0e9, 1_000.0).expect("valid");
+    assert_eq!(shift.carrier_frequency(), 12.0e9);
+    assert_eq!(shift.range_rate(), 1_000.0);
+    assert!(shift.is_receding());
+}
+
+#[test]
+fn the_shift_type_converts_between_frequencies() {
+    let shift = DopplerShift::from_range_rate(12.0e9, 1_000.0).expect("valid");
+    let received = shift.shifted_frequency();
+    assert!(received < 12.0e9);
+    assert!((received - (12.0e9 + shift.shift())).abs() < 1e-6);
+}
+
+#[test]
+fn an_approaching_terminal_has_a_higher_received_frequency() {
+    let shift = DopplerShift::from_range_rate(12.0e9, -1_000.0).expect("valid");
+    assert!(!shift.is_receding());
+    assert!(shift.shifted_frequency() > 12.0e9);
+}
+
+#[test]
+fn the_shift_in_parts_per_million_is_consistent() {
+    let shift = DopplerShift::from_range_rate(12.0e9, 7_000.0).expect("valid");
+    let ppm = shift.shift_in_ppm();
+    assert!((ppm + 7_000.0 / SPEED_OF_LIGHT * 1.0e6).abs() < 1e-9);
+    assert!((ppm + 23.349).abs() < 0.01, "ppm {ppm}");
+}
+
+#[test]
+fn the_relativistic_shift_is_close_to_the_classical_one() {
+    let shift = DopplerShift::from_range_rate(12.0e9, 7_000.0).expect("valid");
+    let classical = shift.shift();
+    let total = shift.shift_including_relativity().expect("valid");
+    // Both are negative for a receding terminal and agree in sign.
+    assert!(classical < 0.0);
+    assert!(total < 0.0);
+    // The second order difference is a few hertz on a shift of a few
+    // hundred kilohertz.
+    assert!((total - classical).abs() > 1.0);
+    assert!((total - classical).abs() < 10.0);
+}
+
+#[test]
+fn invalid_construction_is_rejected() {
+    assert_eq!(DopplerShift::from_range_rate(0.0, 1_000.0), None);
+    assert_eq!(DopplerShift::from_range_rate(12.0e9, f64::NAN), None);
+    assert_eq!(DopplerShift::from_range_rate(-12.0e9, 1_000.0), None);
+}
+
+#[test]
+fn the_shift_rate_follows_the_range_acceleration() {
+    let shift = DopplerShift::from_range_rate(12.0e9, 1_000.0).expect("valid");
+    // A range acceleration of 10 m/s^2 changes the shift by
+    // 12 GHz * 10 / c, about 0.4 Hz per second.
+    let rate = shift.shift_rate(10.0).expect("valid");
+    assert!((rate + 400.3).abs() < 0.1, "rate {rate}");
+    assert_eq!(shift.shift_rate(f64::NAN), None);
+}
+
+#[test]
+fn the_angle_from_a_range_rate_is_the_angle_to_the_velocity() {
+    let velocity = Vec3::new(0.0, 1_000.0, 0.0);
+    // The full speed along the line of sight is zero degrees.
+    let along = angle_from_range_rate(velocity, 1_000.0).expect("valid");
+    assert!(along.abs() < 1e-9);
+    // The full speed across the line of sight is a right angle.
+    let across = angle_from_range_rate(velocity, 0.0).expect("valid");
+    assert!((across - PI / 2.0).abs() < TOLERANCE);
+    // The opposite direction is a straight angle.
+    let opposite = angle_from_range_rate(velocity, -1_000.0).expect("valid");
+    assert!((opposite - PI).abs() < 1e-9);
+}
+
+#[test]
+fn the_angle_from_a_range_rate_rejects_invalid_input() {
+    assert_eq!(angle_from_range_rate(Vec3::ZERO, 1.0), None);
+    assert_eq!(angle_from_range_rate(Vec3::X, f64::NAN), None);
+}
+
+#[test]
+fn the_maximum_shift_covers_the_line_of_sight_case() {
+    let shift = maximum_shift(12.0e9, 7_000.0).expect("valid");
+    assert!((shift - 12.0e9 * 7_000.0 / SPEED_OF_LIGHT).abs() < 1e-3);
+    // A range rate beyond the speed is not physical, so the speed must be
+    // non-negative and finite.
+    assert_eq!(maximum_shift(12.0e9, -1.0), None);
+    assert_eq!(maximum_shift(-1.0, 7_000.0), None);
+    assert_eq!(maximum_shift(12.0e9, f64::INFINITY), None);
+}
+
+#[test]
+fn the_acquisition_span_is_twice_the_maximum_shift() {
+    let span = acquisition_span(12.0e9, 7_000.0).expect("valid");
+    let maximum = maximum_shift(12.0e9, 7_000.0).expect("valid");
+    assert!((span - 2.0 * maximum).abs() < 1e-3);
+    // For a low Earth orbit at 12 GHz the span is about 560 kHz.
+    assert!((span - 560_388.0).abs() < 2.0, "span {span}");
+}
+
+#[test]
+fn display_reports_the_shift_and_the_carrier() {
+    let shift = DopplerShift::from_range_rate(12.0e9, 1_000.0).expect("valid");
+    let text = shift.to_string();
+    assert!(text.contains("Hz"));
+    assert!(text.contains("ppm"));
+}