tempoch-core 0.6.0

Core astronomical time primitives for tempoch.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

// SPDX-License-Identifier: AGPL-3.0-only
// Copyright (C) 2026 Vallés Puig, Ramon

use super::utc_tai::time_data_tai_minus_utc_mjd_extrapolated;
use crate::earth::delta_t::delta_t_seconds_from_modern_points;
use crate::earth::eop::EopValues;
use crate::foundation::error::ConversionError;
use qtty::Day as DayQuantity;
use qtty::Second;
use tempoch_time_data::{EopPoint, TimeDataBundle};

pub(crate) fn time_data_delta_t(
    data: &TimeDataBundle,
    jd_ut: DayQuantity,
) -> Result<Second, ConversionError> {
    delta_t_seconds_from_modern_points(jd_ut, data.modern_delta_t_points())
}

pub(crate) fn time_data_eop_at(data: &TimeDataBundle, mjd_utc: DayQuantity) -> Option<EopValues> {
    let points = data.eop_points();
    let first = points.first()?.mjd;
    let last = points.last()?.mjd;
    let mjd_f = mjd_utc.value();
    let lo_i = mjd_f.floor() as i32;
    let hi_i = lo_i + 1;
    if lo_i < first || lo_i > last {
        return None;
    }
    let lo = find_eop_point(points, lo_i)?;
    let hi = if hi_i > last {
        lo
    } else {
        find_eop_point(points, hi_i)?
    };

    let frac = if lo.mjd == hi.mjd {
        0.0
    } else {
        mjd_f - lo_i as f64
    };
    let lerp = |a: f64, b: f64| a + frac * (b - a);
    let lerp_opt = |a: Option<f64>, b: Option<f64>| match (a, b) {
        (Some(a), Some(b)) => Some(lerp(a, b)),
        _ => None,
    };
    let lod_milliseconds = match (lo.lod_milliseconds, hi.lod_milliseconds) {
        (Some(a), Some(b)) => Some(lerp(a, b)),
        _ => None,
    };

    let ut1_minus_utc = {
        // Allow extrapolation here: these calls are for internal ΔT bookkeeping
        // (correcting EOP-derived UT1-UTC to the actual UTC-TAI offset), not for
        // validating UTC representations. Pre-1961 EOP data is rare but valid.
        let lo_offset =
            time_data_tai_minus_utc_mjd_extrapolated(data, DayQuantity::new(lo_i as f64));
        let hi_offset =
            time_data_tai_minus_utc_mjd_extrapolated(data, DayQuantity::new(hi_i as f64));
        let query_offset = time_data_tai_minus_utc_mjd_extrapolated(data, mjd_utc);
        match (lo_offset, hi_offset, query_offset) {
            (Some(lo_tmu), Some(hi_tmu), Some(query_tmu)) => {
                let lo_cont = lo.ut1_minus_utc_seconds - lo_tmu.value();
                let hi_cont = hi.ut1_minus_utc_seconds - hi_tmu.value();
                Second::new(lerp(lo_cont, hi_cont) + query_tmu.value())
            }
            _ => Second::new(lerp(lo.ut1_minus_utc_seconds, hi.ut1_minus_utc_seconds)),
        }
    };

    Some(EopValues {
        mjd_utc,
        pm_xp_arcsec: lerp_opt(lo.pm_xp_arcsec, hi.pm_xp_arcsec),
        pm_yp_arcsec: lerp_opt(lo.pm_yp_arcsec, hi.pm_yp_arcsec),
        ut1_minus_utc,
        lod_milliseconds,
        dx_milliarcsec: lerp_opt(lo.dx_milliarcsec, hi.dx_milliarcsec),
        dy_milliarcsec: lerp_opt(lo.dy_milliarcsec, hi.dy_milliarcsec),
        ut1_observed: lo.ut1_observed && hi.ut1_observed,
    })
}

fn find_eop_point(points: &[EopPoint], mjd: i32) -> Option<EopPoint> {
    let idx = points.partition_point(|point| point.mjd < mjd);
    let point = *points.get(idx)?;
    (point.mjd == mjd).then_some(point)
}