siderust 0.9.0

High-precision astronomy and satellite mechanics in Rust.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

// SPDX-License-Identifier: AGPL-3.0-or-later
// Copyright (C) 2026 Vallés Puig, Ramon

//! Built-in Earth gravity-field providers.

use principia::{GravityFieldProvider, PrincipiaError};

use crate::astro::dynamics::units::GravitationalParameter;
use crate::qtty::Kilometers;

use crate::archive::gravity::tables::EGM2008_COEFFS;

/// Earth GM (WGS-84), km³/s².
const GM_EARTH_KM3_S2: f64 = 398_600.441_8;
/// Earth equatorial radius (WGS-84 / EGM2008), km.
const R_EARTH_KM: f64 = 6_378.137;

/// Trivial degree-0 Earth gravity field for two-body propagation.
#[derive(Debug, Clone, Copy)]
pub struct TwoBodyEarth;

impl GravityFieldProvider for TwoBodyEarth {
    fn mu(&self) -> GravitationalParameter {
        GravitationalParameter::new(GM_EARTH_KM3_S2)
    }

    fn reference_radius(&self) -> Kilometers {
        Kilometers::new(R_EARTH_KM)
    }

    fn max_degree(&self) -> usize {
        0
    }

    fn max_order(&self) -> usize {
        0
    }

    fn c_normalized(&self, n: usize, _m: usize) -> Result<f64, PrincipiaError> {
        Ok(if n == 0 { 1.0 } else { 0.0 })
    }

    fn s_normalized(&self, _n: usize, _m: usize) -> Result<f64, PrincipiaError> {
        Ok(0.0)
    }
}

/// Earth geopotential through degree/order 4 using EGM2008 normalized coefficients.
#[derive(Debug, Clone, Copy)]
pub struct LowDegreeEarth;

impl LowDegreeEarth {
    fn lookup(n: usize, m: usize) -> (f64, f64) {
        if n == 0 && m == 0 {
            return (1.0, 0.0);
        }
        for &(tn, tm, c, s) in EGM2008_COEFFS {
            if tn == n && tm == m {
                return (c, s);
            }
        }
        (0.0, 0.0)
    }
}

impl GravityFieldProvider for LowDegreeEarth {
    fn mu(&self) -> GravitationalParameter {
        GravitationalParameter::new(GM_EARTH_KM3_S2)
    }

    fn reference_radius(&self) -> Kilometers {
        Kilometers::new(R_EARTH_KM)
    }

    fn max_degree(&self) -> usize {
        4
    }

    fn max_order(&self) -> usize {
        4
    }

    fn c_normalized(&self, n: usize, m: usize) -> Result<f64, PrincipiaError> {
        Ok(if n > self.max_degree() || m > n {
            0.0
        } else {
            Self::lookup(n, m).0
        })
    }

    fn s_normalized(&self, n: usize, m: usize) -> Result<f64, PrincipiaError> {
        Ok(if n > self.max_degree() || m > n {
            0.0
        } else {
            Self::lookup(n, m).1
        })
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn two_body_earth_c00() {
        assert_eq!(TwoBodyEarth.c_normalized(0, 0).unwrap(), 1.0);
        assert_eq!(TwoBodyEarth.c_normalized(2, 0).unwrap(), 0.0);
    }

    #[test]
    fn low_degree_earth_c20_consistent_with_j2() {
        let c20 = LowDegreeEarth.c_normalized(2, 0).unwrap();
        let j2_recovered = -c20 * 5.0_f64.sqrt();
        let j2_expected = 1.082_626_68e-3;
        assert!((j2_recovered - j2_expected).abs() / j2_expected < 1e-7);
    }
}