sidereon-core 0.8.0

//! Precise time scale conversions: UTC -> TAI -> TT -> TDB -> UT1.
//!
//! Mirrors Skyfield's `_utc()` path for bit-exact parity. The delta-T numerics,
//! summation order, and transcendental sequence are preserved EXACTLY: this
//! module is parity-critical and must not be refactored in any way that
//! perturbs a single last bit.
//!
//! The only change from the `orbis_nif` original is visibility: the formerly
//! `pub(crate)` `TimeScales` internals are promoted to a clean public API so a
//! Rust-only consumer of `sidereon-core` can reach the precise time scales
//! without pulling in Rustler or the BEAM.

use crate::astro::constants::time::{
    DAYS_PER_JULIAN_CENTURY, J2000_JD, SECONDS_PER_DAY, TT_MINUS_TAI_S,
};
use crate::astro::data::iers::UT1_DATA;
use crate::astro::time::eop::{
    check_ut1_coverage, CoverageError, LeapSecondTable, TimeScaleInputErrorKind, Ut1Provenance,
    Validated, ValidityMode,
};
use crate::validate::{self, FieldError};

const ROUND_1E7: f64 = 10_000_000.0;

/// Resolved set of Julian-date split time scales for one UTC instant.
///
/// All fields use the Skyfield split convention: `jd_whole` carries the integer
/// (and TAI-aligned) part of the day, and the per-scale `*_fraction` fields carry
/// the residual so that `jd_<scale> == jd_whole + <scale>_fraction` reproduces
/// the full Julian date without catastrophic cancellation.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct TimeScales {
    /// Integer Julian day boundary (TAI-aligned), shared by all scales.
    pub jd_whole: f64,
    /// UT1 day fraction relative to `jd_whole`.
    pub ut1_fraction: f64,
    /// TT day fraction relative to `jd_whole`.
    pub tt_fraction: f64,
    /// TDB day fraction relative to `jd_whole`.
    pub tdb_fraction: f64,
    /// Full UT1 Julian date.
    pub jd_ut1: f64,
    /// Full TT Julian date.
    pub jd_tt: f64,
    /// Full TDB Julian date.
    pub jd_tdb: f64,
}

struct LeapSecondEntry {
    mjd: i32,
    tai_utc: f64,
}

static LEAP_SECONDS: &[LeapSecondEntry] = &[
    LeapSecondEntry {
        mjd: 41317,
        tai_utc: 10.0,
    },
    LeapSecondEntry {
        mjd: 41499,
        tai_utc: 11.0,
    },
    LeapSecondEntry {
        mjd: 41683,
        tai_utc: 12.0,
    },
    LeapSecondEntry {
        mjd: 42048,
        tai_utc: 13.0,
    },
    LeapSecondEntry {
        mjd: 42413,
        tai_utc: 14.0,
    },
    LeapSecondEntry {
        mjd: 42778,
        tai_utc: 15.0,
    },
    LeapSecondEntry {
        mjd: 43144,
        tai_utc: 16.0,
    },
    LeapSecondEntry {
        mjd: 43509,
        tai_utc: 17.0,
    },
    LeapSecondEntry {
        mjd: 43874,
        tai_utc: 18.0,
    },
    LeapSecondEntry {
        mjd: 44239,
        tai_utc: 19.0,
    },
    LeapSecondEntry {
        mjd: 44786,
        tai_utc: 20.0,
    },
    LeapSecondEntry {
        mjd: 45151,
        tai_utc: 21.0,
    },
    LeapSecondEntry {
        mjd: 45516,
        tai_utc: 22.0,
    },
    LeapSecondEntry {
        mjd: 46247,
        tai_utc: 23.0,
    },
    LeapSecondEntry {
        mjd: 47161,
        tai_utc: 24.0,
    },
    LeapSecondEntry {
        mjd: 47892,
        tai_utc: 25.0,
    },
    LeapSecondEntry {
        mjd: 48257,
        tai_utc: 26.0,
    },
    LeapSecondEntry {
        mjd: 48804,
        tai_utc: 27.0,
    },
    LeapSecondEntry {
        mjd: 49169,
        tai_utc: 28.0,
    },
    LeapSecondEntry {
        mjd: 49534,
        tai_utc: 29.0,
    },
    LeapSecondEntry {
        mjd: 50083,
        tai_utc: 30.0,
    },
    LeapSecondEntry {
        mjd: 50448,
        tai_utc: 31.0,
    },
    LeapSecondEntry {
        mjd: 50813,
        tai_utc: 32.0,
    },
    LeapSecondEntry {
        mjd: 53736,
        tai_utc: 33.0,
    },
    LeapSecondEntry {
        mjd: 54832,
        tai_utc: 34.0,
    },
    LeapSecondEntry {
        mjd: 56109,
        tai_utc: 35.0,
    },
    LeapSecondEntry {
        mjd: 57204,
        tai_utc: 36.0,
    },
    LeapSecondEntry {
        mjd: 57754,
        tai_utc: 37.0,
    },
];

impl TimeScales {
    /// Resolve the split-Julian-date time scales for a UTC calendar instant.
    ///
    /// Validates the public boundary, then runs the exact Skyfield `_utc()` path.
    pub fn from_utc(
        year: i32,
        month: i32,
        day: i32,
        hour: i32,
        minute: i32,
        second: f64,
    ) -> Result<Self, CoverageError> {
        validate::finite(second, "second").map_err(map_time_scale_field_error)?;
        validate::civil_datetime_with_second_policy(
            i64::from(year),
            i64::from(month),
            i64::from(day),
            i64::from(hour),
            i64::from(minute),
            second,
            validate::CivilSecondPolicy::UtcLike,
        )
        .map_err(map_time_scale_field_error)?;
        if second >= 60.0 && !is_positive_leap_second_label(year, month, day, hour, minute) {
            return Err(CoverageError::InvalidInput {
                field: "civil datetime",
                kind: TimeScaleInputErrorKind::InvalidCivilTime,
            });
        }
        Ok(Self::from_utc_unchecked(
            year, month, day, hour, minute, second,
        ))
    }

    /// Exact Skyfield `_utc()` path. The arithmetic order below is load-bearing
    /// for 0-ULP parity and MUST NOT be reordered.
    fn from_utc_unchecked(
        year: i32,
        month: i32,
        day: i32,
        hour: i32,
        minute: i32,
        second: f64,
    ) -> Self {
        let jd_day = julian_day_number(year, month, day);
        let jd1 = jd_day as f64 - 0.5;
        let utc_seconds_of_day = hour as f64 * 3600.0 + minute as f64 * 60.0 + second;
        let leap_lookup_second = if second >= 60.0 { 59.0 } else { second };
        let jd2 =
            (leap_lookup_second + minute as f64 * 60.0 + hour as f64 * 3600.0) / SECONDS_PER_DAY;
        let jd_utc_total = jd1 + jd2;

        let leap_seconds = find_leap_seconds(jd_utc_total);
        let utc_seconds_at_midnight = jd1 * SECONDS_PER_DAY;

        let utc_whole_seconds = utc_seconds_of_day.trunc();
        let utc_subsecond = utc_seconds_of_day.fract();

        // Mirror Skyfield's _utc() path.
        let tai_seconds = utc_seconds_at_midnight + leap_seconds + utc_whole_seconds;
        let jd_whole = (tai_seconds / SECONDS_PER_DAY).floor();
        let tai_fraction =
            (tai_seconds - jd_whole * SECONDS_PER_DAY + utc_subsecond) / SECONDS_PER_DAY;
        let tt_offset_days = TT_MINUS_TAI_S / SECONDS_PER_DAY;

        let tt_fraction = tai_fraction + tt_offset_days;
        let jd_tt = jd_whole + tt_fraction;

        let delta_t = interpolate_delta_t(jd_tt);
        let ut1_fraction = tt_fraction - delta_t / SECONDS_PER_DAY;
        let jd_ut1 = jd_whole + ut1_fraction;

        let t = (jd_whole - J2000_JD + tt_fraction) / DAYS_PER_JULIAN_CENTURY;
        let tdb_minus_tt_seconds = 0.001657 * (628.3076 * t + 6.2401).sin()
            + 0.000022 * (575.3385 * t + 4.2970).sin()
            + 0.000014 * (1256.6152 * t + 6.1969).sin()
            + 0.000005 * (606.9777 * t + 4.0212).sin()
            + 0.000005 * (52.9691 * t + 0.4444).sin()
            + 0.000002 * (21.3299 * t + 5.5431).sin()
            + 0.000010 * t * (628.3076 * t + 4.2490).sin();

        let tdb_fraction = tt_fraction + tdb_minus_tt_seconds / SECONDS_PER_DAY;
        let jd_tdb = jd_whole + tdb_fraction;

        Self {
            jd_whole,
            ut1_fraction,
            tt_fraction,
            tdb_fraction,
            jd_ut1,
            jd_tt,
            jd_tdb,
        }
    }

    /// Coverage-policy-enforced variant of [`TimeScales::from_utc`].
    ///
    /// The numerics are produced by [`TimeScales::from_utc`] **unchanged** (the
    /// delta-T / UT1-UTC lookups still clamp at the embedded table edges exactly
    /// as before, preserving Skyfield 0-ULP parity). The only addition is that
    /// the resulting TT instant is classified against the embedded UT1/EOP
    /// coverage interval under the requested [`ValidityMode`]:
    ///
    /// - [`ValidityMode::Strict`]: an instant outside `[first_jd_tt, last_jd_tt]`
    ///   (where the delta-T table would have silently clamped/extrapolated)
    ///   returns [`CoverageError::OutsideCoverage`]. Nothing degraded is ever
    ///   returned.
    /// - [`ValidityMode::Permissive`]: the clamped value is returned, paired with
    ///   a [`crate::astro::time::eop::DegradeReason`] when the instant fell outside
    ///   coverage. This is the historical (parity) behaviour, now made explicit.
    ///
    /// In-coverage results are bit-identical to [`TimeScales::from_utc`] and are
    /// flagged not-degraded.
    pub fn from_utc_validated(
        year: i32,
        month: i32,
        day: i32,
        hour: i32,
        minute: i32,
        second: f64,
        mode: ValidityMode,
    ) -> Result<Validated<Self>, CoverageError> {
        // Numerics first, exactly as the parity path produces them.
        let scales = Self::from_utc(year, month, day, hour, minute, second)?;
        // Classify the (already-clamped) instant against UT1 coverage. We
        // classify at jd_tt because the delta-T table axis is in TT (see
        // `ut1_coverage`), and jd_tt is independent of the clamped delta-T.
        let prov = ut1_coverage();
        let degraded = check_ut1_coverage(&prov, scales.jd_tt, mode)?;
        Ok(Validated {
            value: scales,
            degraded,
        })
    }
}

pub(crate) fn is_positive_leap_second_label(
    year: i32,
    month: i32,
    day: i32,
    hour: i32,
    minute: i32,
) -> bool {
    if hour != 23 || minute != 59 {
        return false;
    }
    let jd1 = julian_day_number(year, month, day) as f64 - 0.5;
    find_leap_seconds(jd1 + 1.0) > find_leap_seconds(jd1)
}

impl From<&FieldError> for TimeScaleInputErrorKind {
    fn from(error: &FieldError) -> Self {
        match error {
            FieldError::Missing { .. } => Self::Missing,
            FieldError::NonFinite { .. } => Self::NonFinite,
            FieldError::NotPositive { .. } => Self::NotPositive,
            FieldError::Negative { .. } => Self::Negative,
            FieldError::OutOfRange { .. } => Self::OutOfRange,
            FieldError::FloatParse { .. } => Self::FloatParse,
            FieldError::IntParse { .. } => Self::IntParse,
            FieldError::InvalidCivilDate { .. } => Self::InvalidCivilDate,
            FieldError::InvalidCivilTime { .. } => Self::InvalidCivilTime,
        }
    }
}

fn map_time_scale_field_error(error: FieldError) -> CoverageError {
    CoverageError::InvalidInput {
        field: error.field(),
        kind: TimeScaleInputErrorKind::from(&error),
    }
}

/// Civil calendar -> Julian day number (Fliegel-style, integer arithmetic).
pub fn julian_day_number(year: i32, month: i32, day: i32) -> i64 {
    let year = i64::from(year);
    let month = i64::from(month);
    let day = i64::from(day);
    let janfeb = month <= 2;
    let g = year + 4716 - if janfeb { 1 } else { 0 };
    let f = (month + 9) % 12;
    let e = 1461 * g / 4 + day - 1402;
    let j = e + (153 * f + 2) / 5;
    j + 38 - ((g + 184) / 100) * 3 / 4
}

/// TAI-UTC (cumulative leap seconds) for a UTC Julian date, from the embedded
/// IERS leap-second table. The lookup clamps to the table edges exactly as the
/// original `orbis_nif` implementation did (parity-preserving).
pub fn find_leap_seconds(jd_utc: f64) -> f64 {
    let mjd = (jd_utc - 2400000.5) as i32;
    let mut ls = 10.0;
    for entry in LEAP_SECONDS {
        if mjd >= entry.mjd {
            ls = entry.tai_utc;
        } else {
            break;
        }
    }
    ls
}

/// Provenance + coverage descriptor for the embedded leap-second table.
///
/// Exposed so a precision pipeline can interrogate table coverage and apply
/// strict-vs-permissive policy (see [`crate::astro::time::eop`]).
pub fn leap_second_table() -> LeapSecondTable {
    LeapSecondTable {
        source: "IERS Bulletin C (TAI-UTC), bundled in sidereon-core",
        first_mjd: LEAP_SECONDS.first().map(|e| e.mjd).unwrap_or(0),
        last_mjd: LEAP_SECONDS.last().map(|e| e.mjd).unwrap_or(0),
        entries: LEAP_SECONDS.len(),
    }
}

fn interpolate_delta_t(jd_tt: f64) -> f64 {
    // Build delta-T table on first call (matching C++ lazy static pattern).
    use std::sync::LazyLock;

    struct DeltaTRow {
        jd_tt: f64,
        delta_t: f64,
    }

    static TABLE: LazyLock<Vec<DeltaTRow>> = LazyLock::new(|| {
        UT1_DATA
            .iter()
            .map(|entry| {
                let jd_utc = entry.mjd as f64 + 2400000.5;
                let leap_seconds = find_leap_seconds(jd_utc);
                let tt_minus_utc = leap_seconds + TT_MINUS_TAI_S;
                let delta_t = ((tt_minus_utc - entry.ut1_utc) * ROUND_1E7).round() / ROUND_1E7;
                DeltaTRow {
                    jd_tt: jd_utc + tt_minus_utc / SECONDS_PER_DAY,
                    delta_t,
                }
            })
            .collect()
    });

    // Binary search for the bracketing entries.
    match TABLE.binary_search_by(|row| row.jd_tt.partial_cmp(&jd_tt).unwrap()) {
        Ok(i) => TABLE[i].delta_t,
        Err(0) => TABLE[0].delta_t,
        Err(i) if i >= TABLE.len() => TABLE.last().unwrap().delta_t,
        Err(i) => {
            let p1 = &TABLE[i - 1];
            let p2 = &TABLE[i];
            p1.delta_t + (jd_tt - p1.jd_tt) * (p2.delta_t - p1.delta_t) / (p2.jd_tt - p1.jd_tt)
        }
    }
}

/// UT1 coverage interval for the embedded EOP table, in TT Julian dates.
///
/// Outside this interval the delta-T interpolation clamps to the nearest table
/// edge (parity-preserving Skyfield behaviour). Strict-mode callers should treat
/// instants outside `[first_jd_tt, last_jd_tt]` as degraded; see
/// [`crate::astro::time::eop`].
pub fn ut1_coverage() -> Ut1Provenance {
    let first = UT1_DATA.first();
    let last = UT1_DATA.last();
    let to_jd_tt = |mjd: i32| -> f64 {
        let jd_utc = mjd as f64 + 2400000.5;
        let tt_minus_utc = find_leap_seconds(jd_utc) + TT_MINUS_TAI_S;
        jd_utc + tt_minus_utc / SECONDS_PER_DAY
    };
    Ut1Provenance {
        source: "IERS Earth Orientation Parameters (UT1-UTC), bundled",
        first_mjd: first.map(|e| e.mjd).unwrap_or(0),
        last_mjd: last.map(|e| e.mjd).unwrap_or(0),
        first_jd_tt: first.map(|e| to_jd_tt(e.mjd)).unwrap_or(0.0),
        last_jd_tt: last.map(|e| to_jd_tt(e.mjd)).unwrap_or(0.0),
        entries: UT1_DATA.len(),
    }
}

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

    #[test]
    fn julian_day_number_widens_extreme_inputs_before_arithmetic() {
        let _ = julian_day_number(i32::MIN, i32::MAX, i32::MAX);
        let _ = julian_day_number(i32::MAX, i32::MIN, i32::MIN);
    }
}