astrodynamics 0.13.0

Numerical astrodynamics engine for orbit propagation, force models, and flight-dynamics primitives
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

//! Golden validation of [`sun_moon_ecef`] against a Skyfield + JPL DE440 ITRS
//! reference (see `tests/fixtures/bodies/gen_sun_moon_golden.py`).
//!
//! `sun_moon_ecef` is the RTKLIB-style low-precision analytic series rotated
//! GCRS -> ITRS by the crate's IAU 2000A/2006 transform, so the bar is the
//! model's own accuracy (sub-degree direction), not bit-exactness. The value of
//! the golden is catching frame / units / series-coefficient regressions: a
//! wrong ECI-frame interpretation or a flipped/duplicated series term shows up
//! as degree-to-huge ITRS deviations that this test gates out.

use std::path::PathBuf;

use serde_json::Value;

use super::sun_moon_ecef;
use crate::time::scales::TimeScales;

fn fixture_path() -> PathBuf {
    PathBuf::from(env!("CARGO_MANIFEST_DIR"))
        .join("tests/fixtures/bodies/sun_moon_skyfield_golden.json")
        .canonicalize()
        .expect("locate sun_moon_skyfield_golden.json")
}

fn vec3(v: &Value) -> [f64; 3] {
    let a = v.as_array().expect("vec3 array");
    [
        a[0].as_f64().unwrap(),
        a[1].as_f64().unwrap(),
        a[2].as_f64().unwrap(),
    ]
}

fn norm(v: &[f64; 3]) -> f64 {
    (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt()
}

fn angle_deg(a: &[f64; 3], b: &[f64; 3]) -> f64 {
    let dot = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
    let c = (dot / (norm(a) * norm(b))).clamp(-1.0, 1.0);
    c.acos().to_degrees()
}

#[test]
fn sun_moon_ecef_matches_skyfield_de440_golden() {
    // Tolerances justified by the analytic model accuracy (sub-degree direction,
    // sub-percent distance) plus margin; tight enough that a frame/units/series
    // regression (degrees+) fails, loose enough that the faithful model passes.
    const SUN_DIR_TOL_DEG: f64 = 0.10;
    const SUN_DIST_TOL: f64 = 0.003; // 0.3%
    const MOON_DIR_TOL_DEG: f64 = 0.60;
    const MOON_DIST_TOL: f64 = 0.012; // 1.2%

    let raw = std::fs::read_to_string(fixture_path()).expect("read golden");
    let doc: Value = serde_json::from_str(&raw).expect("parse golden");
    let cases = doc["cases"].as_array().expect("cases");

    let (mut max_sun_dir, mut max_sun_dist) = (0.0_f64, 0.0_f64);
    let (mut max_moon_dir, mut max_moon_dist) = (0.0_f64, 0.0_f64);
    let mut failures: Vec<String> = Vec::new();

    for case in cases {
        let utc = &case["utc"];
        let ts = TimeScales::from_utc(
            utc["year"].as_i64().unwrap() as i32,
            utc["month"].as_i64().unwrap() as i32,
            utc["day"].as_i64().unwrap() as i32,
            utc["hour"].as_i64().unwrap() as i32,
            utc["minute"].as_i64().unwrap() as i32,
            utc["second"].as_f64().unwrap(),
        );
        let sm = sun_moon_ecef(&ts);

        let sun_ref = vec3(&case["sun_itrs_m"]);
        let moon_ref = vec3(&case["moon_itrs_m"]);

        let sun_dir = angle_deg(&sm.sun, &sun_ref);
        let moon_dir = angle_deg(&sm.moon, &moon_ref);
        let sun_dist = (norm(&sm.sun) - norm(&sun_ref)).abs() / norm(&sun_ref);
        let moon_dist = (norm(&sm.moon) - norm(&moon_ref)).abs() / norm(&moon_ref);

        max_sun_dir = max_sun_dir.max(sun_dir);
        max_moon_dir = max_moon_dir.max(moon_dir);
        max_sun_dist = max_sun_dist.max(sun_dist);
        max_moon_dist = max_moon_dist.max(moon_dist);

        if sun_dir > SUN_DIR_TOL_DEG
            || sun_dist > SUN_DIST_TOL
            || moon_dir > MOON_DIR_TOL_DEG
            || moon_dist > MOON_DIST_TOL
        {
            failures.push(format!(
                "{}: sun_dir={sun_dir:.4}deg sun_dist={sun_dist:.4e} moon_dir={moon_dir:.4}deg moon_dist={moon_dist:.4e}",
                case["utc"]
            ));
        }
    }

    println!(
        "sun_moon golden max dev: sun_dir={max_sun_dir:.4}deg sun_dist={max_sun_dist:.4e} moon_dir={max_moon_dir:.4}deg moon_dist={max_moon_dist:.4e}"
    );

    assert!(
        failures.is_empty(),
        "sun_moon_ecef golden out of tolerance:\n{}",
        failures.join("\n")
    );
}