satkit 0.16.2

Satellite Toolkit
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

use crate::mathtypes::*;

// Equation 3.37 in Montenbruck & Gill
pub fn point_gravity(
    r: &Vector3, // object
    s: &Vector3, // distant attractor
    mu: f64,
) -> Vector3 {
    let sr = s - r;
    let srnorm2 = sr.norm_squared();
    let srnorm = srnorm2.sqrt();
    let snorm2 = s.norm_squared();
    let snorm = snorm2.sqrt();
    mu * (sr / (srnorm * srnorm2) - s / (snorm * snorm2))
}

// Return tuple with point gravity force and
// point gravity partial (da/dr)
// Equation 3.37 in Montenbruck & Gill for point gravity
// Equation 7.75 in Montenbruck & Gill for partials

pub fn point_gravity_and_partials(
    r: &Vector3, // object
    s: &Vector3, // distant attractor
    mu: f64,
) -> (Vector3, Matrix3) {
    let rs = r - s;
    let rsnorm2 = rs.norm_squared();
    let rsnorm = rsnorm2.sqrt();
    let snorm2 = s.norm_squared();
    let snorm = snorm2.sqrt();
    let rsnorm3 = rsnorm2 * rsnorm;
    (
        -mu * (rs / rsnorm3 + s / (snorm * snorm2)),
        -mu * (Matrix3::eye() / rsnorm3
            - 3.0 * rs * rs.transpose() / (rsnorm2 * rsnorm3)),
    )
}

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

    #[test]
    fn test_point_gravity_known() {
        // Moon at ~384,400 km from Earth center, satellite at GEO (~42,164 km)
        let mu_moon = 4.9048695e12; // m³/s²
        let s_moon = numeris::vector![384_400.0e3, 0.0, 0.0]; // Moon position
        let r_sat = numeris::vector![42_164.0e3, 0.0, 0.0]; // GEO satellite

        let accel = point_gravity(&r_sat, &s_moon, mu_moon);
        // Lunar perturbation at GEO should be ~1e-6 m/s² order of magnitude
        let mag = accel.norm();
        assert!(
            mag > 1.0e-7 && mag < 1.0e-4,
            "Lunar perturbation at GEO = {:.3e}, expected ~1e-6",
            mag
        );
    }

    #[test]
    fn test_point_gravity_partials() {
        let mu = 4.9048695e12;
        let s = numeris::vector![384_400.0e3, 50_000.0e3, 20_000.0e3];
        let r = numeris::vector![42_164.0e3, 1000.0e3, 500.0e3];
        let dr = numeris::vector![10.0, 20.0, -15.0];

        let (accel0, partials) = point_gravity_and_partials(&r, &s, mu);
        let accel1 = point_gravity(&(r + dr), &s, mu);

        // accel(r+dr) ≈ accel(r) + (da/dr)*dr
        let predicted = accel0 + partials * dr;
        let err = (accel1 - predicted).norm() / accel1.norm();
        assert!(
            err < 1.0e-6,
            "Partial derivative error = {:.3e}, expected < 1e-6",
            err
        );
    }
}