siderust 0.7.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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255

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

//! # Angular Separation Audit Tests
//!
//! Verifies the `angular_separation` function across all frame variants:
//!
//! - **Formula**: all frame variants share the same underlying
//!   `affn::spherical::angular_separation_impl` which uses the Vincenty
//!   formula `atan2(|n×m|, n·m)` — NOT naive `acos(n·m)`.
//! - **Return type**: typed `Degrees` (`qtty::angular::Degrees`), not raw `f64`.
//! - **Identity**: `a.angular_separation(&a) == Degrees::new(0.0)` exactly.
//! - **Antipode**: `a.angular_separation(&antipode(a)) == 180°`.
//! - **Symmetry**: `sep(a,b) ≈ sep(b,a)` within a small tolerance
//!   (the computation paths differ in FP order when swapping arguments).
//! - **Cross-frame**: angular separation is a geometric invariant; equatorial
//!   and ecliptic frames must agree to 1 × 10⁻¹² rad.
//! - **Precision near 0**: two points 10⁻⁹ rad apart must yield a result > 0
//!   and within 1 % of the expected value.  Naive `acos` would lose all
//!   significant digits here.
//! - **Precision near π**: two near-antipodal points must yield a result
//!   within 10⁻¹² rad of the expected value.

use std::f64::consts::PI;

use siderust::coordinates::spherical::direction;
use siderust::coordinates::transform::TransformFrame;
use siderust::qtty::*;

// ---------------------------------------------------------------------------
// Internal helper: Lehmer LCG for deterministic pseudo-random test angles.
// ---------------------------------------------------------------------------
fn lcg_next(state: &mut u64) -> f64 {
    *state = state
        .wrapping_mul(6_364_136_223_846_793_005)
        .wrapping_add(1_442_695_040_888_963_407);
    // Upper 31 bits → [0, 1)
    ((*state >> 33) as f64) / (0x8000_0000_u64 as f64)
}

// ---------------------------------------------------------------------------
// 1. Identity – identical directions must yield exactly 0°
// ---------------------------------------------------------------------------

#[test]
fn angular_separation_identity_is_exact_zero() {
    let cases: &[(f64, f64)] = &[
        (0.0, 0.0),
        (37.9546, 89.2641),  // near Polaris
        (279.2347, 38.7836), // near Vega
        (180.0, -45.0),
        (359.9999, -89.9999),
    ];

    for &(ra, dec) in cases {
        let d = direction::EquatorialMeanJ2000::new(Degrees::new(ra), Degrees::new(dec));
        let sep = d.angular_separation(&d).value();
        assert_eq!(
            sep, 0.0,
            "identity failed for (RA={ra}, Dec={dec}): got {sep}"
        );
    }
}

// ---------------------------------------------------------------------------
// 2. Antipode – must return exactly 180° (π rad) within floating-point limits
// ---------------------------------------------------------------------------

#[test]
fn angular_separation_poles_are_antipodal() {
    let north = direction::EquatorialMeanJ2000::new(0.0 * DEG, 90.0 * DEG);
    let south = direction::EquatorialMeanJ2000::new(0.0 * DEG, -90.0 * DEG);
    let sep = north.angular_separation(&south).value();
    assert!(
        (sep - 180.0).abs() < 1e-12,
        "pole antipode: expected 180°, got {sep}"
    );
}

#[test]
fn angular_separation_general_antipode() {
    // For a direction at (ra, dec), its antipode is (ra + 180°, -dec).
    let cases: &[(f64, f64)] = &[
        (45.0, 30.0),
        (270.0, -60.0),
        (120.0, 0.0),
        (0.0, 45.0),
        (10.0, 89.0), // near pole
    ];

    for &(ra, dec) in cases {
        let a = direction::EquatorialMeanJ2000::new(Degrees::new(ra), Degrees::new(dec));
        let antipode = direction::EquatorialMeanJ2000::new(
            Degrees::new((ra + 180.0) % 360.0),
            Degrees::new(-dec),
        );
        let sep = a.angular_separation(&antipode).value();
        assert!(
            (sep - 180.0).abs() < 1e-11,
            "antipode ({ra},{dec}): expected 180°, got {sep}"
        );
    }
}

// ---------------------------------------------------------------------------
// 3. Symmetry – sep(a,b) ≈ sep(b,a) over 100 seeded random pairs.
//    The Vincenty formula is mathematically symmetric; in floating-point the
//    two evaluation paths differ slightly.  We accept up to 4 ulps.
// ---------------------------------------------------------------------------

#[test]
fn angular_separation_symmetry_100_random_pairs() {
    let mut state: u64 = 0xdead_beef_cafe_f00d;

    for i in 0u64..100 {
        let ra1 = Degrees::new(lcg_next(&mut state) * 360.0);
        let dec1 = Degrees::new(lcg_next(&mut state) * 180.0 - 90.0);
        let ra2 = Degrees::new(lcg_next(&mut state) * 360.0);
        let dec2 = Degrees::new(lcg_next(&mut state) * 180.0 - 90.0);

        let a = direction::EquatorialMeanJ2000::new(ra1, dec1);
        let b = direction::EquatorialMeanJ2000::new(ra2, dec2);

        let sep_ab = a.angular_separation(&b).value();
        let sep_ba = b.angular_separation(&a).value();

        // Angular separation is always non-negative; compare raw bit patterns.
        let ba = sep_ab.to_bits();
        let bb = sep_ba.to_bits();
        let ulp_diff = ba.abs_diff(bb);

        assert!(
            ulp_diff <= 4,
            "pair {i}: symmetry violated: sep(a,b)={sep_ab:.15e} sep(b,a)={sep_ba:.15e} ulps={ulp_diff}"
        );
    }
}

// ---------------------------------------------------------------------------
// 4. Cross-frame agreement – Vega vs Polaris: equatorial == ecliptic
//    Angular separation is a geometric invariant under any rigid rotation.
// ---------------------------------------------------------------------------

#[test]
fn angular_separation_cross_frame_equatorial_vs_ecliptic() {
    // Vega  J2000: RA=279.2347°, Dec=+38.7836°
    // Polaris J2000: RA=37.9546°, Dec=+89.2641°
    let vega_eq =
        direction::EquatorialMeanJ2000::new(Degrees::new(279.2347), Degrees::new(38.7836));
    let polaris_eq =
        direction::EquatorialMeanJ2000::new(Degrees::new(37.9546), Degrees::new(89.2641));

    let sep_eq_deg = vega_eq.angular_separation(&polaris_eq).value();

    // Convert both directions to ecliptic mean J2000.
    let vega_ecl: direction::EclipticMeanJ2000 = TransformFrame::to_frame(&vega_eq);
    let polaris_ecl: direction::EclipticMeanJ2000 = TransformFrame::to_frame(&polaris_eq);

    let sep_ecl_deg = vega_ecl.angular_separation(&polaris_ecl).value();

    let diff_rad = (sep_eq_deg - sep_ecl_deg).abs().to_radians();
    assert!(
        diff_rad < 1e-12,
        "cross-frame disagreement: equatorial={sep_eq_deg:.15}° ecliptic={sep_ecl_deg:.15}° diff={diff_rad:.2e} rad"
    );
}

// ---------------------------------------------------------------------------
// 5. Precision near 0 – two points 1×10⁻⁹ rad apart.
//    Naive acos(dot) loses all significant figures here; Vincenty does not.
// ---------------------------------------------------------------------------

#[test]
fn angular_separation_precision_near_zero() {
    let epsilon_rad = 1e-9_f64;
    let epsilon_deg = epsilon_rad.to_degrees(); // ≈ 5.73×10⁻⁸ °

    // Separation entirely in RA along the equator.
    let a = direction::EquatorialMeanJ2000::new(0.0 * DEG, 0.0 * DEG);
    let b = direction::EquatorialMeanJ2000::new(Degrees::new(epsilon_deg), 0.0 * DEG);

    let sep_deg = a.angular_separation(&b).value();

    assert!(
        sep_deg > 0.0,
        "precision near 0: Vincenty returned exactly 0 (expected ~{epsilon_deg:.2e}°)"
    );
    let rel_err = (sep_deg - epsilon_deg).abs() / epsilon_deg;
    assert!(
        rel_err < 0.01,
        "precision near 0: relative error {rel_err:.2e} (expected < 1 %)"
    );
}

// ---------------------------------------------------------------------------
// 6. Precision near π – two near-antipodal points separated by (π − 10⁻⁹) rad.
//    Vincenty handles this without precision loss.
// ---------------------------------------------------------------------------

#[test]
fn angular_separation_precision_near_pi() {
    let delta_rad = 1e-9_f64; // deviation from exact antipode
    let delta_deg = delta_rad.to_degrees();

    // North pole, and a point delta_rad above the south pole (dec = -(90 - δ)°).
    let north = direction::EquatorialMeanJ2000::new(0.0 * DEG, 90.0 * DEG);
    let near_south =
        direction::EquatorialMeanJ2000::new(0.0 * DEG, Degrees::new(-(90.0 - delta_deg)));

    let sep_rad = north.angular_separation(&near_south).value().to_radians();
    let expected_rad = PI - delta_rad;
    let error_rad = (sep_rad - expected_rad).abs();

    assert!(
        error_rad < 1e-12,
        "precision near pi: error={error_rad:.2e} rad (expected {expected_rad:.15}, got {sep_rad:.15})"
    );
}

// ---------------------------------------------------------------------------
// 7. Galactic and Horizontal frames share the same implementation.
//    Quick sanity checks to confirm the formula applies there too.
// ---------------------------------------------------------------------------

#[test]
fn angular_separation_galactic_identity_and_quarter_turn() {
    // Identity
    let gc = direction::Galactic::new(120.0 * DEG, 45.0 * DEG);
    assert_eq!(gc.angular_separation(&gc).value(), 0.0, "galactic identity");

    // Two galactic directions 90° apart in longitude at b=0
    let a = direction::Galactic::new(0.0 * DEG, 0.0 * DEG);
    let b = direction::Galactic::new(90.0 * DEG, 0.0 * DEG);
    let sep = a.angular_separation(&b).value();
    assert!(
        (sep - 90.0).abs() < 1e-12,
        "galactic 90° separation: expected 90°, got {sep}"
    );
}

#[test]
fn angular_separation_horizontal_identity_and_quarter_turn() {
    // Identity
    let h = direction::Horizontal::new(180.0 * DEG, 30.0 * DEG);
    assert_eq!(h.angular_separation(&h).value(), 0.0, "horizontal identity");

    // Zenith (alt=90) and horizon (alt=0) are 90° apart
    let zenith = direction::Horizontal::new(0.0 * DEG, 90.0 * DEG);
    let horizon = direction::Horizontal::new(0.0 * DEG, 0.0 * DEG);
    let sep = zenith.angular_separation(&horizon).value();
    assert!(
        (sep - 90.0).abs() < 1e-12,
        "horizontal 90° separation: expected 90°, got {sep}"
    );
}