fits-well 0.1.4

//! Typed time coordinates (§9).
//!
//! Covers the computational core: ISO-8601 datetimes ↔ Julian Date / MJD
//! (proleptic-Gregorian calendar math), `J`/`B` epochs → JD, time-scale
//! conversions ([`TimeScale`]) among `UTC`/`TAI`/`TT`/`GPS`/`TCG`/`TDB`/`TCB`
//! (UTC↔TAI via an embedded leap-second table; TDB via the standard periodic
//! approximation), and a [`FitsTime`] view over a header's time keywords
//! (`TIMESYS`, `MJDREF*`/`JDREF*`/`DATEREF`, `TIMEUNIT`, `TREFPOS`, and the global
//! `DATE-OBS`/`MJD-OBS`/`TSTART`/… set). Validated against `astropy.time`.

use crate::error::FitsError;
use crate::error::Result;
use crate::header::Header;
use crate::keyword::key;

/// JD of the MJD zero point (1858-11-17T00:00 UTC).
const MJD0: f64 = 2_400_000.5;
/// 1977-01-01T00:00:00 TAI — the defining epoch where TT/TCG/TCB/TDB coincide.
const T1977_JD: f64 = 2_443_144.5;
/// TT − TAI, seconds (exact, by definition).
const TT_TAI: f64 = 32.184;
/// GPS = TAI − 19 s (GPS time is offset from TAI by a fixed 19 s).
const TAI_GPS: f64 = 19.0;
/// TCG rate: `(TCG − TT) = L_G · (TT − 1977.0)` (IAU 2000 Resolution B1.9).
const L_G: f64 = 6.969_290_134e-10;
/// TCB rate: `(TCB − TDB) ≈ L_B · (TDB − 1977.0)` (IAU 2006 Resolution B3).
const L_B: f64 = 1.550_519_768e-8;
/// `TDB_0`: the constant term of the `TDB − TT` relation (IAU 2006 Resolution B3).
const TDB_0: f64 = -6.55e-5;
const SEC_PER_DAY: f64 = 86_400.0;

/// A calendar datetime (proleptic Gregorian, UTC-agnostic). `second` may reach
/// 60.x to represent a leap second; the JD arithmetic rolls it over naturally.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Datetime {
    pub year: i64,
    pub month: u32,
    pub day: u32,
    pub hour: u32,
    pub minute: u32,
    pub second: f64,
}

impl Datetime {
    /// Parse an ISO-8601 `FITS` datetime: `YYYY-MM-DD` or
    /// `YYYY-MM-DDThh:mm:ss[.sss…]` (§9.1.1). No component defaulting; the date is
    /// required, the time part optional (then midnight).
    pub fn parse(s: &str) -> Result<Datetime> {
        let invalid = || FitsError::InvalidValue {
            card: format!("DATE '{s}'"),
        };
        let s = s.trim();
        // §9.1.1: no timezone designator is permitted (`Z` or a numeric offset).
        if s.contains(['Z', 'z']) {
            return Err(invalid());
        }
        let (date, time) = match s.split_once('T') {
            Some((d, t)) => (d, Some(t)),
            None => (s, None),
        };
        // `[±]CCYY-MM-DD`: the year has ≥4 digits and an optional sign; month/day
        // are exactly two digits (§9.1.1 — leading zeros may not be omitted).
        let (sign, rest) = match date.strip_prefix('-') {
            Some(r) => (-1, r),
            None => (1, date.strip_prefix('+').unwrap_or(date)),
        };
        let mut dp = rest.split('-');
        let y_str = dp.next().ok_or_else(invalid)?;
        let m_str = dp.next().ok_or_else(invalid)?;
        let d_str = dp.next().ok_or_else(invalid)?;
        if dp.next().is_some() || y_str.len() < 4 || !all_digits(y_str) {
            return Err(invalid());
        }
        let year = sign * y_str.parse::<i64>().map_err(|_| invalid())?;
        // `to_jd` computes `365·(year + 4800)` in `i64`; bound the year so a hostile
        // DATE can't overflow that. Any real observation/epoch date is far inside ±10⁶.
        if year.unsigned_abs() > 1_000_000 {
            return Err(invalid());
        }
        let month = parse_fixed(m_str, 2).ok_or_else(invalid)?;
        let day = parse_fixed(d_str, 2).ok_or_else(invalid)?;
        if !(1..=12).contains(&month) || day < 1 || day > days_in_month(year, month) {
            return Err(invalid());
        }

        let (mut hour, mut minute, mut second) = (0u32, 0u32, 0.0f64);
        if let Some(t) = time {
            let mut tp = t.split(':');
            hour = parse_fixed(tp.next().ok_or_else(invalid)?, 2).ok_or_else(invalid)?;
            minute = parse_fixed(tp.next().ok_or_else(invalid)?, 2).ok_or_else(invalid)?;
            if let Some(sec) = tp.next() {
                second = parse_seconds(sec).ok_or_else(invalid)?;
            }
            // Second 60 is the leap second; this type is scale-agnostic, so the
            // "only in UTC" restriction is left to the caller.
            if tp.next().is_some() || hour >= 24 || minute >= 60 || !(0.0..61.0).contains(&second) {
                return Err(invalid());
            }
        }
        Ok(Datetime {
            year,
            month,
            day,
            hour,
            minute,
            second,
        })
    }

    /// Julian Date of this datetime, interpreting the fields in their own time
    /// scale (no scale conversion is applied here).
    pub fn to_jd(&self) -> f64 {
        let day_fraction =
            (self.hour as f64 * 3600.0 + self.minute as f64 * 60.0 + self.second) / SEC_PER_DAY;
        gregorian_to_jdn(self.year, self.month as i64, self.day as i64) as f64 - 0.5 + day_fraction
    }

    /// Modified Julian Date (`JD − 2400000.5`).
    pub fn to_mjd(&self) -> f64 {
        self.to_jd() - MJD0
    }

    /// Build a datetime from a JD (inverse of [`Datetime::to_jd`]). A single `f64`
    /// JD at present epochs (~2.46e6) resolves to ~0.1 ms, so the recovered second
    /// carries that much rounding — fine for display, not for sub-ms timing.
    pub fn from_jd(jd: f64) -> Datetime {
        // Split into integer day (JDN at noon) and the fraction past midnight.
        let z = (jd + 0.5).floor();
        let frac = (jd + 0.5) - z;
        let date = jdn_to_gregorian(z as i64);
        let mut secs = frac * SEC_PER_DAY;
        let hour = (secs / 3600.0).floor();
        secs -= hour * 3600.0;
        let minute = (secs / 60.0).floor();
        secs -= minute * 60.0;
        Datetime {
            year: date.year,
            month: date.month,
            day: date.day,
            hour: hour as u32,
            minute: minute as u32,
            second: secs,
        }
    }
}

/// True if `s` is non-empty and all ASCII digits.
fn all_digits(s: &str) -> bool {
    !s.is_empty() && s.bytes().all(|b| b.is_ascii_digit())
}

/// Parse a fixed-width all-digit field (§9.1.1 forbids omitted leading zeros, so
/// the length must be exact).
fn parse_fixed(s: &str, width: usize) -> Option<u32> {
    (s.len() == width && all_digits(s))
        .then(|| s.parse().ok())
        .flatten()
}

/// Parse a `ss[.s…]` seconds field: exactly two integer digits, optional fraction.
fn parse_seconds(s: &str) -> Option<f64> {
    let (int, frac) = s.split_once('.').map_or((s, None), |(i, f)| (i, Some(f)));
    if int.len() != 2 || !all_digits(int) || frac.is_some_and(|f| !all_digits(f)) {
        return None;
    }
    s.parse().ok()
}

/// A reference epoch: Julian (`J2000.0`) or Besselian (`B1950.0`).
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum Epoch {
    /// Julian epoch in years (e.g. `2000.0`).
    Julian(f64),
    /// Besselian epoch in years (e.g. `1950.0`).
    Besselian(f64),
}

impl Epoch {
    /// Parse `J<year>` or `B<year>` (e.g. `'J2000.0'`, `'B1950.0'`).
    pub fn parse(s: &str) -> Result<Epoch> {
        let s = s.trim();
        let invalid = || FitsError::InvalidValue {
            card: format!("epoch '{s}'"),
        };
        let (tag, rest) = s.split_at(
            s.char_indices()
                .next()
                .map(|(_, c)| c.len_utf8())
                .unwrap_or(0),
        );
        let year: f64 = rest.parse().map_err(|_| invalid())?;
        match tag {
            "J" | "j" => Ok(Epoch::Julian(year)),
            "B" | "b" => Ok(Epoch::Besselian(year)),
            _ => Err(invalid()),
        }
    }

    /// Julian Date of the epoch. Julian: `2451545.0 + (y−2000)·365.25`; Besselian:
    /// `2415020.31352 + (y−1900)·365.242198781` (the tropical year).
    pub fn to_jd(self) -> f64 {
        match self {
            Epoch::Julian(y) => 2_451_545.0 + (y - 2000.0) * 365.25,
            Epoch::Besselian(y) => 2_415_020.313_52 + (y - 1900.0) * 365.242_198_781,
        }
    }

    /// Modified Julian Date of the epoch.
    pub fn to_mjd(self) -> f64 {
        self.to_jd() - MJD0
    }
}

/// A FITS time scale (`TIMESYS` / `CTYPEi`).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum TimeScale {
    /// Coordinated Universal Time.
    Utc,
    /// Universal Time (UT1) — treated as UTC here (ΔUT1 needs an external table).
    Ut1,
    /// International Atomic Time.
    Tai,
    /// Terrestrial Time.
    Tt,
    /// Geocentric Coordinate Time.
    Tcg,
    /// Barycentric Dynamical Time.
    Tdb,
    /// Barycentric Coordinate Time.
    Tcb,
    /// GPS time.
    Gps,
    /// An unspecified local clock (`LOCAL`); no conversion is possible.
    Local,
}

impl TimeScale {
    /// Parse a `TIMESYS`/`CTYPE` time-scale string (case-insensitive); accepts the
    /// `TDT`/`ET` → `TT` and `IAT` → `TAI` aliases. Unknown values map to `LOCAL`.
    pub fn parse(s: &str) -> TimeScale {
        // A high-precision value appends a realization in parentheses — `TT(TAI)`,
        // `UTC(NIST)` (§9.2.1); strip it before matching the scale name.
        let base = s.trim().split('(').next().unwrap_or("").trim();
        match base.to_ascii_uppercase().as_str() {
            "UTC" | "GMT" => TimeScale::Utc, // §9.2.1: GMT is continuous with UTC
            "UT1" | "UT" => TimeScale::Ut1,
            "TAI" | "IAT" => TimeScale::Tai,
            "TT" | "TDT" | "ET" => TimeScale::Tt,
            "TCG" => TimeScale::Tcg,
            "TDB" => TimeScale::Tdb,
            "TCB" => TimeScale::Tcb,
            "GPS" => TimeScale::Gps,
            _ => TimeScale::Local,
        }
    }

    /// Convert a Julian Date in this scale to a JD in `target`, treating `UT1` as
    /// `UTC` (ΔUT1 = 0). Use [`TimeScale::convert_dut1`] to supply a real ΔUT1.
    pub fn convert(self, jd: f64, target: TimeScale) -> f64 {
        self.convert_dut1(jd, target, 0.0)
    }

    /// Convert a Julian Date with an explicit `dut1 = UT1 − UTC` (seconds, from
    /// IERS) so conversions to/from `UT1` are exact. `Local` passes through.
    pub fn convert_dut1(self, jd: f64, target: TimeScale, dut1: f64) -> f64 {
        if self == target || self == TimeScale::Local || target == TimeScale::Local {
            return jd;
        }
        from_tt(self.to_tt(jd, dut1), target, dut1)
    }

    /// This scale's JD expressed as TT (the common pivot). `dut1 = UT1 − UTC` (s).
    fn to_tt(self, jd: f64, dut1: f64) -> f64 {
        match self {
            TimeScale::Tt => jd,
            TimeScale::Tai => jd + TT_TAI / SEC_PER_DAY,
            TimeScale::Gps => jd + (TT_TAI + TAI_GPS) / SEC_PER_DAY,
            TimeScale::Utc => jd + (TT_TAI + leap_seconds(jd - MJD0)) / SEC_PER_DAY,
            // UT1 → UTC (subtract ΔUT1) → TT.
            TimeScale::Ut1 => {
                let utc = jd - dut1 / SEC_PER_DAY;
                utc + (TT_TAI + leap_seconds(utc - MJD0)) / SEC_PER_DAY
            }
            TimeScale::Tcg => jd - L_G * (jd - T1977_JD),
            TimeScale::Tdb => jd - tdb_minus_tt(jd) / SEC_PER_DAY,
            TimeScale::Tcb => {
                let tdb = jd - L_B * (jd - T1977_JD) + TDB_0 / SEC_PER_DAY;
                tdb - tdb_minus_tt(tdb) / SEC_PER_DAY
            }
            // convert_dut1 returns early for Local, so it never reaches the pivot.
            TimeScale::Local => unreachable!("convert_dut1 short-circuits Local"),
        }
    }
}

/// TT (as a JD) expressed in `target` — the inverse of [`TimeScale::to_tt`].
fn from_tt(tt: f64, target: TimeScale, dut1: f64) -> f64 {
    match target {
        TimeScale::Tt => tt,
        TimeScale::Tai => tt - TT_TAI / SEC_PER_DAY,
        TimeScale::Gps => tt - (TT_TAI + TAI_GPS) / SEC_PER_DAY,
        TimeScale::Utc => {
            let tai = tt - TT_TAI / SEC_PER_DAY;
            // leap is a function of UTC; one lookup at the TAI date suffices away
            // from the ≤1 s boundary ambiguity inherent to UTC.
            tai - leap_seconds(tai - MJD0) / SEC_PER_DAY
        }
        // TT → UTC → UT1 (add ΔUT1).
        TimeScale::Ut1 => {
            let tai = tt - TT_TAI / SEC_PER_DAY;
            let utc = tai - leap_seconds(tai - MJD0) / SEC_PER_DAY;
            utc + dut1 / SEC_PER_DAY
        }
        TimeScale::Tcg => tt + L_G * (tt - T1977_JD),
        TimeScale::Tdb => tt + tdb_minus_tt(tt) / SEC_PER_DAY,
        TimeScale::Tcb => {
            let tdb = tt + tdb_minus_tt(tt) / SEC_PER_DAY;
            tdb + L_B * (tdb - T1977_JD) - TDB_0 / SEC_PER_DAY
        }
        // convert_dut1 returns early for Local, so it never reaches the pivot.
        TimeScale::Local => unreachable!("convert_dut1 short-circuits Local"),
    }
}

/// `TDB − TT` in seconds — the standard periodic approximation (~10 µs accuracy):
/// `0.001658·sin g + 0.000014·sin 2g`, `g = 357.53° + 0.9856003°·(JD_TT − J2000)`.
fn tdb_minus_tt(jd_tt: f64) -> f64 {
    let g = (357.53 + 0.985_600_3 * (jd_tt - 2_451_545.0)).to_radians();
    0.001_658 * g.sin() + 0.000_014 * (2.0 * g).sin()
}

/// `TAI − UTC` in seconds for a given UTC MJD: the integer leap-second count from
/// the IERS table (1972–2017).
///
/// Outside that range the value is **frozen at the nearest table end**, which the
/// caller should be aware of: past the last entry it returns the latest count
/// (correct until a new leap second is announced), and before 1972 it returns the
/// first entry (10 s) even though pre-1972 UTC used fractional "rubber seconds", so
/// UTC↔TAI for pre-1972 dates is approximate, not exact.
fn leap_seconds(mjd: f64) -> f64 {
    // (year, month, day, TAI−UTC) at each step, 1972 onward.
    const TABLE: &[(i64, i64, i64, f64)] = &[
        (1972, 1, 1, 10.0),
        (1972, 7, 1, 11.0),
        (1973, 1, 1, 12.0),
        (1974, 1, 1, 13.0),
        (1975, 1, 1, 14.0),
        (1976, 1, 1, 15.0),
        (1977, 1, 1, 16.0),
        (1978, 1, 1, 17.0),
        (1979, 1, 1, 18.0),
        (1980, 1, 1, 19.0),
        (1981, 7, 1, 20.0),
        (1982, 7, 1, 21.0),
        (1983, 7, 1, 22.0),
        (1985, 7, 1, 23.0),
        (1988, 1, 1, 24.0),
        (1990, 1, 1, 25.0),
        (1991, 1, 1, 26.0),
        (1992, 7, 1, 27.0),
        (1993, 7, 1, 28.0),
        (1994, 7, 1, 29.0),
        (1996, 1, 1, 30.0),
        (1997, 7, 1, 31.0),
        (1999, 1, 1, 32.0),
        (2006, 1, 1, 33.0),
        (2009, 1, 1, 34.0),
        (2012, 7, 1, 35.0),
        (2015, 7, 1, 36.0),
        (2017, 1, 1, 37.0),
    ];
    let mut leap = TABLE[0].3;
    for &(y, m, d, l) in TABLE {
        let threshold = gregorian_to_jdn(y, m, d) as f64 - 0.5 - MJD0;
        if mjd >= threshold {
            leap = l;
        } else {
            break;
        }
    }
    leap
}

/// Julian Day Number at noon of a proleptic-Gregorian calendar date (the standard
/// integer formula).
/// Whether `year` is a leap year in the proleptic Gregorian calendar.
fn is_leap_year(year: i64) -> bool {
    year % 4 == 0 && (year % 100 != 0 || year % 400 == 0)
}

/// Number of days in `month` (1–12) of `year`; `0` for an out-of-range month.
fn days_in_month(year: i64, month: u32) -> u32 {
    match month {
        1 | 3 | 5 | 7 | 8 | 10 | 12 => 31,
        4 | 6 | 9 | 11 => 30,
        2 if is_leap_year(year) => 29,
        2 => 28,
        _ => 0,
    }
}

fn gregorian_to_jdn(year: i64, month: i64, day: i64) -> i64 {
    let a = (14 - month) / 12;
    let y = year + 4800 - a;
    let m = month + 12 * a - 3;
    day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045
}

/// A proleptic-Gregorian calendar date — the result of [`jdn_to_gregorian`].
#[derive(Debug, Clone, Copy, PartialEq)]
struct CalendarDate {
    year: i64,
    month: u32,
    day: u32,
}

/// Inverse of [`gregorian_to_jdn`]: JDN → calendar date.
fn jdn_to_gregorian(jdn: i64) -> CalendarDate {
    let a = jdn + 32044;
    let b = (4 * a + 3) / 146097;
    let c = a - (146097 * b) / 4;
    let d = (4 * c + 3) / 1461;
    let e = c - (1461 * d) / 4;
    let m = (5 * e + 2) / 153;
    let day = e - (153 * m + 2) / 5 + 1;
    let month = m + 3 - 12 * (m / 10);
    let year = 100 * b + d - 4800 + m / 10;
    CalendarDate {
        year,
        month: month as u32,
        day: day as u32,
    }
}

/// A time from a `JEPOCH`/`BEPOCH` keyword: its MJD and the scale the keyword
/// implies (TDB for `JEPOCH`, ET ≈ TT for `BEPOCH`).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct EpochTime {
    pub mjd: f64,
    pub scale: TimeScale,
}

/// The global bound / duration / error time keywords (§9.4, §9.5, §9.7), as read
/// by [`Header::time_bounds`](crate::Header::time_bounds). Start/end are absolute
/// MJD; the rest are in `TIMEUNIT`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct TimeBounds {
    /// Observation start: `MJD-BEG`, else `DATE-BEG` → MJD.
    pub beg_mjd: Option<f64>,
    /// Observation end: `MJD-END`, else `DATE-END` → MJD.
    pub end_mjd: Option<f64>,
    /// Observation midpoint: `MJD-AVG`, else `DATE-AVG` → MJD (§9.5, Table 35).
    pub avg_mjd: Option<f64>,
    /// `XPOSURE` — effective exposure time.
    pub xposure: Option<f64>,
    /// `TELAPSE` — total elapsed time.
    pub telapse: Option<f64>,
    /// `TIMEDEL` — time resolution / bin width.
    pub timedel: Option<f64>,
    /// `TIMEPIXR` — pixel position within a bin (0–1, default 0.5).
    pub timepixr: f64,
    /// `TIMSYER` — systematic time error.
    pub timsyer: Option<f64>,
    /// `TIMRDER` — random time error.
    pub timrder: Option<f64>,
}

/// A Good Time Interval as absolute MJD (§9.7).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GtiInterval {
    pub start_mjd: f64,
    pub stop_mjd: f64,
}

/// The §9.6 `'PHASE'` axis folding parameters: the zero-phase reference time
/// `CZPHSia` and the period `CPERIia`, in `TIMEUNIT` relative to the time
/// reference (the `TSTART` convention).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PhaseAxis {
    pub zero_phase: f64,
    pub period: f64,
}

impl PhaseAxis {
    /// Fold a relative time (in `TIMEUNIT`) to a phase in `[0, 1)`; `0.0` if the
    /// period is zero.
    pub fn fold(&self, time: f64) -> f64 {
        if self.period == 0.0 {
            return 0.0;
        }
        ((time - self.zero_phase) / self.period).rem_euclid(1.0)
    }
}

/// A header's time coordinate frame (§9): the reference epoch, scale, unit, and
/// the resolved global time keywords.
#[derive(Debug, Clone)]
pub struct FitsTime {
    /// `TIMESYS` time scale (default `UTC`).
    pub scale: TimeScale,
    /// Reference epoch as MJD (from `MJDREF`/`MJDREFI`+`MJDREFF`, `JDREF*`, or
    /// `DATEREF`); `0.0` if none is given.
    pub mjdref: f64,
    /// `TIMEUNIT` (default `'s'`).
    pub timeunit: String,
    /// `TIMEOFFS` (§9.4.1): a uniform additive clock correction in `TIMEUNIT`,
    /// equivalent to shifting the reference time. Default `0.0`.
    pub timeoffs: f64,
    /// `TREFPOS` (reference position, e.g. `'TOPOCENTER'`), if present.
    pub trefpos: Option<String>,
}

impl FitsTime {
    /// Parse the time frame from a header. The public entry point is
    /// [`Header::time`](crate::Header::time), which forwards here.
    pub(crate) fn from_header(header: &Header) -> FitsTime {
        let scale = header
            .get_text("TIMESYS")
            .map(TimeScale::parse)
            .unwrap_or(TimeScale::Utc);
        let timeunit = header.get_text("TIMEUNIT").unwrap_or("s").to_string();
        let trefpos = header.get_text("TREFPOS").map(str::to_string);
        FitsTime {
            scale,
            mjdref: reference_mjd(header),
            timeunit,
            timeoffs: header.get_real("TIMEOFFS").unwrap_or(0.0),
            trefpos,
        }
    }

    /// `TIMEUNIT` expressed in seconds (`s`, `d`/day, `a`/`yr` Julian year).
    pub fn unit_seconds(&self) -> f64 {
        // Table 34. The deprecated tropical/Besselian years use their conventional
        // lengths; a truly unknown unit falls back to seconds (the default).
        match self.timeunit.trim() {
            "s" => 1.0,
            "min" => 60.0,
            "h" => 3600.0,
            "d" | "day" => SEC_PER_DAY,
            "a" | "yr" | "y" => 365.25 * SEC_PER_DAY, // Julian year
            "cy" => 36525.0 * SEC_PER_DAY,            // Julian century = 100 a
            "ta" => 365.24219 * SEC_PER_DAY,          // tropical year (deprecated)
            "Ba" => 365.2421988 * SEC_PER_DAY,        // Besselian year (deprecated)
            _ => 1.0,
        }
    }

    /// Resolve a time value measured *relative* to `MJDREF` (e.g. `TSTART`,
    /// `TSTOP`), in `TIMEUNIT`, to an absolute MJD in the frame's own scale. The
    /// `TIMEOFFS` clock correction (§9.4.1) is added before scaling.
    pub fn relative_to_mjd(&self, value: f64) -> f64 {
        self.mjdref + (value + self.timeoffs) * self.unit_seconds() / SEC_PER_DAY
    }

    /// The observation MJD from `MJD-OBS`, else `DATE-OBS`, else `None`. Reads only
    /// the header (not the parsed frame). The public entry point is
    /// [`Header::obs_mjd`](crate::Header::obs_mjd), which forwards here.
    pub(crate) fn obs_mjd(header: &Header) -> Option<f64> {
        if let Some(mjd) = header.get_real("MJD-OBS") {
            return Some(mjd);
        }
        if let Some(d) = header
            .get_text("DATE-OBS")
            .and_then(|s| Datetime::parse(s).ok())
        {
            return Some(d.to_mjd());
        }
        // §9.5: with no DATE-OBS/MJD-OBS, JEPOCH/BEPOCH stand in for the observation time.
        Self::epoch(header).map(|e| e.mjd)
    }

    /// The Julian (`JEPOCH`, implied scale TDB) or Besselian (`BEPOCH`, implied
    /// scale ET ≈ TT) epoch keyword as an [`EpochTime`], if present (§9.1.2, §9.5).
    /// `JEPOCH` wins if both appear. Reads only the header, so it takes no `self`.
    pub(crate) fn epoch(header: &Header) -> Option<EpochTime> {
        if let Some(j) = header.get_real("JEPOCH") {
            return Some(EpochTime {
                mjd: Epoch::Julian(j).to_mjd(),
                scale: TimeScale::Tdb,
            });
        }
        let b = header.get_real("BEPOCH")?;
        Some(EpochTime {
            mjd: Epoch::Besselian(b).to_mjd(),
            scale: TimeScale::Tt, // ET ≈ TT
        })
    }

    /// The global bound / duration / error time keywords (§9.4, §9.5, §9.7). The
    /// start/end are resolved to absolute MJD (`MJD-BEG`/`-END`, else `DATE-BEG`/
    /// `-END`); durations and errors are returned as stored, in `TIMEUNIT`. Reads
    /// only the header, so it takes no `self`.
    pub(crate) fn bounds(header: &Header) -> TimeBounds {
        let mjd_or_date = |mjd: &str, date: &str| {
            header.get_real(mjd).or_else(|| {
                header
                    .get_text(date)
                    .and_then(|s| Datetime::parse(s).ok())
                    .map(|d| d.to_mjd())
            })
        };
        TimeBounds {
            beg_mjd: mjd_or_date("MJD-BEG", "DATE-BEG"),
            end_mjd: mjd_or_date("MJD-END", "DATE-END"),
            avg_mjd: mjd_or_date("MJD-AVG", "DATE-AVG"),
            xposure: header.get_real("XPOSURE"),
            telapse: header.get_real("TELAPSE"),
            timedel: header.get_real("TIMEDEL"),
            timepixr: header.get_real("TIMEPIXR").unwrap_or(0.5), // §9.4.2 default
            timsyer: header.get_real("TIMSYER"),
            timrder: header.get_real("TIMRDER"),
        }
    }

    /// Convert Good Time Interval `START`/`STOP` column values (relative to
    /// `MJDREF`, in `TIMEUNIT`) to absolute-MJD intervals (§9.7). Pairs are taken
    /// element-wise up to the shorter slice.
    pub fn gti_intervals(&self, starts: &[f64], stops: &[f64]) -> Vec<GtiInterval> {
        starts
            .iter()
            .zip(stops)
            .map(|(&s, &e)| GtiInterval {
                start_mjd: self.relative_to_mjd(s),
                stop_mjd: self.relative_to_mjd(e),
            })
            .collect()
    }

    /// If WCS axis `axis` (1-based) is a time axis (`CTYPEi = 'TIME'` or a
    /// time-scale name, §9.2.3), convert a 1-based pixel coordinate along it to an
    /// absolute MJD in the frame's scale: the linear axis value (elapsed time in
    /// `TIMEUNIT` from `MJDREF`) plus the reference. `None` if not a time axis.
    pub fn time_axis_mjd(&self, header: &Header, axis: usize, pixel: f64) -> Option<f64> {
        let ctype = header.get_text(key!("CTYPE{axis}").as_str())?;
        if !is_time_ctype(ctype) {
            return None;
        }
        let crpix = header.get_real(key!("CRPIX{axis}").as_str()).unwrap_or(0.0);
        let crval = header.get_real(key!("CRVAL{axis}").as_str()).unwrap_or(0.0);
        let cdelt = header.get_real(key!("CDELT{axis}").as_str()).unwrap_or(1.0);
        Some(self.relative_to_mjd(crval + cdelt * (pixel - crpix)))
    }

    /// The §9.6 `'PHASE'` axis parameters (`CZPHSia` zero-phase time, `CPERIia`
    /// period) for WCS axis `axis` (1-based), or `None` if it is not a phase axis.
    /// Reads only the header, so it takes no `self`.
    pub(crate) fn phase_axis(header: &Header, axis: usize) -> Option<PhaseAxis> {
        let ctype = header.get_text(key!("CTYPE{axis}").as_str())?;
        if TimeAxisKind::from_ctype(ctype) != Some(TimeAxisKind::Phase) {
            return None;
        }
        Some(PhaseAxis {
            zero_phase: header.get_real(key!("CZPHS{axis}").as_str()).unwrap_or(0.0),
            period: header.get_real(key!("CPERI{axis}").as_str()).unwrap_or(0.0),
        })
    }
}

/// A time-related WCS axis type (§9.6).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum TimeAxisKind {
    /// `'TIME'` or a time-scale name — an absolute time axis (→ MJD).
    Time,
    /// `'PHASE'` — phase folded on a period (`CPERIia`, zero `CZPHSia`).
    Phase,
    /// `'TIMELAG'` — a correlation/cross-spectral time lag.
    Timelag,
    /// `'FREQUENCY'` — a frequency axis.
    Frequency,
}

impl TimeAxisKind {
    /// Classify a `CTYPE` as a time-related axis (§9.6), or `None` if it is not one.
    pub fn from_ctype(ctype: &str) -> Option<TimeAxisKind> {
        let head = ctype.split('-').next().unwrap_or("").trim();
        match head {
            "TIME" => Some(TimeAxisKind::Time),
            "PHASE" => Some(TimeAxisKind::Phase),
            "TIMELAG" => Some(TimeAxisKind::Timelag),
            "FREQUENCY" => Some(TimeAxisKind::Frequency),
            _ if !matches!(TimeScale::parse(head), TimeScale::Local) => Some(TimeAxisKind::Time),
            _ => None,
        }
    }
}

/// True if a `CTYPE` denotes an absolute *time* axis (`'TIME'` or a time-scale
/// name) — the kind [`FitsTime::time_axis_mjd`] resolves to an MJD.
fn is_time_ctype(ctype: &str) -> bool {
    TimeAxisKind::from_ctype(ctype) == Some(TimeAxisKind::Time)
}

/// The reference epoch as MJD: `MJDREF` (or `MJDREFI`+`MJDREFF`), else `JDREF`
/// (or `JDREFI`+`JDREFF`), else `DATEREF`, else `0.0`.
fn reference_mjd(header: &Header) -> f64 {
    if let Some(mjd) = resolve_split_ref(header, "MJDREF", "MJDREFI", "MJDREFF") {
        return mjd;
    }
    if let Some(jd) = resolve_split_ref(header, "JDREF", "JDREFI", "JDREFF") {
        return jd - MJD0;
    }
    header
        .get_text("DATEREF")
        .and_then(|s| Datetime::parse(s).ok())
        .map(|d| d.to_mjd())
        .unwrap_or(0.0)
}

/// Resolve a reference epoch from its single (`MJDREF`) and split-precision
/// (`MJDREFI`+`MJDREFF`) keywords. Per §9.2.2 a *full* integer+fractional split
/// takes precedence over the single value; otherwise the single value is used,
/// falling back to a lone split part.
fn resolve_split_ref(header: &Header, single: &str, int: &str, frac: &str) -> Option<f64> {
    let i = header.get_real(int);
    let f = header.get_real(frac);
    match (i, f) {
        (Some(i), Some(f)) => Some(i + f),
        _ => header.get_real(single).or_else(|| match (i, f) {
            (None, None) => None,
            _ => Some(i.unwrap_or(0.0) + f.unwrap_or(0.0)),
        }),
    }
}

#[cfg(test)]
mod tests;