sidereon_core/reduced_orbit/

time.rs

//! Time-scale bridge for reduced-orbit fitting/evaluation.

use crate::astro::constants::time::{SECONDS_PER_DAY, TT_MINUS_TAI_S};
use crate::astro::time::model::TimeScale;
use crate::astro::time::scales::{find_leap_seconds, julian_day_number, TimeScales};

use crate::constants::{BDT_MINUS_TAI_S, GPST_MINUS_TAI_S};

/// A UTC calendar instant `(year, month, day, hour, minute, second)`, the form
/// the core [`TimeScales::from_utc`] consumes. The Elixir layer produces these
/// from each sample/query epoch; no `Instant`->`TimeScales` bridge exists in the
/// core crate, so the calendar tuple is carried explicitly to the boundary.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct CalendarEpoch {
    /// Calendar year.
    pub year: i32,
    /// Calendar month, 1-12.
    pub month: i32,
    /// Calendar day of month, 1-31.
    pub day: i32,
    /// Hour of day, 0-23.
    pub hour: i32,
    /// Minute of hour, 0-59.
    pub minute: i32,
    /// Second of minute, fractional.
    pub second: f64,
}

impl CalendarEpoch {
    /// Construct a calendar epoch from its components.
    pub const fn new(year: i32, month: i32, day: i32, hour: i32, minute: i32, second: f64) -> Self {
        Self {
            year,
            month,
            day,
            hour,
            minute,
            second,
        }
    }

    /// Build the core [`TimeScales`] for this instant, interpreted in `scale`.
    ///
    /// Non-UTC scales (GPST/GST/BDT/TAI/TT) are converted to UTC before the
    /// Skyfield split is built, so the Earth orientation used by the frame
    /// transforms is correct rather than offset by the scale's leap-second gap.
    pub(crate) fn time_scales(self, scale: TimeScale) -> TimeScales {
        let utc = scale_calendar_to_utc(self, scale);
        TimeScales::from_utc(
            utc.year, utc.month, utc.day, utc.hour, utc.minute, utc.second,
        )
        .expect("calendar epoch has a finite second")
    }
}

fn scale_calendar_to_utc(cal: CalendarEpoch, scale: TimeScale) -> CalendarEpoch {
    if matches!(scale, TimeScale::Utc) {
        return cal;
    }
    let tai = normalize_calendar_seconds(cal, cal.second + tai_minus_scale_seconds(scale));
    tai_calendar_to_utc(tai)
}

fn tai_minus_scale_seconds(scale: TimeScale) -> f64 {
    match scale {
        TimeScale::Utc => 0.0,
        TimeScale::Tai => 0.0,
        TimeScale::Tt | TimeScale::Tdb => -TT_MINUS_TAI_S,
        TimeScale::Gpst | TimeScale::Gst => GPST_MINUS_TAI_S,
        TimeScale::Bdt => BDT_MINUS_TAI_S,
    }
}

fn tai_calendar_to_utc(tai: CalendarEpoch) -> CalendarEpoch {
    if let Some(utc) = positive_leap_second_utc_label(tai) {
        return utc;
    }

    let mut leap = leap_seconds_at_utc_label(tai);
    let mut utc = normalize_calendar_seconds(tai, tai.second - leap);
    for _ in 0..3 {
        let next_leap = leap_seconds_at_utc_label(utc);
        if next_leap == leap {
            return utc;
        }
        leap = next_leap;
        utc = normalize_calendar_seconds(tai, tai.second - leap);
    }
    utc
}

fn positive_leap_second_utc_label(tai: CalendarEpoch) -> Option<CalendarEpoch> {
    let tai_sod = seconds_of_day(tai);
    let utc_midnight = CalendarEpoch::new(tai.year, tai.month, tai.day, 0, 0, 0.0);
    let previous_second = normalize_calendar_seconds(utc_midnight, -1.0);
    let old_leap = leap_seconds_at_utc_label(previous_second);
    let new_leap = leap_seconds_at_utc_label(utc_midnight);
    if new_leap <= old_leap || !(old_leap..new_leap).contains(&tai_sod) {
        return None;
    }

    let mut utc = previous_second;
    utc.second = 60.0 + (tai_sod - old_leap);
    Some(utc)
}

fn leap_seconds_at_utc_label(cal: CalendarEpoch) -> f64 {
    let jd1 = julian_day_number(cal.year, cal.month, cal.day) as f64 - 0.5;
    let lookup_second = if cal.second >= 60.0 { 59.0 } else { cal.second };
    let jd2 =
        (cal.hour as f64 * 3600.0 + cal.minute as f64 * 60.0 + lookup_second) / SECONDS_PER_DAY;
    find_leap_seconds(jd1 + jd2)
}

fn seconds_of_day(cal: CalendarEpoch) -> f64 {
    cal.hour as f64 * 3600.0 + cal.minute as f64 * 60.0 + cal.second
}

fn normalize_calendar_seconds(mut cal: CalendarEpoch, second: f64) -> CalendarEpoch {
    cal.second = second;
    while cal.second < 0.0 {
        cal.second += 60.0;
        cal.minute -= 1;
    }
    while cal.second >= 60.0 {
        cal.second -= 60.0;
        cal.minute += 1;
    }
    while cal.minute < 0 {
        cal.minute += 60;
        cal.hour -= 1;
    }
    while cal.minute > 59 {
        cal.minute -= 60;
        cal.hour += 1;
    }
    while cal.hour < 0 {
        cal.hour += 24;
        cal.day -= 1;
    }
    while cal.hour > 23 {
        cal.hour -= 24;
        cal.day += 1;
    }
    while cal.day < 1 {
        cal.month -= 1;
        if cal.month < 1 {
            cal.month = 12;
            cal.year -= 1;
        }
        cal.day = days_in_month(cal.year, cal.month);
    }
    loop {
        let month_days = days_in_month(cal.year, cal.month);
        if cal.day <= month_days {
            break;
        }
        cal.day -= month_days;
        cal.month += 1;
        if cal.month > 12 {
            cal.month = 1;
            cal.year += 1;
        }
    }
    cal
}

fn days_in_month(year: i32, month: i32) -> i32 {
    match month {
        1 | 3 | 5 | 7 | 8 | 10 | 12 => 31,
        4 | 6 | 9 | 11 => 30,
        2 if is_leap_year(year) => 29,
        2 => 28,
        _ => 31,
    }
}

fn is_leap_year(year: i32) -> bool {
    (year % 4 == 0 && year % 100 != 0) || year % 400 == 0
}

/// Seconds between two calendar epochs via their J2000-TT split day numbers.
pub(crate) fn dt_seconds(t0: &TimeScales, t: &TimeScales) -> f64 {
    let dwhole = t.jd_whole - t0.jd_whole;
    let dfrac = t.tt_fraction - t0.tt_fraction;
    (dwhole + dfrac) * SECONDS_PER_DAY
}