siderust 0.6.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

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

//! # Shared Horizontal Coordinate Helpers
//!
//! Common routines for converting equatorial coordinates to horizontal
//! (altitude/azimuth), shared across Sun, Moon, and stellar modules.
//!
//! These helpers factor out the duplicated equatorial→horizontal transform
//! that previously existed independently in `solar::sun_equations` and
//! `lunar::moon_equations`.

use crate::astro::earth_rotation::jd_ut1_from_tt_eop;
use crate::astro::eop::EopProvider;
use crate::astro::nutation::nutation_iau2000b;
use crate::astro::precession;
use crate::astro::sidereal::gast_iau2006;
use crate::coordinates::centers::Geodetic;
use crate::coordinates::frames::ECEF;
use crate::coordinates::transform::AstroContext;
use crate::coordinates::{cartesian, centers::*, frames, spherical};
use crate::time::JulianDate;
use qtty::{AstronomicalUnits, Degree, LengthUnit, Meter, Quantity, Radian, Radians};

// =============================================================================
// Apparent Topocentric Pipeline
// =============================================================================

/// Transforms a **geocentric, mean-J2000 equatorial** Cartesian position into
/// **topocentric, true-of-date equatorial** spherical coordinates.
///
/// This pipeline is shared between Sun and Moon; the only prerequisite is
/// that the caller provides the body's geocentric J2000 Cartesian position
/// (obtained from VSOP87 for the Sun or ELP2000 for the Moon).
///
/// ## Pipeline Steps
///
/// 1. **Topocentric parallax**: Geocentric → Topocentric in J2000 frame
/// 2. **Precession**: J2000 → Mean-of-Date rotation
/// 3. **Nutation**: Mean-of-Date → True-of-Date rotation
/// 4. **Spherical conversion**: Cartesian → Spherical (RA, Dec, distance)
pub fn geocentric_j2000_to_apparent_topocentric<U: LengthUnit>(
    geo_cart_j2000: &cartesian::Position<Geocentric, frames::EquatorialMeanJ2000, U>,
    site: Geodetic<frames::ECEF>,
    jd: JulianDate,
) -> spherical::Position<Topocentric, frames::EquatorialTrueOfDate, U>
where
    Quantity<U>: From<Quantity<Meter>> + From<AstronomicalUnits>,
{
    let ctx: AstroContext = AstroContext::default();
    geocentric_j2000_to_apparent_topocentric_with_ctx(geo_cart_j2000, site, jd, &ctx)
}

/// Context-aware variant of [`geocentric_j2000_to_apparent_topocentric`].
pub fn geocentric_j2000_to_apparent_topocentric_with_ctx<
    U: LengthUnit,
    Eph,
    Eop: EopProvider,
    Nut,
>(
    geo_cart_j2000: &cartesian::Position<Geocentric, frames::EquatorialMeanJ2000, U>,
    site: Geodetic<frames::ECEF>,
    jd: JulianDate,
    ctx: &AstroContext<Eph, Eop, Nut>,
) -> spherical::Position<Topocentric, frames::EquatorialTrueOfDate, U>
where
    Quantity<U>: From<Quantity<Meter>> + From<AstronomicalUnits>,
{
    use crate::coordinates::transform::to_topocentric_with_ctx;

    // 1) Translate geocentric → topocentric in J2000 frame (applies parallax)
    let topo_cart_j2000 = to_topocentric_with_ctx(geo_cart_j2000, site, jd, ctx);

    // 2) Rotate J2000 → mean-of-date using IAU 2006 precession
    let rot_prec = precession::precession_matrix_iau2006(jd);
    let [x_m, y_m, z_m] = rot_prec
        * [
            topo_cart_j2000.x(),
            topo_cart_j2000.y(),
            topo_cart_j2000.z(),
        ];

    let topo_cart_mod =
        cartesian::Position::<Topocentric, frames::EquatorialMeanOfDate, U>::new_with_params(
            *topo_cart_j2000.center_params(),
            x_m,
            y_m,
            z_m,
        );

    // 3) Apply nutation rotation (mean-of-date → true-of-date) — IAU 2000B
    let rot_nut = crate::astro::nutation::nutation_rotation_iau2000b(jd);
    let [x_t, y_t, z_t] = rot_nut * [topo_cart_mod.x(), topo_cart_mod.y(), topo_cart_mod.z()];

    let topo_cart_true =
        cartesian::Position::<Topocentric, frames::EquatorialTrueOfDate, U>::new_with_params(
            *topo_cart_mod.center_params(),
            x_t,
            y_t,
            z_t,
        );

    // 4) Convert to spherical topocentric equatorial (true of date)
    spherical::Position::from_cartesian(&topo_cart_true)
}

// =============================================================================
// Equatorial → Horizontal Transform
// =============================================================================

/// Converts an **apparent topocentric equatorial** position to **horizontal**
/// coordinates (altitude, azimuth, distance).
///
/// This is the standard equatorial-to-horizontal transform using the observer's
/// latitude and the local apparent sidereal time (GAST → LAST → HA).
///
/// Both Sun and Moon `get_horizontal` methods delegate here after computing
/// their body-specific apparent topocentric equatorial coordinates.
pub fn equatorial_to_horizontal<U: LengthUnit>(
    eq_position: &spherical::Position<Topocentric, frames::EquatorialTrueOfDate, U>,
    site: Geodetic<frames::ECEF>,
    jd: JulianDate,
) -> spherical::Position<Topocentric, frames::Horizontal, U> {
    equatorial_to_horizontal_true_of_date(eq_position, site, jd)
}

/// Converts apparent topocentric **true-of-date** equatorial coordinates to
/// horizontal coordinates using IAU 2006 GAST.
pub fn equatorial_to_horizontal_true_of_date<U: LengthUnit>(
    eq_position: &spherical::Position<Topocentric, frames::EquatorialTrueOfDate, U>,
    site: Geodetic<frames::ECEF>,
    jd: JulianDate,
) -> spherical::Position<Topocentric, frames::Horizontal, U> {
    let ctx: AstroContext = AstroContext::default();
    equatorial_to_horizontal_true_of_date_with_ctx(eq_position, site, jd, &ctx)
}

/// Context-aware variant of [`equatorial_to_horizontal_true_of_date`].
///
/// Uses the context EOP provider to derive UT1 for GAST.
pub fn equatorial_to_horizontal_true_of_date_with_ctx<U: LengthUnit, Eph, Eop: EopProvider, Nut>(
    eq_position: &spherical::Position<Topocentric, frames::EquatorialTrueOfDate, U>,
    site: Geodetic<ECEF>,
    jd: JulianDate,
    ctx: &AstroContext<Eph, Eop, Nut>,
) -> spherical::Position<Topocentric, frames::Horizontal, U> {
    // Extract RA, Dec, distance from the equatorial position
    let ra = eq_position.azimuth;
    let dec = eq_position.polar;
    let distance = eq_position.distance;

    // True-of-date RA/Dec must use apparent sidereal time (GAST).
    let eop = ctx.eop_at(jd);
    let jd_ut1 = jd_ut1_from_tt_eop(jd, &eop);
    let nut = nutation_iau2000b(jd);
    let gast = gast_iau2006(jd_ut1, jd, nut.dpsi, nut.true_obliquity());
    let lst = gast + site.lon.to::<Radian>();
    let ha = (lst - ra.to::<Radian>())
        .value()
        .rem_euclid(std::f64::consts::TAU);
    let ha = Radians::new(ha);

    // Convert equatorial to horizontal using standard spherical trig
    let lat = site.lat.to::<Radian>();
    let dec_rad = dec.to::<Radian>();

    // Altitude: sin(alt) = sin(dec)·sin(lat) + cos(dec)·cos(lat)·cos(HA)
    let sin_alt = dec_rad.sin() * lat.sin() + dec_rad.cos() * lat.cos() * ha.cos();
    let alt = Radians::new(sin_alt.asin()).to::<Degree>();

    // Azimuth: tan(az) = -sin(HA) / (cos(lat)·tan(dec) - sin(lat)·cos(HA))
    let az_rad = (-dec_rad.cos() * ha.sin())
        .atan2(dec_rad.sin() * lat.cos() - dec_rad.cos() * ha.cos() * lat.sin());
    let az = Radians::new(az_rad).normalize().to::<Degree>();

    Topocentric::horizontal(site, alt, az, distance)
}

/// Backward-compatible context-aware wrapper.
#[inline]
pub fn equatorial_to_horizontal_with_ctx<U: LengthUnit, Eph, Eop: EopProvider, Nut>(
    eq_position: &spherical::Position<Topocentric, frames::EquatorialTrueOfDate, U>,
    site: Geodetic<frames::ECEF>,
    jd: JulianDate,
    ctx: &AstroContext<Eph, Eop, Nut>,
) -> spherical::Position<Topocentric, frames::Horizontal, U> {
    equatorial_to_horizontal_true_of_date_with_ctx(eq_position, site, jd, ctx)
}

/// Computes the **horizontal direction** (altitude, azimuth) for a fixed-star
/// RA/Dec (J2000) from an observer site at a given Julian Date.
///
/// This is the stellar equivalent of `Sun::get_horizontal` / `Moon::get_horizontal`.
/// Since stars lack a distance, only a `Direction<Horizontal>` is returned.
///
/// ## Pipeline
///
/// 1. Precess J2000 → mean-of-date
/// 2. Apply nutation correction to RA
/// 3. Compute GAST → LAST → HA
/// 4. Standard equatorial→horizontal altitude/azimuth formula
pub fn star_horizontal(
    ra_j2000: qtty::Degrees,
    dec_j2000: qtty::Degrees,
    site: &Geodetic<frames::ECEF>,
    jd: JulianDate,
) -> spherical::Direction<frames::Horizontal> {
    // Full IAU 2006/2000B NPB matrix-based approach:
    // 1. Convert J2000 direction to Cartesian unit vector
    // 2. Apply full precession-nutation-bias (NPB) matrix → true-of-date
    // 3. Convert back to spherical → RA/Dec of date
    // 4. Compute GAST (IAU 2006) → LAST → HA → altitude/azimuth

    let ra_rad = ra_j2000.to::<Radian>();
    let dec_rad = dec_j2000.to::<Radian>();

    // J2000 direction as unit vector
    let (sin_ra, cos_ra) = ra_rad.sin_cos();
    let (sin_dec, cos_dec) = dec_rad.sin_cos();
    let x0 = cos_dec * cos_ra;
    let y0 = cos_dec * sin_ra;
    let z0 = sin_dec;

    // Full NPB matrix: GCRS → true equator/equinox of date
    let nut = nutation_iau2000b(jd);
    let npb = precession::precession_nutation_matrix(jd, nut.dpsi, nut.deps);
    let [x_t, y_t, z_t] = npb.apply_array([x0, y0, z0]);

    // True-of-date RA/Dec
    let ra_tod = Radians::new(y_t.atan2(x_t));
    let dec_tod = Radians::new(z_t.asin());

    // True-of-date RA/Dec requires apparent sidereal time (GAST).
    let ctx: AstroContext = AstroContext::default();
    let eop = ctx.eop_at(jd);
    let jd_ut1 = jd_ut1_from_tt_eop(jd, &eop);
    let gast = gast_iau2006(jd_ut1, jd, nut.dpsi, nut.true_obliquity());
    let lst = gast + site.lon.to::<Radian>();
    let ha = Radians::new((lst - ra_tod).value().rem_euclid(std::f64::consts::TAU));

    // Equatorial → horizontal
    let lat = site.lat.to::<Radian>();
    let sin_alt = dec_tod.sin() * lat.sin() + dec_tod.cos() * lat.cos() * ha.cos();
    let alt = Radians::new(sin_alt.asin()).to::<Degree>();

    let az_rad = (-dec_tod.cos() * ha.sin())
        .atan2(dec_tod.sin() * lat.cos() - dec_tod.cos() * ha.cos() * lat.sin());
    let az = Radians::new(az_rad).normalize().to::<Degree>();

    spherical::Direction::<frames::Horizontal>::new(alt, az)
}