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

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

//! Pluto ephemeris computation using the Meeus/Williams abbreviated series.
//!
//! This module provides a single public function [`Pluto::get_heliocentric`] that
//! returns Pluto's heliocentric ecliptic position for a given Julian Day.
//!
//! The algorithm is taken from Jean Meeus, *Astronomical Algorithms*, 2nd ed.
//! (1998), chapter 36, which in turn reproduces the numerical series derived
//! by T. G. Williams (1991). Accuracy is on the order of **0.5 arc‑seconds in
//! longitude** and **0.2 arc‑seconds in latitude** over the interval 1885 – 2099,
//! while being dramatically faster than a full numerical integration or the
//! DE ephemerides.
//!
//! Coefficient tables live in the crate-private `archive::pluto_tables` module.
//!
//! ```text
//! Step outline (see code for details):
//!  1. Convert the supplied Julian Day to Julian centuries T w.r.t. J2000.0.
//!  2. Compute mean longitudes λJ, λS, λP of Jupiter, Saturn and Pluto using
//!     linear expressions in T.
//!  3. For each of the 42 / 43 periodic terms:
//!       A. Form the argument  Ai = j·λJ + s·λS + p·λP  (degrees).
//!       B. Evaluate sin Ai and cos Ai.
//!       C. Accumulate contributions to ΣL, ΣB, ΣR via  A·sin Ai + B·cos Ai.
//!  4. Scale the sums and add constant/base terms to obtain L, B, R.
//!  5. Convert spherical ⟨L,B,R⟩ to Cartesian ⟨x,y,z⟩ in the ecliptic frame.
//! ```
//!
//! # References
//! * Meeus, J. (1998). *Astronomical Algorithms* (2nd ed.). Willmann‑Bell.
//! * Williams, T. G. (1991). "An optimized algorithm for Pluto", *Mem. Brit.
//!   Astron. Assoc.* **99** (2), 75–82.

use crate::archive::pluto_tables::{ARGUMENTS, LATITUDE_TERMS, LONGITUDE_TERMS, RADIUS_TERMS};
use crate::coordinates::{cartesian, centers::Heliocentric, frames::EclipticMeanJ2000, spherical};
use crate::qtty::{AstronomicalUnit, Degrees, Radian, AU};
use crate::time::JulianDate;

pub struct Pluto;

impl Pluto {
    pub fn get_heliocentric(
        jd: JulianDate,
    ) -> cartesian::Position<Heliocentric, EclipticMeanJ2000, AstronomicalUnit> {
        let t = jd.julian_centuries();

        // 2. Calculate mean longitudes (in degrees) for Jupiter, Saturn, and Pluto.
        let jupiter_lon = Degrees::new(34.35 + 3034.9057 * t);
        let saturn_lon = Degrees::new(50.08 + 1222.1138 * t);
        let pluto_lon = Degrees::new(238.96 + 144.9600 * t);

        // 3. Initialize sums for the periodic terms.
        let mut sum_longitude = Degrees::new(0.0);
        let mut sum_latitude = Degrees::new(0.0);
        let mut sum_radius = 0.0;

        // 4. Loop over all periodic terms.
        for i in 0..ARGUMENTS.len() {
            // Calculate the argument:
            let a = jupiter_lon * ARGUMENTS[i].j
                + saturn_lon * ARGUMENTS[i].s
                + pluto_lon * ARGUMENTS[i].p;

            // Convert 'a' from degrees to radians.
            let (sin_a, cos_a) = a.to::<Radian>().sin_cos();

            // Add periodic corrections for longitude, latitude, and radius.
            sum_longitude +=
                Degrees::new(LONGITUDE_TERMS[i].a * sin_a + LONGITUDE_TERMS[i].b * cos_a);
            sum_latitude += Degrees::new(LATITUDE_TERMS[i].a * sin_a + LATITUDE_TERMS[i].b * cos_a);
            sum_radius += RADIUS_TERMS[i].a * sin_a + RADIUS_TERMS[i].b * cos_a;
        }

        // 5. Calculate the final heliocentric spherical coordinates.
        // These base values and scale factors come from the chosen model (e.g., Meeus's data).
        let lon = sum_longitude * 0.000001 + Degrees::new(238.958116 + 144.96 * t);
        let lat = sum_latitude * 0.000001 - Degrees::new(3.908239);
        let rad = sum_radius * 0.0000001 + 40.7241346;

        // EclipticMeanJ2000: lon = azimuth, lat = polar
        spherical::Position::<Heliocentric, EclipticMeanJ2000, AstronomicalUnit>::new(
            lon,
            lat,
            rad * AU,
        )
        .to_cartesian()
    }
}

#[cfg(test)]
mod tests {
    use crate::calculus::pluto::Pluto;
    use crate::qtty::AstronomicalUnits;
    use crate::time::JulianDate;

    #[test]
    fn pluto_heliocentric_position_j2000() {
        let pos = Pluto::get_heliocentric(JulianDate::J2000);
        assert!(
            (pos.x() - AstronomicalUnits::new(-9.875333629852145))
                .abs()
                .value()
                < 1e-6
        );
        assert!(
            (pos.y() - AstronomicalUnits::new(-27.958786187190157))
                .abs()
                .value()
                < 1e-6
        );
        assert!(
            (pos.z() - AstronomicalUnits::new(5.850444258527083))
                .abs()
                .value()
                < 1e-6
        );
    }
}