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

rf: add Doppler shift and range rate models

This commit is contained in:
ApsisRS contributors 2026-09-26 14:40:12 +08:00
parent 7d08f900a9
commit b3d3071778
5 changed files with 573 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Cargo.lock generated
View File

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

2
crates/rf/apsis-doppler/Cargo.toml
View File

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

45
crates/rf/apsis-doppler/src/lib.rs
View File

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

251
crates/rf/apsis-doppler/src/shift.rs Normal file
View File

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

271
crates/rf/apsis-doppler/tests/shift.rs Normal file
View File

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