fits_well/time/
mod.rs

1//! Typed time coordinates (§9).
2//!
3//! Covers the computational core: ISO-8601 datetimes ↔ Julian Date / MJD
4//! (proleptic-Gregorian calendar math), `J`/`B` epochs → JD, time-scale
5//! conversions ([`TimeScale`]) among `UTC`/`TAI`/`TT`/`GPS`/`TCG`/`TDB`/`TCB`
6//! (UTC↔TAI via an embedded leap-second table; TDB via the standard periodic
7//! approximation), and a [`FitsTime`] view over a header's time keywords
8//! (`TIMESYS`, `MJDREF*`/`JDREF*`/`DATEREF`, `TIMEUNIT`, `TREFPOS`, and the global
9//! `DATE-OBS`/`MJD-OBS`/`TSTART`/… set). Validated against `astropy.time`.
10
11use crate::error::FitsError;
12use crate::error::Result;
13use crate::header::Header;
14use crate::keyword::key;
15
16/// JD of the MJD zero point (1858-11-17T00:00 UTC).
17const MJD0: f64 = 2_400_000.5;
18/// 1977-01-01T00:00:00 TAI — the defining epoch where TT/TCG/TCB/TDB coincide.
19const T1977_JD: f64 = 2_443_144.5;
20/// TT − TAI, seconds (exact, by definition).
21const TT_TAI: f64 = 32.184;
22/// GPS = TAI − 19 s (GPS time is offset from TAI by a fixed 19 s).
23const TAI_GPS: f64 = 19.0;
24/// TCG rate: `(TCG − TT) = L_G · (TT − 1977.0)` (IAU 2000 Resolution B1.9).
25const L_G: f64 = 6.969_290_134e-10;
26/// TCB rate: `(TCB − TDB) ≈ L_B · (TDB − 1977.0)` (IAU 2006 Resolution B3).
27const L_B: f64 = 1.550_519_768e-8;
28/// `TDB_0`: the constant term of the `TDB − TT` relation (IAU 2006 Resolution B3).
29const TDB_0: f64 = -6.55e-5;
30const SEC_PER_DAY: f64 = 86_400.0;
31
32/// A calendar datetime (proleptic Gregorian, UTC-agnostic). `second` may reach
33/// 60.x to represent a leap second; the JD arithmetic rolls it over naturally.
34#[derive(Debug, Clone, Copy, PartialEq)]
35pub struct Datetime {
36    pub year: i64,
37    pub month: u32,
38    pub day: u32,
39    pub hour: u32,
40    pub minute: u32,
41    pub second: f64,
42}
43
44impl Datetime {
45    /// Parse an ISO-8601 `FITS` datetime: `YYYY-MM-DD` or
46    /// `YYYY-MM-DDThh:mm:ss[.sss…]` (§9.1.1). No component defaulting; the date is
47    /// required, the time part optional (then midnight).
48    pub fn parse(s: &str) -> Result<Datetime> {
49        let invalid = || FitsError::InvalidValue {
50            card: format!("DATE '{s}'"),
51        };
52        let s = s.trim();
53        // §9.1.1: no timezone designator is permitted (`Z` or a numeric offset).
54        if s.contains(['Z', 'z']) {
55            return Err(invalid());
56        }
57        let (date, time) = match s.split_once('T') {
58            Some((d, t)) => (d, Some(t)),
59            None => (s, None),
60        };
61        // `[±]CCYY-MM-DD`: the year has ≥4 digits and an optional sign; month/day
62        // are exactly two digits (§9.1.1 — leading zeros may not be omitted).
63        let (sign, rest) = match date.strip_prefix('-') {
64            Some(r) => (-1, r),
65            None => (1, date.strip_prefix('+').unwrap_or(date)),
66        };
67        let mut dp = rest.split('-');
68        let y_str = dp.next().ok_or_else(invalid)?;
69        let m_str = dp.next().ok_or_else(invalid)?;
70        let d_str = dp.next().ok_or_else(invalid)?;
71        if dp.next().is_some() || y_str.len() < 4 || !all_digits(y_str) {
72            return Err(invalid());
73        }
74        let year = sign * y_str.parse::<i64>().map_err(|_| invalid())?;
75        // `to_jd` computes `365·(year + 4800)` in `i64`; bound the year so a hostile
76        // DATE can't overflow that. Any real observation/epoch date is far inside ±10⁶.
77        if year.unsigned_abs() > 1_000_000 {
78            return Err(invalid());
79        }
80        let month = parse_fixed(m_str, 2).ok_or_else(invalid)?;
81        let day = parse_fixed(d_str, 2).ok_or_else(invalid)?;
82        if !(1..=12).contains(&month) || day < 1 || day > days_in_month(year, month) {
83            return Err(invalid());
84        }
85
86        let (mut hour, mut minute, mut second) = (0u32, 0u32, 0.0f64);
87        if let Some(t) = time {
88            let mut tp = t.split(':');
89            hour = parse_fixed(tp.next().ok_or_else(invalid)?, 2).ok_or_else(invalid)?;
90            minute = parse_fixed(tp.next().ok_or_else(invalid)?, 2).ok_or_else(invalid)?;
91            if let Some(sec) = tp.next() {
92                second = parse_seconds(sec).ok_or_else(invalid)?;
93            }
94            // Second 60 is the leap second; this type is scale-agnostic, so the
95            // "only in UTC" restriction is left to the caller.
96            if tp.next().is_some() || hour >= 24 || minute >= 60 || !(0.0..61.0).contains(&second) {
97                return Err(invalid());
98            }
99        }
100        Ok(Datetime {
101            year,
102            month,
103            day,
104            hour,
105            minute,
106            second,
107        })
108    }
109
110    /// Julian Date of this datetime, interpreting the fields in their own time
111    /// scale (no scale conversion is applied here).
112    pub fn to_jd(&self) -> f64 {
113        let day_fraction =
114            (self.hour as f64 * 3600.0 + self.minute as f64 * 60.0 + self.second) / SEC_PER_DAY;
115        gregorian_to_jdn(self.year, self.month as i64, self.day as i64) as f64 - 0.5 + day_fraction
116    }
117
118    /// Modified Julian Date (`JD − 2400000.5`).
119    pub fn to_mjd(&self) -> f64 {
120        self.to_jd() - MJD0
121    }
122
123    /// Build a datetime from a JD (inverse of [`Datetime::to_jd`]). A single `f64`
124    /// JD at present epochs (~2.46e6) resolves to ~0.1 ms, so the recovered second
125    /// carries that much rounding — fine for display, not for sub-ms timing.
126    pub fn from_jd(jd: f64) -> Datetime {
127        // Split into integer day (JDN at noon) and the fraction past midnight.
128        let z = (jd + 0.5).floor();
129        let frac = (jd + 0.5) - z;
130        let date = jdn_to_gregorian(z as i64);
131        let mut secs = frac * SEC_PER_DAY;
132        let hour = (secs / 3600.0).floor();
133        secs -= hour * 3600.0;
134        let minute = (secs / 60.0).floor();
135        secs -= minute * 60.0;
136        Datetime {
137            year: date.year,
138            month: date.month,
139            day: date.day,
140            hour: hour as u32,
141            minute: minute as u32,
142            second: secs,
143        }
144    }
145}
146
147/// True if `s` is non-empty and all ASCII digits.
148fn all_digits(s: &str) -> bool {
149    !s.is_empty() && s.bytes().all(|b| b.is_ascii_digit())
150}
151
152/// Parse a fixed-width all-digit field (§9.1.1 forbids omitted leading zeros, so
153/// the length must be exact).
154fn parse_fixed(s: &str, width: usize) -> Option<u32> {
155    (s.len() == width && all_digits(s))
156        .then(|| s.parse().ok())
157        .flatten()
158}
159
160/// Parse a `ss[.s…]` seconds field: exactly two integer digits, optional fraction.
161fn parse_seconds(s: &str) -> Option<f64> {
162    let (int, frac) = s.split_once('.').map_or((s, None), |(i, f)| (i, Some(f)));
163    if int.len() != 2 || !all_digits(int) || frac.is_some_and(|f| !all_digits(f)) {
164        return None;
165    }
166    s.parse().ok()
167}
168
169/// A reference epoch: Julian (`J2000.0`) or Besselian (`B1950.0`).
170#[derive(Debug, Clone, Copy, PartialEq)]
171pub enum Epoch {
172    /// Julian epoch in years (e.g. `2000.0`).
173    Julian(f64),
174    /// Besselian epoch in years (e.g. `1950.0`).
175    Besselian(f64),
176}
177
178impl Epoch {
179    /// Parse `J<year>` or `B<year>` (e.g. `'J2000.0'`, `'B1950.0'`).
180    pub fn parse(s: &str) -> Result<Epoch> {
181        let s = s.trim();
182        let invalid = || FitsError::InvalidValue {
183            card: format!("epoch '{s}'"),
184        };
185        let (tag, rest) = s.split_at(
186            s.char_indices()
187                .next()
188                .map(|(_, c)| c.len_utf8())
189                .unwrap_or(0),
190        );
191        let year: f64 = rest.parse().map_err(|_| invalid())?;
192        match tag {
193            "J" | "j" => Ok(Epoch::Julian(year)),
194            "B" | "b" => Ok(Epoch::Besselian(year)),
195            _ => Err(invalid()),
196        }
197    }
198
199    /// Julian Date of the epoch. Julian: `2451545.0 + (y−2000)·365.25`; Besselian:
200    /// `2415020.31352 + (y−1900)·365.242198781` (the tropical year).
201    pub fn to_jd(self) -> f64 {
202        match self {
203            Epoch::Julian(y) => 2_451_545.0 + (y - 2000.0) * 365.25,
204            Epoch::Besselian(y) => 2_415_020.313_52 + (y - 1900.0) * 365.242_198_781,
205        }
206    }
207
208    /// Modified Julian Date of the epoch.
209    pub fn to_mjd(self) -> f64 {
210        self.to_jd() - MJD0
211    }
212}
213
214/// A FITS time scale (`TIMESYS` / `CTYPEi`).
215#[derive(Debug, Clone, Copy, PartialEq, Eq)]
216pub enum TimeScale {
217    /// Coordinated Universal Time.
218    Utc,
219    /// Universal Time (UT1) — treated as UTC here (ΔUT1 needs an external table).
220    Ut1,
221    /// International Atomic Time.
222    Tai,
223    /// Terrestrial Time.
224    Tt,
225    /// Geocentric Coordinate Time.
226    Tcg,
227    /// Barycentric Dynamical Time.
228    Tdb,
229    /// Barycentric Coordinate Time.
230    Tcb,
231    /// GPS time.
232    Gps,
233    /// An unspecified local clock (`LOCAL`); no conversion is possible.
234    Local,
235}
236
237impl TimeScale {
238    /// Parse a `TIMESYS`/`CTYPE` time-scale string (case-insensitive); accepts the
239    /// `TDT`/`ET` → `TT` and `IAT` → `TAI` aliases. Unknown values map to `LOCAL`.
240    pub fn parse(s: &str) -> TimeScale {
241        // A high-precision value appends a realization in parentheses — `TT(TAI)`,
242        // `UTC(NIST)` (§9.2.1); strip it before matching the scale name.
243        let base = s.trim().split('(').next().unwrap_or("").trim();
244        match base.to_ascii_uppercase().as_str() {
245            "UTC" | "GMT" => TimeScale::Utc, // §9.2.1: GMT is continuous with UTC
246            "UT1" | "UT" => TimeScale::Ut1,
247            "TAI" | "IAT" => TimeScale::Tai,
248            "TT" | "TDT" | "ET" => TimeScale::Tt,
249            "TCG" => TimeScale::Tcg,
250            "TDB" => TimeScale::Tdb,
251            "TCB" => TimeScale::Tcb,
252            "GPS" => TimeScale::Gps,
253            _ => TimeScale::Local,
254        }
255    }
256
257    /// Convert a Julian Date in this scale to a JD in `target`, treating `UT1` as
258    /// `UTC` (ΔUT1 = 0). Use [`TimeScale::convert_dut1`] to supply a real ΔUT1.
259    pub fn convert(self, jd: f64, target: TimeScale) -> f64 {
260        self.convert_dut1(jd, target, 0.0)
261    }
262
263    /// Convert a Julian Date with an explicit `dut1 = UT1 − UTC` (seconds, from
264    /// IERS) so conversions to/from `UT1` are exact. `Local` passes through.
265    pub fn convert_dut1(self, jd: f64, target: TimeScale, dut1: f64) -> f64 {
266        if self == target || self == TimeScale::Local || target == TimeScale::Local {
267            return jd;
268        }
269        from_tt(self.to_tt(jd, dut1), target, dut1)
270    }
271
272    /// This scale's JD expressed as TT (the common pivot). `dut1 = UT1 − UTC` (s).
273    fn to_tt(self, jd: f64, dut1: f64) -> f64 {
274        match self {
275            TimeScale::Tt => jd,
276            TimeScale::Tai => jd + TT_TAI / SEC_PER_DAY,
277            TimeScale::Gps => jd + (TT_TAI + TAI_GPS) / SEC_PER_DAY,
278            TimeScale::Utc => jd + (TT_TAI + leap_seconds(jd - MJD0)) / SEC_PER_DAY,
279            // UT1 → UTC (subtract ΔUT1) → TT.
280            TimeScale::Ut1 => {
281                let utc = jd - dut1 / SEC_PER_DAY;
282                utc + (TT_TAI + leap_seconds(utc - MJD0)) / SEC_PER_DAY
283            }
284            TimeScale::Tcg => jd - L_G * (jd - T1977_JD),
285            TimeScale::Tdb => jd - tdb_minus_tt(jd) / SEC_PER_DAY,
286            TimeScale::Tcb => {
287                let tdb = jd - L_B * (jd - T1977_JD) + TDB_0 / SEC_PER_DAY;
288                tdb - tdb_minus_tt(tdb) / SEC_PER_DAY
289            }
290            // convert_dut1 returns early for Local, so it never reaches the pivot.
291            TimeScale::Local => unreachable!("convert_dut1 short-circuits Local"),
292        }
293    }
294}
295
296/// TT (as a JD) expressed in `target` — the inverse of [`TimeScale::to_tt`].
297fn from_tt(tt: f64, target: TimeScale, dut1: f64) -> f64 {
298    match target {
299        TimeScale::Tt => tt,
300        TimeScale::Tai => tt - TT_TAI / SEC_PER_DAY,
301        TimeScale::Gps => tt - (TT_TAI + TAI_GPS) / SEC_PER_DAY,
302        TimeScale::Utc => {
303            let tai = tt - TT_TAI / SEC_PER_DAY;
304            // leap is a function of UTC; one lookup at the TAI date suffices away
305            // from the ≤1 s boundary ambiguity inherent to UTC.
306            tai - leap_seconds(tai - MJD0) / SEC_PER_DAY
307        }
308        // TT → UTC → UT1 (add ΔUT1).
309        TimeScale::Ut1 => {
310            let tai = tt - TT_TAI / SEC_PER_DAY;
311            let utc = tai - leap_seconds(tai - MJD0) / SEC_PER_DAY;
312            utc + dut1 / SEC_PER_DAY
313        }
314        TimeScale::Tcg => tt + L_G * (tt - T1977_JD),
315        TimeScale::Tdb => tt + tdb_minus_tt(tt) / SEC_PER_DAY,
316        TimeScale::Tcb => {
317            let tdb = tt + tdb_minus_tt(tt) / SEC_PER_DAY;
318            tdb + L_B * (tdb - T1977_JD) - TDB_0 / SEC_PER_DAY
319        }
320        // convert_dut1 returns early for Local, so it never reaches the pivot.
321        TimeScale::Local => unreachable!("convert_dut1 short-circuits Local"),
322    }
323}
324
325/// `TDB − TT` in seconds — the standard periodic approximation (~10 µs accuracy):
326/// `0.001658·sin g + 0.000014·sin 2g`, `g = 357.53° + 0.9856003°·(JD_TT − J2000)`.
327fn tdb_minus_tt(jd_tt: f64) -> f64 {
328    let g = (357.53 + 0.985_600_3 * (jd_tt - 2_451_545.0)).to_radians();
329    0.001_658 * g.sin() + 0.000_014 * (2.0 * g).sin()
330}
331
332/// `TAI − UTC` in seconds for a given UTC MJD: the integer leap-second count from
333/// the IERS table (1972–2017).
334///
335/// Outside that range the value is **frozen at the nearest table end**, which the
336/// caller should be aware of: past the last entry it returns the latest count
337/// (correct until a new leap second is announced), and before 1972 it returns the
338/// first entry (10 s) even though pre-1972 UTC used fractional "rubber seconds", so
339/// UTC↔TAI for pre-1972 dates is approximate, not exact.
340fn leap_seconds(mjd: f64) -> f64 {
341    // (year, month, day, TAI−UTC) at each step, 1972 onward.
342    const TABLE: &[(i64, i64, i64, f64)] = &[
343        (1972, 1, 1, 10.0),
344        (1972, 7, 1, 11.0),
345        (1973, 1, 1, 12.0),
346        (1974, 1, 1, 13.0),
347        (1975, 1, 1, 14.0),
348        (1976, 1, 1, 15.0),
349        (1977, 1, 1, 16.0),
350        (1978, 1, 1, 17.0),
351        (1979, 1, 1, 18.0),
352        (1980, 1, 1, 19.0),
353        (1981, 7, 1, 20.0),
354        (1982, 7, 1, 21.0),
355        (1983, 7, 1, 22.0),
356        (1985, 7, 1, 23.0),
357        (1988, 1, 1, 24.0),
358        (1990, 1, 1, 25.0),
359        (1991, 1, 1, 26.0),
360        (1992, 7, 1, 27.0),
361        (1993, 7, 1, 28.0),
362        (1994, 7, 1, 29.0),
363        (1996, 1, 1, 30.0),
364        (1997, 7, 1, 31.0),
365        (1999, 1, 1, 32.0),
366        (2006, 1, 1, 33.0),
367        (2009, 1, 1, 34.0),
368        (2012, 7, 1, 35.0),
369        (2015, 7, 1, 36.0),
370        (2017, 1, 1, 37.0),
371    ];
372    let mut leap = TABLE[0].3;
373    for &(y, m, d, l) in TABLE {
374        let threshold = gregorian_to_jdn(y, m, d) as f64 - 0.5 - MJD0;
375        if mjd >= threshold {
376            leap = l;
377        } else {
378            break;
379        }
380    }
381    leap
382}
383
384/// Julian Day Number at noon of a proleptic-Gregorian calendar date (the standard
385/// integer formula).
386/// Whether `year` is a leap year in the proleptic Gregorian calendar.
387fn is_leap_year(year: i64) -> bool {
388    year % 4 == 0 && (year % 100 != 0 || year % 400 == 0)
389}
390
391/// Number of days in `month` (1–12) of `year`; `0` for an out-of-range month.
392fn days_in_month(year: i64, month: u32) -> u32 {
393    match month {
394        1 | 3 | 5 | 7 | 8 | 10 | 12 => 31,
395        4 | 6 | 9 | 11 => 30,
396        2 if is_leap_year(year) => 29,
397        2 => 28,
398        _ => 0,
399    }
400}
401
402fn gregorian_to_jdn(year: i64, month: i64, day: i64) -> i64 {
403    let a = (14 - month) / 12;
404    let y = year + 4800 - a;
405    let m = month + 12 * a - 3;
406    day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045
407}
408
409/// A proleptic-Gregorian calendar date — the result of [`jdn_to_gregorian`].
410#[derive(Debug, Clone, Copy, PartialEq)]
411struct CalendarDate {
412    year: i64,
413    month: u32,
414    day: u32,
415}
416
417/// Inverse of [`gregorian_to_jdn`]: JDN → calendar date.
418fn jdn_to_gregorian(jdn: i64) -> CalendarDate {
419    let a = jdn + 32044;
420    let b = (4 * a + 3) / 146097;
421    let c = a - (146097 * b) / 4;
422    let d = (4 * c + 3) / 1461;
423    let e = c - (1461 * d) / 4;
424    let m = (5 * e + 2) / 153;
425    let day = e - (153 * m + 2) / 5 + 1;
426    let month = m + 3 - 12 * (m / 10);
427    let year = 100 * b + d - 4800 + m / 10;
428    CalendarDate {
429        year,
430        month: month as u32,
431        day: day as u32,
432    }
433}
434
435/// A time from a `JEPOCH`/`BEPOCH` keyword: its MJD and the scale the keyword
436/// implies (TDB for `JEPOCH`, ET ≈ TT for `BEPOCH`).
437#[derive(Debug, Clone, Copy, PartialEq)]
438pub struct EpochTime {
439    pub mjd: f64,
440    pub scale: TimeScale,
441}
442
443/// The global bound / duration / error time keywords (§9.4, §9.5, §9.7), as read
444/// by [`Header::time_bounds`](crate::Header::time_bounds). Start/end are absolute
445/// MJD; the rest are in `TIMEUNIT`.
446#[derive(Debug, Clone, Copy, PartialEq)]
447pub struct TimeBounds {
448    /// Observation start: `MJD-BEG`, else `DATE-BEG` → MJD.
449    pub beg_mjd: Option<f64>,
450    /// Observation end: `MJD-END`, else `DATE-END` → MJD.
451    pub end_mjd: Option<f64>,
452    /// Observation midpoint: `MJD-AVG`, else `DATE-AVG` → MJD (§9.5, Table 35).
453    pub avg_mjd: Option<f64>,
454    /// `XPOSURE` — effective exposure time.
455    pub xposure: Option<f64>,
456    /// `TELAPSE` — total elapsed time.
457    pub telapse: Option<f64>,
458    /// `TIMEDEL` — time resolution / bin width.
459    pub timedel: Option<f64>,
460    /// `TIMEPIXR` — pixel position within a bin (0–1, default 0.5).
461    pub timepixr: f64,
462    /// `TIMSYER` — systematic time error.
463    pub timsyer: Option<f64>,
464    /// `TIMRDER` — random time error.
465    pub timrder: Option<f64>,
466}
467
468/// A Good Time Interval as absolute MJD (§9.7).
469#[derive(Debug, Clone, Copy, PartialEq)]
470pub struct GtiInterval {
471    pub start_mjd: f64,
472    pub stop_mjd: f64,
473}
474
475/// The §9.6 `'PHASE'` axis folding parameters: the zero-phase reference time
476/// `CZPHSia` and the period `CPERIia`, in `TIMEUNIT` relative to the time
477/// reference (the `TSTART` convention).
478#[derive(Debug, Clone, Copy, PartialEq)]
479pub struct PhaseAxis {
480    pub zero_phase: f64,
481    pub period: f64,
482}
483
484impl PhaseAxis {
485    /// Fold a relative time (in `TIMEUNIT`) to a phase in `[0, 1)`; `0.0` if the
486    /// period is zero.
487    pub fn fold(&self, time: f64) -> f64 {
488        if self.period == 0.0 {
489            return 0.0;
490        }
491        ((time - self.zero_phase) / self.period).rem_euclid(1.0)
492    }
493}
494
495/// A header's time coordinate frame (§9): the reference epoch, scale, unit, and
496/// the resolved global time keywords.
497#[derive(Debug, Clone)]
498pub struct FitsTime {
499    /// `TIMESYS` time scale (default `UTC`).
500    pub scale: TimeScale,
501    /// Reference epoch as MJD (from `MJDREF`/`MJDREFI`+`MJDREFF`, `JDREF*`, or
502    /// `DATEREF`); `0.0` if none is given.
503    pub mjdref: f64,
504    /// `TIMEUNIT` (default `'s'`).
505    pub timeunit: String,
506    /// `TIMEOFFS` (§9.4.1): a uniform additive clock correction in `TIMEUNIT`,
507    /// equivalent to shifting the reference time. Default `0.0`.
508    pub timeoffs: f64,
509    /// `TREFPOS` (reference position, e.g. `'TOPOCENTER'`), if present.
510    pub trefpos: Option<String>,
511}
512
513impl FitsTime {
514    /// Parse the time frame from a header. The public entry point is
515    /// [`Header::time`](crate::Header::time), which forwards here.
516    pub(crate) fn from_header(header: &Header) -> FitsTime {
517        let scale = header
518            .get_text("TIMESYS")
519            .map(TimeScale::parse)
520            .unwrap_or(TimeScale::Utc);
521        let timeunit = header.get_text("TIMEUNIT").unwrap_or("s").to_string();
522        let trefpos = header.get_text("TREFPOS").map(str::to_string);
523        FitsTime {
524            scale,
525            mjdref: reference_mjd(header),
526            timeunit,
527            timeoffs: header.get_real("TIMEOFFS").unwrap_or(0.0),
528            trefpos,
529        }
530    }
531
532    /// `TIMEUNIT` expressed in seconds (`s`, `d`/day, `a`/`yr` Julian year).
533    pub fn unit_seconds(&self) -> f64 {
534        // Table 34. The deprecated tropical/Besselian years use their conventional
535        // lengths; a truly unknown unit falls back to seconds (the default).
536        match self.timeunit.trim() {
537            "s" => 1.0,
538            "min" => 60.0,
539            "h" => 3600.0,
540            "d" | "day" => SEC_PER_DAY,
541            "a" | "yr" | "y" => 365.25 * SEC_PER_DAY, // Julian year
542            "cy" => 36525.0 * SEC_PER_DAY,            // Julian century = 100 a
543            "ta" => 365.24219 * SEC_PER_DAY,          // tropical year (deprecated)
544            "Ba" => 365.2421988 * SEC_PER_DAY,        // Besselian year (deprecated)
545            _ => 1.0,
546        }
547    }
548
549    /// Resolve a time value measured *relative* to `MJDREF` (e.g. `TSTART`,
550    /// `TSTOP`), in `TIMEUNIT`, to an absolute MJD in the frame's own scale. The
551    /// `TIMEOFFS` clock correction (§9.4.1) is added before scaling.
552    pub fn relative_to_mjd(&self, value: f64) -> f64 {
553        self.mjdref + (value + self.timeoffs) * self.unit_seconds() / SEC_PER_DAY
554    }
555
556    /// The observation MJD from `MJD-OBS`, else `DATE-OBS`, else `None`. Reads only
557    /// the header (not the parsed frame). The public entry point is
558    /// [`Header::obs_mjd`](crate::Header::obs_mjd), which forwards here.
559    pub(crate) fn obs_mjd(header: &Header) -> Option<f64> {
560        if let Some(mjd) = header.get_real("MJD-OBS") {
561            return Some(mjd);
562        }
563        if let Some(d) = header
564            .get_text("DATE-OBS")
565            .and_then(|s| Datetime::parse(s).ok())
566        {
567            return Some(d.to_mjd());
568        }
569        // §9.5: with no DATE-OBS/MJD-OBS, JEPOCH/BEPOCH stand in for the observation time.
570        Self::epoch(header).map(|e| e.mjd)
571    }
572
573    /// The Julian (`JEPOCH`, implied scale TDB) or Besselian (`BEPOCH`, implied
574    /// scale ET ≈ TT) epoch keyword as an [`EpochTime`], if present (§9.1.2, §9.5).
575    /// `JEPOCH` wins if both appear. Reads only the header, so it takes no `self`.
576    pub(crate) fn epoch(header: &Header) -> Option<EpochTime> {
577        if let Some(j) = header.get_real("JEPOCH") {
578            return Some(EpochTime {
579                mjd: Epoch::Julian(j).to_mjd(),
580                scale: TimeScale::Tdb,
581            });
582        }
583        let b = header.get_real("BEPOCH")?;
584        Some(EpochTime {
585            mjd: Epoch::Besselian(b).to_mjd(),
586            scale: TimeScale::Tt, // ET ≈ TT
587        })
588    }
589
590    /// The global bound / duration / error time keywords (§9.4, §9.5, §9.7). The
591    /// start/end are resolved to absolute MJD (`MJD-BEG`/`-END`, else `DATE-BEG`/
592    /// `-END`); durations and errors are returned as stored, in `TIMEUNIT`. Reads
593    /// only the header, so it takes no `self`.
594    pub(crate) fn bounds(header: &Header) -> TimeBounds {
595        let mjd_or_date = |mjd: &str, date: &str| {
596            header.get_real(mjd).or_else(|| {
597                header
598                    .get_text(date)
599                    .and_then(|s| Datetime::parse(s).ok())
600                    .map(|d| d.to_mjd())
601            })
602        };
603        TimeBounds {
604            beg_mjd: mjd_or_date("MJD-BEG", "DATE-BEG"),
605            end_mjd: mjd_or_date("MJD-END", "DATE-END"),
606            avg_mjd: mjd_or_date("MJD-AVG", "DATE-AVG"),
607            xposure: header.get_real("XPOSURE"),
608            telapse: header.get_real("TELAPSE"),
609            timedel: header.get_real("TIMEDEL"),
610            timepixr: header.get_real("TIMEPIXR").unwrap_or(0.5), // §9.4.2 default
611            timsyer: header.get_real("TIMSYER"),
612            timrder: header.get_real("TIMRDER"),
613        }
614    }
615
616    /// Convert Good Time Interval `START`/`STOP` column values (relative to
617    /// `MJDREF`, in `TIMEUNIT`) to absolute-MJD intervals (§9.7). Pairs are taken
618    /// element-wise up to the shorter slice.
619    pub fn gti_intervals(&self, starts: &[f64], stops: &[f64]) -> Vec<GtiInterval> {
620        starts
621            .iter()
622            .zip(stops)
623            .map(|(&s, &e)| GtiInterval {
624                start_mjd: self.relative_to_mjd(s),
625                stop_mjd: self.relative_to_mjd(e),
626            })
627            .collect()
628    }
629
630    /// If WCS axis `axis` (1-based) is a time axis (`CTYPEi = 'TIME'` or a
631    /// time-scale name, §9.2.3), convert a 1-based pixel coordinate along it to an
632    /// absolute MJD in the frame's scale: the linear axis value (elapsed time in
633    /// `TIMEUNIT` from `MJDREF`) plus the reference. `None` if not a time axis.
634    pub fn time_axis_mjd(&self, header: &Header, axis: usize, pixel: f64) -> Option<f64> {
635        let ctype = header.get_text(key!("CTYPE{axis}").as_str())?;
636        if !is_time_ctype(ctype) {
637            return None;
638        }
639        let crpix = header.get_real(key!("CRPIX{axis}").as_str()).unwrap_or(0.0);
640        let crval = header.get_real(key!("CRVAL{axis}").as_str()).unwrap_or(0.0);
641        let cdelt = header.get_real(key!("CDELT{axis}").as_str()).unwrap_or(1.0);
642        Some(self.relative_to_mjd(crval + cdelt * (pixel - crpix)))
643    }
644
645    /// The §9.6 `'PHASE'` axis parameters (`CZPHSia` zero-phase time, `CPERIia`
646    /// period) for WCS axis `axis` (1-based), or `None` if it is not a phase axis.
647    /// Reads only the header, so it takes no `self`.
648    pub(crate) fn phase_axis(header: &Header, axis: usize) -> Option<PhaseAxis> {
649        let ctype = header.get_text(key!("CTYPE{axis}").as_str())?;
650        if TimeAxisKind::from_ctype(ctype) != Some(TimeAxisKind::Phase) {
651            return None;
652        }
653        Some(PhaseAxis {
654            zero_phase: header.get_real(key!("CZPHS{axis}").as_str()).unwrap_or(0.0),
655            period: header.get_real(key!("CPERI{axis}").as_str()).unwrap_or(0.0),
656        })
657    }
658}
659
660/// A time-related WCS axis type (§9.6).
661#[derive(Debug, Clone, Copy, PartialEq, Eq)]
662pub enum TimeAxisKind {
663    /// `'TIME'` or a time-scale name — an absolute time axis (→ MJD).
664    Time,
665    /// `'PHASE'` — phase folded on a period (`CPERIia`, zero `CZPHSia`).
666    Phase,
667    /// `'TIMELAG'` — a correlation/cross-spectral time lag.
668    Timelag,
669    /// `'FREQUENCY'` — a frequency axis.
670    Frequency,
671}
672
673impl TimeAxisKind {
674    /// Classify a `CTYPE` as a time-related axis (§9.6), or `None` if it is not one.
675    pub fn from_ctype(ctype: &str) -> Option<TimeAxisKind> {
676        let head = ctype.split('-').next().unwrap_or("").trim();
677        match head {
678            "TIME" => Some(TimeAxisKind::Time),
679            "PHASE" => Some(TimeAxisKind::Phase),
680            "TIMELAG" => Some(TimeAxisKind::Timelag),
681            "FREQUENCY" => Some(TimeAxisKind::Frequency),
682            _ if !matches!(TimeScale::parse(head), TimeScale::Local) => Some(TimeAxisKind::Time),
683            _ => None,
684        }
685    }
686}
687
688/// True if a `CTYPE` denotes an absolute *time* axis (`'TIME'` or a time-scale
689/// name) — the kind [`FitsTime::time_axis_mjd`] resolves to an MJD.
690fn is_time_ctype(ctype: &str) -> bool {
691    TimeAxisKind::from_ctype(ctype) == Some(TimeAxisKind::Time)
692}
693
694/// The reference epoch as MJD: `MJDREF` (or `MJDREFI`+`MJDREFF`), else `JDREF`
695/// (or `JDREFI`+`JDREFF`), else `DATEREF`, else `0.0`.
696fn reference_mjd(header: &Header) -> f64 {
697    if let Some(mjd) = resolve_split_ref(header, "MJDREF", "MJDREFI", "MJDREFF") {
698        return mjd;
699    }
700    if let Some(jd) = resolve_split_ref(header, "JDREF", "JDREFI", "JDREFF") {
701        return jd - MJD0;
702    }
703    header
704        .get_text("DATEREF")
705        .and_then(|s| Datetime::parse(s).ok())
706        .map(|d| d.to_mjd())
707        .unwrap_or(0.0)
708}
709
710/// Resolve a reference epoch from its single (`MJDREF`) and split-precision
711/// (`MJDREFI`+`MJDREFF`) keywords. Per §9.2.2 a *full* integer+fractional split
712/// takes precedence over the single value; otherwise the single value is used,
713/// falling back to a lone split part.
714fn resolve_split_ref(header: &Header, single: &str, int: &str, frac: &str) -> Option<f64> {
715    let i = header.get_real(int);
716    let f = header.get_real(frac);
717    match (i, f) {
718        (Some(i), Some(f)) => Some(i + f),
719        _ => header.get_real(single).or_else(|| match (i, f) {
720            (None, None) => None,
721            _ => Some(i.unwrap_or(0.0) + f.unwrap_or(0.0)),
722        }),
723    }
724}
725
726#[cfg(test)]
727mod tests;
fits_well/time/mod.rs

fits_well/time/
mod.rs
