kshana 0.37.0

// SPDX-License-Identifier: AGPL-3.0-only
//! Two-part (high-precision) Julian dates.
//!
//! A Julian date held in a single `f64` carries ~15–16 significant digits, so near JD
//! 2 451 545 (J2000) the least significant bit is ~50 µs — too coarse for sub-µs timing,
//! phase, and frequency-transfer work. Splitting the date into an integer **day** part and
//! a fractional **frac** part (the SOFA / hifitime two-part convention) keeps the full
//! precision of the fraction regardless of the size of the day count, so differences of
//! nearby epochs are exact to the `f64` floor.

/// Seconds in a day.
const SEC_PER_DAY: f64 = 86_400.0;

/// A Julian date as an integer day plus a fractional remainder in `[0, 1)`.
/// `JD = day + frac`, but arithmetic keeps `frac` precise independent of `day`.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct Jd2 {
    /// Integer Julian day (the bulk of the magnitude).
    pub day: f64,
    /// Fractional day in `[0, 1)`.
    pub frac: f64,
}

impl Jd2 {
    /// Split a single-`f64` Julian date into its two-part form.
    pub fn new(jd: f64) -> Self {
        let day = jd.floor();
        Self {
            day,
            frac: jd - day,
        }
    }

    /// Construct from explicit parts, renormalising so `day` is an integer and
    /// `frac ∈ [0, 1)`: a fractional part of `day` is moved into `frac` (the time conversions
    /// read `day` as the integer day). With an integer `day`, as every engine caller passes,
    /// the parts are kept bit for bit.
    pub fn from_parts(day: f64, frac: f64) -> Self {
        let mut j = Self { day, frac };
        j.normalize();
        j
    }

    fn normalize(&mut self) {
        let whole = self.day.floor();
        self.frac += self.day - whole;
        self.day = whole;
        let carry = self.frac.floor();
        self.day += carry;
        self.frac -= carry;
    }

    /// Advance by `seconds` (which may be negative), keeping full precision.
    pub fn add_seconds(self, seconds: f64) -> Self {
        Self::from_parts(self.day, self.frac + seconds / SEC_PER_DAY)
    }

    /// The full Julian date as a single `f64` (loses sub-`f64`-floor precision for large
    /// day counts — use [`Jd2::diff_seconds`] when precision matters).
    pub fn total(self) -> f64 {
        self.day + self.frac
    }

    /// Difference `self − other` in seconds, computed part-by-part so a small difference of
    /// two large dates does not lose precision (the integer days cancel exactly).
    pub fn diff_seconds(self, other: Jd2) -> f64 {
        ((self.day - other.day) + (self.frac - other.frac)) * SEC_PER_DAY
    }
}

/// Gregorian calendar date `(year, month, day)` of the integer Julian day number `jdn`
/// (the Julian date of that day's noon), by the Fliegel and Van Flandern (1968) integer
/// algorithm.
fn calendar_of_jdn(jdn: i64) -> (i32, u32, u32) {
    let l = jdn + 68_569;
    let n = 4 * l / 146_097;
    let l = l - (146_097 * n + 3) / 4;
    let i = 4_000 * (l + 1) / 1_461_001;
    let l = l - 1_461 * i / 4 + 31;
    let j = 80 * l / 2_447;
    let day = l - 2_447 * j / 80;
    let l = j / 11;
    let month = j + 2 - 12 * l;
    let year = 100 * (n - 49) + i + l;
    (year as i32, month as u32, day as u32)
}

/// Integer Julian day number (the Julian date of noon) of a Gregorian calendar date, by the
/// Fliegel and Van Flandern (1968) integer algorithm.
fn jdn_of_calendar(year: i32, month: u32, day: u32) -> i64 {
    let (y, m, d) = (year as i64, month as i64, day as i64);
    let a = (m - 14) / 12;
    (1_461 * (y + 4_800 + a)) / 4 + (367 * (m - 2 - 12 * a)) / 12
        - (3 * ((y + 4_900 + a) / 100)) / 4
        + d
        - 32_075
}

/// TAI − UTC (s) at 00:00 UTC of the day whose noon is Julian day number `jdn` (SOFA
/// `iauDat`: the integer leap-second table from 1972-01-01, the drifting-rate expressions
/// from 1960-01-01, 0 before; [`crate::timescales::tai_minus_utc_on_date`]).
fn tai_minus_utc_on(jdn: i64) -> f64 {
    tai_minus_utc_at(jdn, 0.0)
}

/// TAI − UTC (s) at the fraction `fd` of the UTC day whose noon is Julian day number `jdn`.
fn tai_minus_utc_at(jdn: i64, fd: f64) -> f64 {
    let (y, m, d) = calendar_of_jdn(jdn);
    crate::timescales::tai_minus_utc_on_date(y, m, d, fd)
}

/// The drift of TAI − UTC over the UTC day `jdn` (s/day), `dlod = 2·(ΔAT(12h) − ΔAT(0h))`:
/// zero from 1972; before, the rate at which a UTC second differed from an SI second.
fn utc_day_drift_s(jdn: i64) -> f64 {
    let dat0 = tai_minus_utc_on(jdn);
    2.0 * (tai_minus_utc_at(jdn, 0.5) - dat0)
}

/// Length in UTC seconds of the UTC day whose noon is Julian day number `jdn`: 86 400 s, plus
/// the leap second inserted at its end (86 401 s) or minus one removed, or, before 1972, the
/// fractional step at its end (SOFA `iauDtf2d`: `dleap = ΔAT(next 0h) − (2·ΔAT(12h) − ΔAT(0h))`,
/// which takes the day's own drift out of the step).
fn utc_day_length_s(jdn: i64) -> f64 {
    let dat0 = tai_minus_utc_on(jdn);
    let dat12 = tai_minus_utc_at(jdn, 0.5);
    SEC_PER_DAY + (tai_minus_utc_on(jdn + 1) - (2.0 * dat12 - dat0))
}

/// The UTC day of a two-part quasi Julian date: its Julian day number and the fraction of
/// that day elapsed since 00:00 UTC, in `[0, 1)`.
fn utc_day_and_fraction(jd: Jd2) -> (i64, f64) {
    if jd.frac >= 0.5 {
        (jd.day as i64 + 1, jd.frac - 0.5)
    } else {
        (jd.day as i64, jd.frac + 0.5)
    }
}

impl Jd2 {
    /// A UTC calendar date and time as a two-part UTC quasi Julian date, with the
    /// convention of SOFA (Standards of Fundamental Astronomy) `iauDtf2d` for days that end in
    /// a leap second: such a day is 86 401 s long, `23:59:60.x` is a valid time, and the
    /// fraction of the day is counted in units of its own length. Scope: UTC from 1960-01-01:
    /// before 1972 the day carries the fractional step of the drifting-rate era (SOFA
    /// `iauDat`); before 1960, where UTC is undefined, TAI − UTC is taken as 0, as in SOFA.
    ///
    /// Errors on a month, day, hour or minute out of range, or a second outside `[0, 60)`
    /// (`[0, 61)` in the last minute of a day that ends in a leap second).
    pub fn from_utc_calendar(
        year: i32,
        month: u32,
        day: u32,
        hour: u32,
        minute: u32,
        second: f64,
    ) -> Result<Self, String> {
        if !(1..=12).contains(&month) || day == 0 || day > 31 || hour > 23 || minute > 59 {
            return Err(format!(
                "UTC calendar date out of range: {year}-{month}-{day} {hour}:{minute}"
            ));
        }
        let jdn = jdn_of_calendar(year, month, day);
        if calendar_of_jdn(jdn) != (year, month, day) {
            return Err(format!("no such calendar date: {year}-{month}-{day}"));
        }
        let length = utc_day_length_s(jdn);
        let last_minute = hour == 23 && minute == 59;
        let second_limit = if last_minute {
            60.0 + (length - SEC_PER_DAY)
        } else {
            60.0
        };
        if !(0.0..second_limit).contains(&second) {
            return Err(format!(
                "UTC second {second} out of range at {hour}:{minute}"
            ));
        }
        let elapsed = 60.0 * (60.0 * hour as f64 + minute as f64) + second;
        // 00:00 UTC of the day is Julian date (jdn - 1) + 0.5.
        Ok(Self::from_parts((jdn - 1) as f64, 0.5 + elapsed / length))
    }
}

/// UTC (a two-part quasi Julian date, the [`Jd2::from_utc_calendar`] convention) to TAI
/// (International Atomic Time), the convention of SOFA `iauUtctai`: the fraction of a
/// leap-second day is scaled by that day's 86 401 s length before TAI − UTC is added.
pub fn utc_to_tai(utc: Jd2) -> Jd2 {
    let (jdn, fd) = utc_day_and_fraction(utc);
    let mut elapsed = fd * utc_day_length_s(jdn);
    // Before 1972 a UTC second was not an SI second: scale by (86 400 + dlod)/86 400, as SOFA
    // `iauUtctai` does. From 1972 dlod is zero and this is skipped, so the arithmetic is the
    // integer-leap-second era's unchanged.
    let dlod = utc_day_drift_s(jdn);
    if dlod != 0.0 {
        elapsed *= (SEC_PER_DAY + dlod) / SEC_PER_DAY;
    }
    Jd2::from_parts(
        (jdn - 1) as f64,
        0.5 + (elapsed + tai_minus_utc_on(jdn)) / SEC_PER_DAY,
    )
}

/// TAI to the UTC quasi Julian date, inverting [`utc_to_tai`] by the fixed-point iteration
/// of SOFA `iauTaiutc` (three corrections, each the TAI residual of the current guess).
pub fn tai_to_utc(tai: Jd2) -> Jd2 {
    let mut guess = tai;
    for _ in 0..3 {
        let back = utc_to_tai(guess);
        guess = Jd2::from_parts(
            guess.day + (tai.day - back.day),
            guess.frac + (tai.frac - back.frac),
        );
    }
    guess
}

/// UT1 from the UTC quasi Julian date and UT1 − UTC (s), the convention of SOFA `iauUtcut1`:
/// through TAI, `UT1 = TAI + (UT1 − UTC) − (TAI − UTC)`, with TAI − UTC of the UTC day. Adding
/// UT1 − UTC to the quasi Julian date directly would be wrong by up to a second near the end of
/// a day that ends in a leap second, whose fraction counts 86 401 s.
pub fn utc_to_ut1(utc: Jd2, ut1_minus_utc_s: f64) -> Jd2 {
    let (jdn, _) = utc_day_and_fraction(utc);
    utc_to_tai(utc).add_seconds(ut1_minus_utc_s - tai_minus_utc_on(jdn))
}

/// TT (Terrestrial Time) from TAI: TT = TAI + 32.184 s exactly.
pub fn tai_to_tt(tai: Jd2) -> Jd2 {
    tai.add_seconds(crate::timescales::TT_MINUS_TAI)
}

/// The Earth rotation angle (rad, in `[0, 2π)`) at a two-part UT1 Julian date, the IAU 2000
/// expression of SOFA `iauEra00` evaluated on the two parts so the angle keeps the
/// microsecond a single-`f64` Julian date would lose (about 3e-9 rad of rotation).
pub fn earth_rotation_angle(ut1: Jd2) -> f64 {
    let t = (ut1.day - crate::timescales::JD_J2000) + ut1.frac;
    let turns = ut1.frac + crate::timescales::ERA_TURNS_AT_J2000 + 0.002_737_811_911_354_48 * t;
    std::f64::consts::TAU * turns.rem_euclid(1.0)
}

/// Greenwich mean sidereal time (rad, in `[0, 2π)`), the IAU 1982 model (the SGP4 / TEME
/// convention of [`crate::sgp4::gstime`]) at a two-part UT1 Julian date, evaluated as SOFA
/// `iauGmst82` does on the two parts: the polynomial in Julian centuries from J2000 plus the
/// fraction of the day in seconds, so the sub-day part keeps the microsecond.
pub fn gmst82(ut1: Jd2) -> f64 {
    const A: f64 = 24_110.548_41 - SEC_PER_DAY / 2.0;
    const B: f64 = 8_640_184.812_866;
    const C: f64 = 0.093_104;
    const D: f64 = -6.2e-6;
    const DS2R: f64 = 7.272_205_216_643_04e-5; // 2π / 86 400, ERFA DS2R to f64
    let (d1, d2) = if ut1.day < ut1.frac {
        (ut1.day, ut1.frac)
    } else {
        (ut1.frac, ut1.day)
    };
    let t = (d1 + (d2 - crate::timescales::JD_J2000)) / 36_525.0;
    let f = SEC_PER_DAY * ((d1 % 1.0) + (d2 % 1.0));
    (DS2R * ((A + (B + (C + D * t) * t) * t) + f)).rem_euclid(std::f64::consts::TAU)
}

/// [`crate::timescales::TwoPartJd`] (fraction in `[−½, ½]`) to [`Jd2`] (fraction in `[0, 1)`):
/// the same instant, renormalised. The sum is preserved to the rounding of one addition to the
/// fraction (about 1e-11 s).
impl From<crate::timescales::TwoPartJd> for Jd2 {
    fn from(t: crate::timescales::TwoPartJd) -> Self {
        Jd2::from_parts(t.day, t.frac)
    }
}

/// [`Jd2`] to [`crate::timescales::TwoPartJd`]: the same instant, renormalised to that type's
/// convention by its own constructor.
impl From<Jd2> for crate::timescales::TwoPartJd {
    fn from(j: Jd2) -> Self {
        crate::timescales::TwoPartJd::new(j.day, j.frac)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn gmst82_matches_the_single_f64_gstime_away_from_the_floor() {
        // The same IAU 1982 model; the f64 form differs only by its date rounding.
        for jd in [2_451_545.0, 2_455_000.25, 2_461_326.731] {
            let a = gmst82(Jd2::new(jd));
            let b = crate::sgp4::gstime(jd);
            let d = (a - b + std::f64::consts::PI).rem_euclid(std::f64::consts::TAU)
                - std::f64::consts::PI;
            assert!(d.abs() < 1e-8, "{jd}: {d}");
        }
    }

    #[test]
    fn two_part_types_convert_without_moving_the_instant() {
        let t = crate::timescales::TwoPartJd::new(2_461_326.0, -0.3);
        let j: Jd2 = t.into();
        assert_eq!(j.day, 2_461_325.0);
        assert!(((j.day - t.day) + (j.frac - t.frac)).abs() * 86_400.0 < 1e-10);
        let back: crate::timescales::TwoPartJd = j.into();
        assert!(back.diff_seconds(t).abs() < 1e-10);
    }

    #[test]
    fn round_trips_a_single_f64_date() {
        let jd = 2_451_545.523_4;
        let j = Jd2::new(jd);
        assert_eq!(j.day, 2_451_545.0);
        // The fraction only carries the input's own f64 precision (~1e-9 at JD 2.45e6) —
        // which is exactly the coarseness this two-part form exists to escape.
        assert!((j.frac - 0.523_4).abs() < 1e-9);
        assert!((j.total() - jd).abs() < 1e-9);
    }

    #[test]
    fn from_parts_normalizes_the_fraction() {
        let j = Jd2::from_parts(2_451_545.0, 1.25); // 1.25 days of fraction
        assert_eq!(j.day, 2_451_546.0);
        assert!((j.frac - 0.25).abs() < 1e-15);
        let k = Jd2::from_parts(2_451_545.0, -0.25);
        assert_eq!(k.day, 2_451_544.0);
        assert!((k.frac - 0.75).abs() < 1e-15);
    }

    #[test]
    fn preserves_microsecond_precision_a_single_f64_loses() {
        // One microsecond at J2000.
        let t0 = Jd2::new(2_451_545.0);
        let t1 = t0.add_seconds(1.0e-6);
        // The two-part difference recovers the microsecond exactly…
        assert!(
            (t1.diff_seconds(t0) - 1.0e-6).abs() < 1e-15,
            "Δ = {}",
            t1.diff_seconds(t0)
        );
        // …whereas the naive single-f64 JD round-trip cannot resolve it near 2.45e6.
        let naive = (2_451_545.0_f64 + 1.0e-6 / SEC_PER_DAY) - 2_451_545.0;
        assert!(
            (naive * SEC_PER_DAY - 1.0e-6).abs() > 1.0e-7,
            "single-f64 unexpectedly kept the µs: {}",
            naive * SEC_PER_DAY
        );
    }

    #[test]
    fn add_seconds_is_additive_and_reversible() {
        let t = Jd2::new(2_460_000.0);
        let forward = t.add_seconds(3600.0).add_seconds(-3600.0);
        assert!(forward.diff_seconds(t).abs() < 1e-9);
        // A day of seconds advances the day count by one.
        let plus_day = t.add_seconds(SEC_PER_DAY);
        assert_eq!(plus_day.day, 2_460_001.0);
        assert!(plus_day.frac.abs() < 1e-9);
    }

    #[test]
    fn a_leap_second_day_is_86401_seconds_long_and_tai_steps_once() {
        // 2016-12-31 ended in a leap second: 23:59:60.5 is valid, 00:00:60 is not.
        let a = Jd2::from_utc_calendar(2016, 12, 31, 23, 59, 60.5).unwrap();
        assert!(Jd2::from_utc_calendar(2016, 12, 30, 23, 59, 60.5).is_err());
        assert!(Jd2::from_utc_calendar(2017, 2, 29, 0, 0, 0.0).is_err());
        let b = Jd2::from_utc_calendar(2017, 1, 1, 0, 0, 0.0).unwrap();
        // TAI advances 0.5 s from 23:59:60.5 to the next midnight; TAI - UTC goes 36 -> 37 s.
        let (ta, tb) = (utc_to_tai(a), utc_to_tai(b));
        assert!((tb.diff_seconds(ta) - 0.5).abs() < 1e-9);
        let noon = Jd2::from_utc_calendar(2017, 1, 1, 12, 0, 0.0).unwrap();
        assert_eq!((noon.day, noon.frac), (2_457_755.0, 0.0));
        assert!((utc_to_tai(noon).diff_seconds(noon) - 37.0).abs() < 1e-9);
    }

    #[test]
    fn tai_to_utc_inverts_utc_to_tai_to_the_microsecond_across_a_leap_second() {
        for s in [59.0, 59.999_999, 60.0, 60.5, 60.999_999] {
            let u = Jd2::from_utc_calendar(2015, 6, 30, 23, 59, s).unwrap();
            let back = tai_to_utc(utc_to_tai(u));
            assert!(
                back.diff_seconds(u).abs() < 1e-6,
                "{s}: {}",
                back.diff_seconds(u)
            );
        }
    }

    #[test]
    fn the_two_part_earth_rotation_angle_matches_the_single_f64_one() {
        let ut1 = Jd2::from_parts(2_460_000.0, 0.273_4);
        let a = earth_rotation_angle(ut1);
        let b = crate::timescales::earth_rotation_angle(ut1.total());
        let d =
            (a - b + std::f64::consts::PI).rem_euclid(std::f64::consts::TAU) - std::f64::consts::PI;
        assert!(d.abs() < 1e-8, "{d}");
    }

    #[test]
    fn ut1_runs_through_tai_across_a_leap_second_day() {
        // Across the 2016-12-31 leap second UT1 - UTC steps from about -0.6 s to +0.4 s, and
        // UT1 runs on with the Earth: 1.5 s of TAI separate 23:59:59.5 from the next UTC
        // midnight, and so 1.5 s of UT1 (the quasi Julian date would be off by about 1 s).
        let utc = Jd2::from_utc_calendar(2016, 12, 31, 23, 59, 59.5).unwrap();
        let ut1 = utc_to_ut1(utc, -0.6);
        let midnight = Jd2::from_utc_calendar(2017, 1, 1, 0, 0, 0.0).unwrap();
        let ut1_mid = utc_to_ut1(midnight, 0.4);
        assert!(
            (ut1_mid.diff_seconds(ut1) - 1.5).abs() < 1e-9,
            "{}",
            ut1_mid.diff_seconds(ut1)
        );
        // On an ordinary day it is the plain sum.
        let d = Jd2::from_utc_calendar(2020, 5, 5, 12, 0, 0.0).unwrap();
        assert!((utc_to_ut1(d, 0.25).diff_seconds(d) - 0.25).abs() < 1e-9);
    }

    #[test]
    fn a_fractional_day_part_is_moved_into_the_fraction() {
        // from_parts used to keep a fractional `day`, which utc_to_tai and the Earth rotation
        // read as an integer day, dropping the fraction (15 deg of rotation for one hour).
        for jd in [
            2_461_199.041_666_666_7,
            2_460_310.5,
            2_461_199.25,
            2_461_199.499,
            2_461_199.75,
            2_451_545.0,
        ] {
            let split = Jd2::from_parts(jd, 0.0);
            assert_eq!(split, Jd2::new(jd), "{jd}");
            assert_eq!(utc_to_tai(split), utc_to_tai(Jd2::new(jd)), "{jd}");
            assert_eq!(
                crate::sgp4::teme_to_itrs_matrix(split),
                crate::sgp4::teme_to_itrs_matrix(Jd2::new(jd)),
                "{jd}"
            );
        }
        // A split across both parts lands on the same instant, normalised.
        let j = Jd2::from_parts(2_460_310.5, 0.75);
        assert_eq!((j.day, j.frac), (2_460_311.0, 0.25));
        // Integer day parts are untouched bit for bit.
        let k = Jd2::from_parts(2_433_281.0, 0.5 + 26_157.610_558_16);
        assert_eq!(k.day.fract(), 0.0);
        assert_eq!(k.total(), 2_433_281.0 + (0.5 + 26_157.610_558_16));
    }
}