astrodynamics 0.11.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
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

//! Analytic low-precision Sun and Moon positions in ECI and ECEF.
//!
//! The Earth-centered-inertial series below are the standard low-precision
//! analytic solar/lunar formulae (Montenbruck & Gill, "Satellite Orbits",
//! sec. 3.3.2 / 3.3.3) as carried by common open GNSS code. They take the time
//! argument `t` in Julian centuries of Terrestrial Time since the J2000.0 epoch
//! (2000-01-01T12:00:00 TT). The five fundamental (Delaunay) arguments
//! `{l, l', F, D, Omega}` are evaluated by [`ast_args`] from the IAU 1980
//! nutation series coefficients.
//!
//! The series are referred to the mean equator and equinox of date (their mean
//! longitude and obliquity are of-date), so the resulting vectors are rotated
//! into the Earth-fixed (ITRS) frame with
//! [`crate::frames::transforms::mean_of_date_to_itrs_matrix`] (IAU 2000A
//! nutation + GAST). Precession is already implicit in the of-date series and is
//! deliberately NOT re-applied; the full GCRS->ITRS transform would double-count
//! it (about 0.37 deg of Sun-direction error at epoch 2026). This mirrors the
//! GMST/GAST + nutation rotation the analytic references consume the series with.

use crate::frames::transforms::{mat3_vec3_mul, mean_of_date_to_itrs_matrix};
use crate::time::scales::TimeScales;

/// 1 astronomical unit in metres (DE405/IAU value used by the analytic series).
const AU_M: f64 = 149_597_870_691.0;
/// WGS84 Earth semi-major axis (m); the lunar distance series is scaled by it.
const RE_WGS84_M: f64 = 6_378_137.0;
/// Degrees -> radians.
const D2R: f64 = std::f64::consts::PI / 180.0;
/// Arcseconds -> radians.
const AS2R: f64 = D2R / 3600.0;
/// Julian centuries -> J2000 reference Julian Date (TT).
const JD_J2000: f64 = 2_451_545.0;
const DAYS_PER_CENTURY: f64 = 36525.0;

/// Sun and Moon positions (metres) in a single frame.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SunMoon {
    /// Sun position vector (m).
    pub sun: [f64; 3],
    /// Moon position vector (m).
    pub moon: [f64; 3],
}

/// Fundamental (Delaunay) arguments `f = {l, l', F, D, Omega}` in radians for
/// the given time `t` (Julian centuries TT since J2000). IAU 1980 series.
fn ast_args(t: f64) -> [f64; 5] {
    // Coefficients for IAU 1980 nutation: l, l', F, D, OMG.
    const FC: [[f64; 5]; 5] = [
        [
            134.96340251,
            1717915923.2178,
            31.8792,
            0.051635,
            -0.00024470,
        ],
        [357.52910918, 129596581.0481, -0.5532, 0.000136, -0.00001149],
        [
            93.27209062,
            1739527262.8478,
            -12.7512,
            -0.001037,
            0.00000417,
        ],
        [
            297.85019547,
            1602961601.2090,
            -6.3706,
            0.006593,
            -0.00003169,
        ],
        [125.04455501, -6962890.2665, 7.4722, 0.007702, -0.00005939],
    ];

    let mut tt = [0.0_f64; 4];
    tt[0] = t;
    for i in 1..4 {
        tt[i] = tt[i - 1] * t;
    }

    let mut f = [0.0_f64; 5];
    for i in 0..5 {
        let mut v = FC[i][0] * 3600.0;
        for j in 0..4 {
            v += FC[i][j + 1] * tt[j];
        }
        f[i] = (v * AS2R).rem_euclid(2.0 * std::f64::consts::PI);
    }
    f
}

/// Analytic Sun and Moon position in ECI (metres) for time `t` in Julian
/// centuries of TT since J2000.0.
pub fn sun_moon_eci(t: f64) -> SunMoon {
    let f = ast_args(t);

    // Obliquity of the ecliptic (deg -> rad).
    let eps = 23.439291 - 0.0130042 * t;
    let sine = (eps * D2R).sin();
    let cose = (eps * D2R).cos();

    // Sun position in ECI.
    let ms = 357.5277233 + 35999.05034 * t;
    let ls = 280.460
        + 36000.770 * t
        + 1.914666471 * (ms * D2R).sin()
        + 0.019994643 * (2.0 * ms * D2R).sin();
    let rs = AU_M
        * (1.000140612 - 0.016708617 * (ms * D2R).cos() - 0.000139589 * (2.0 * ms * D2R).cos());
    let sinl = (ls * D2R).sin();
    let cosl = (ls * D2R).cos();
    let sun = [rs * cosl, rs * cose * sinl, rs * sine * sinl];

    // Moon position in ECI.
    let lm = 218.32 + 481267.883 * t + 6.29 * f[0].sin() - 1.27 * (f[0] - 2.0 * f[3]).sin()
        + 0.66 * (2.0 * f[3]).sin()
        + 0.21 * (2.0 * f[0]).sin()
        - 0.19 * f[1].sin()
        - 0.11 * (2.0 * f[2]).sin();
    let pm = 5.13 * f[2].sin() + 0.28 * (f[0] + f[2]).sin()
        - 0.28 * (f[2] - f[0]).sin()
        - 0.17 * (f[2] - 2.0 * f[3]).sin();
    let rm = RE_WGS84_M
        / ((0.9508
            + 0.0518 * f[0].cos()
            + 0.0095 * (f[0] - 2.0 * f[3]).cos()
            + 0.0078 * (2.0 * f[3]).cos()
            + 0.0028 * (2.0 * f[0]).cos())
            * D2R)
            .sin();
    let sinlm = (lm * D2R).sin();
    let coslm = (lm * D2R).cos();
    let sinp = (pm * D2R).sin();
    let cosp = (pm * D2R).cos();
    let moon = [
        rm * cosp * coslm,
        rm * (cose * cosp * sinlm - sine * sinp),
        rm * (sine * cosp * sinlm + cose * sinp),
    ];

    SunMoon { sun, moon }
}

/// Analytic Sun and Moon geocentric positions in the Earth-fixed (ITRS) frame
/// (metres) for the given UTC instant.
///
/// The analytic series ([`sun_moon_eci`]) are referred to the **mean equator and
/// equinox of date** (of-date mean longitude and obliquity), so they are rotated
/// to ITRS with [`mean_of_date_to_itrs_matrix`] (nutation + GAST). Precession is
/// already implicit in the of-date series and must NOT be applied a second time;
/// using the full GCRS->ITRS transform here would double-count precession (about
/// 0.37 deg of Sun-direction error at epoch 2026).
pub fn sun_moon_ecef(ts: &TimeScales) -> SunMoon {
    let jd_tt = ts.jd_tt;
    let t = (jd_tt - JD_J2000) / DAYS_PER_CENTURY;

    let eci = sun_moon_eci(t);
    let r = mean_of_date_to_itrs_matrix(ts);
    let sun = mat3_vec3_mul(&r, &eci.sun);
    let moon = mat3_vec3_mul(&r, &eci.moon);
    SunMoon { sun, moon }
}