From 85f3d8198c2914628f6431bb584d5cb2d5920ff1 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Sat, 26 Sep 2026 14:40:12 +0800 Subject: [PATCH] rf: add Doppler shift and range rate models --- Cargo.lock | 4 + crates/rf/apsis-doppler/Cargo.toml | 2 + crates/rf/apsis-doppler/src/lib.rs | 45 ++++ crates/rf/apsis-doppler/src/shift.rs | 251 +++++++++++++++++++++++ crates/rf/apsis-doppler/tests/shift.rs | 271 +++++++++++++++++++++++++ 5 files changed, 573 insertions(+) create mode 100644 crates/rf/apsis-doppler/src/shift.rs create mode 100644 crates/rf/apsis-doppler/tests/shift.rs diff --git a/Cargo.lock b/Cargo.lock index 916569c..48100c8 100644 --- a/Cargo.lock +++ b/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" diff --git a/crates/rf/apsis-doppler/Cargo.toml b/crates/rf/apsis-doppler/Cargo.toml index a8e0f3f..cb9ad9f 100644 --- a/crates/rf/apsis-doppler/Cargo.toml +++ b/crates/rf/apsis-doppler/Cargo.toml @@ -11,3 +11,5 @@ publish = false workspace = true [dependencies] +apsis-core.workspace = true +apsis-math.workspace = true diff --git a/crates/rf/apsis-doppler/src/lib.rs b/crates/rf/apsis-doppler/src/lib.rs index 87d95e8..0f6eb63 100644 --- a/crates/rf/apsis-doppler/src/lib.rs +++ b/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; diff --git a/crates/rf/apsis-doppler/src/shift.rs b/crates/rf/apsis-doppler/src/shift.rs new file mode 100644 index 0000000..0511635 --- /dev/null +++ b/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 { + 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 { + 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 { + 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 { + 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 { + 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 { + 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 { + 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 { + 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 { + maximum_shift(carrier_frequency, speed).map(|shift| 2.0 * shift) +} diff --git a/crates/rf/apsis-doppler/tests/shift.rs b/crates/rf/apsis-doppler/tests/shift.rs new file mode 100644 index 0000000..97fc399 --- /dev/null +++ b/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")); +}