celestial-time 0.1.1-alpha.2

Pure Rust astronomical time scales
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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
//! Conversions between Coordinated Universal Time (UTC) and Universal Time (UT1).
//!
//! UT1 is the principal form of Universal Time, directly tied to Earth's rotation angle.
//! UTC is the civil time standard maintained by atomic clocks. The difference between
//! them, called DUT1 (= UT1 - UTC), is published by the IERS in Bulletin A.
//!
//! # The DUT1 Offset
//!
//! DUT1 measures how much Earth's actual rotation deviates from the uniform UTC clock:
//!
//! ```text
//! UT1 = UTC + DUT1
//! ```
//!
//! The IERS keeps |DUT1| < 0.9 seconds by inserting leap seconds into UTC. DUT1 values
//! are published weekly in IERS Bulletin A with ~1 ms precision. For high-precision
//! applications (astrometry, VLBI, satellite tracking), the correct DUT1 must be
//! obtained from IERS data for the specific date.
//!
//! # Conversion Paths
//!
//! This module provides two paths from UT1 to UTC:
//!
//! ```text
//! Direct:      UT1 ←→ UTC    (via DUT1 offset)
//! Via TAI:     UT1 → TAI → UTC    (for verification)
//! ```
//!
//! The direct path (`ToUTCWithDUT1`) handles leap second boundaries correctly by
//! adjusting the effective DUT1 value near discontinuities. The TAI path
//! (`ToUTCViaTAI`) chains through intermediate time scales and serves as a
//! verification mechanism.
//!
//! # Leap Second Handling
//!
//! The UT1→UTC conversion is complicated by leap seconds: when UTC inserts a leap
//! second, there's a discontinuity in the UTC-TAI offset. The `adjust_dut1_for_leap_second`
//! function scans nearby days for offset changes and smoothly interpolates the
//! correction across the leap second boundary.
//!
//! # Precision
//!
//! Round-trip conversions (UTC → UT1 → UTC or UT1 → UTC → UT1) achieve ~1 picosecond
//! accuracy when using consistent DUT1 values. The two-part Julian Date representation
//! preserves full f64 precision throughout.
//!
//! # Usage
//!
//! ```
//! use celestial_time::scales::{UTC, UT1};
//! use celestial_time::scales::conversions::{ToUT1WithDUT1, ToUTCWithDUT1};
//! use celestial_time::julian::JulianDate;
//! use celestial_core::constants::J2000_JD;
//!
//! // DUT1 for 2000-01-01 was approximately +0.3 seconds (from IERS Bulletin A)
//! let dut1 = 0.3;
//!
//! let utc = UTC::from_julian_date(JulianDate::new(J2000_JD, 0.0));
//! let ut1 = utc.to_ut1_with_dut1(dut1).unwrap();
//!
//! // UT1 should be DUT1 seconds ahead of UTC
//! let diff_days = ut1.to_julian_date().to_f64() - utc.to_julian_date().to_f64();
//! let diff_seconds = diff_days * 86400.0;
//! assert!((diff_seconds - dut1).abs() < 0.01);
//!
//! // Round-trip preserves the original time
//! let utc_back = ut1.to_utc_with_dut1(dut1).unwrap();
//! let round_trip_diff = (utc.to_julian_date().to_f64()
//!     - utc_back.to_julian_date().to_f64()).abs();
//! assert!(round_trip_diff < 1e-14);  // ~1 picosecond
//! ```
//!
//! # References
//!
//! - IERS Bulletin A: Weekly publication of UT1-UTC values
//! - IERS Conventions (2010): Chapter 5, Earth Rotation
//! - USNO Earth Orientation Parameters

use super::super::common::get_tai_utc_offset; // Direct import from common
use super::ut1_tai::{ToTAIWithOffset, ToUT1WithOffset};
use super::utc_tai::{calendar_to_julian, julian_to_calendar};
use super::{ToTAI, ToUTC};
use crate::julian::JulianDate;
use crate::scales::{UT1, UTC};
use crate::TimeResult;
use celestial_core::constants::SECONDS_PER_DAY_F64;

/// Convert to UT1 using a known DUT1 (UT1-UTC) offset.
///
/// DUT1 values must be obtained from IERS Bulletin A for the specific date.
/// The offset is typically in the range -0.9 to +0.9 seconds.
pub trait ToUT1WithDUT1 {
    /// Convert to UT1 given the DUT1 offset in seconds.
    ///
    /// # Arguments
    ///
    /// * `dut1_seconds` - The UT1-UTC offset in seconds (from IERS Bulletin A)
    ///
    /// # Returns
    ///
    /// The corresponding UT1 instant. The conversion chains through TAI:
    /// UTC → TAI → UT1, computing the UT1-TAI offset from DUT1 and the
    /// TAI-UTC offset for the date.
    fn to_ut1_with_dut1(&self, dut1_seconds: f64) -> TimeResult<UT1>;
}

/// Convert to UTC using a known DUT1 (UT1-UTC) offset.
///
/// This is the inverse of `ToUT1WithDUT1`. Given a UT1 instant and the
/// DUT1 offset, computes the corresponding UTC instant.
pub trait ToUTCWithDUT1 {
    /// Convert to UTC given the DUT1 offset in seconds.
    ///
    /// # Arguments
    ///
    /// * `dut1_seconds` - The UT1-UTC offset in seconds (from IERS Bulletin A)
    ///
    /// # Returns
    ///
    /// The corresponding UTC instant. Handles leap second boundaries by
    /// adjusting the effective DUT1 value near discontinuities.
    fn to_utc_with_dut1(&self, dut1_seconds: f64) -> TimeResult<UTC>;
}

impl ToUT1WithDUT1 for UTC {
    /// Convert UTC to UT1 by computing UT1-TAI from DUT1 and TAI-UTC.
    ///
    /// The conversion uses the relationship:
    ///
    /// ```text
    /// UT1 - TAI = DUT1 - (TAI - UTC) = DUT1 - TAI_UTC_offset
    /// ```
    ///
    /// The TAI-UTC offset is looked up from the leap second table for the
    /// specific date. The result chains: UTC → TAI → UT1.
    fn to_ut1_with_dut1(&self, dut1_seconds: f64) -> TimeResult<UT1> {
        let tai = self.to_tai()?;

        let utc_jd = self.to_julian_date();
        let (year, month, day, day_fraction) = julian_to_calendar(utc_jd.jd1(), utc_jd.jd2())?;
        let tai_utc_seconds = get_tai_utc_offset(year, month, day, day_fraction);
        let ut1_tai_offset = dut1_seconds - tai_utc_seconds;

        tai.to_ut1_with_offset(ut1_tai_offset)
    }
}

/// Adjust DUT1 for leap second discontinuities near the given Julian Date.
///
/// When converting UT1 to UTC near a leap second boundary, the naive subtraction
/// of DUT1 can place the result on the wrong side of the discontinuity. This
/// function detects nearby leap seconds and adjusts the effective DUT1 value
/// to smoothly interpolate across the boundary.
///
/// # Algorithm
///
/// 1. Scan days from (JD - 1) to (JD + 3) looking for TAI-UTC offset changes
/// 2. If a change > 0.5 seconds is found, a leap second occurred
/// 3. If the leap second and DUT1 have the same sign, subtract the leap from DUT1
/// 4. Compute the fraction of the way past the leap second boundary
/// 5. Gradually add back the leap second contribution based on that fraction
///
/// The range [-1, +3] days ensures leap seconds are detected whether the input
/// time is just before, during, or just after the discontinuity.
///
/// # Arguments
///
/// * `jd_big` - Larger magnitude component of the Julian Date
/// * `jd_small` - Smaller magnitude component of the Julian Date
/// * `dut1` - The raw DUT1 offset in seconds
///
/// # Returns
///
/// The adjusted DUT1 value that accounts for any nearby leap second.
fn adjust_dut1_for_leap_second(jd_big: f64, jd_small: f64, dut1: f64) -> TimeResult<f64> {
    let mut duts = dut1;
    let mut prev_offset = 0.0;

    for i in -1..=3 {
        let jd_frac = jd_small + i as f64;
        let (year, month, day, _) = julian_to_calendar(jd_big, jd_frac)?;
        let curr_offset = get_tai_utc_offset(year, month, day, 0.0);

        if i == -1 {
            prev_offset = curr_offset;
            continue;
        }

        let delta = curr_offset - prev_offset;
        if delta.abs() < 0.5 {
            prev_offset = curr_offset;
            continue;
        }

        // Found leap second boundary
        if delta * duts >= 0.0 {
            duts -= delta;
        }

        let (leap_d1, leap_d2) = calendar_to_julian(year, month, day);
        let time_past_leap =
            (jd_big - leap_d1) + (jd_small - (leap_d2 - 1.0 + duts / SECONDS_PER_DAY_F64));

        if time_past_leap > 0.0 {
            let fraction =
                (time_past_leap * SECONDS_PER_DAY_F64 / (SECONDS_PER_DAY_F64 + delta)).min(1.0);
            duts += delta * fraction;
        }
        break;
    }

    Ok(duts)
}

impl ToUTCWithDUT1 for UT1 {
    /// Convert UT1 to UTC by subtracting the adjusted DUT1 offset.
    ///
    /// The conversion:
    ///
    /// 1. Determines which JD component has larger magnitude (for precision)
    /// 2. Adjusts DUT1 for any nearby leap second boundaries
    /// 3. Subtracts the adjusted DUT1 from the smaller-magnitude component
    /// 4. Preserves the original JD component ordering
    ///
    /// The leap second adjustment ensures correct behavior at discontinuities
    /// where the UTC scale gains an extra second.
    fn to_utc_with_dut1(&self, dut1_seconds: f64) -> TimeResult<UTC> {
        let ut1_jd = self.to_julian_date();
        let (big, small, big_first) = if ut1_jd.jd1().abs() >= ut1_jd.jd2().abs() {
            (ut1_jd.jd1(), ut1_jd.jd2(), true)
        } else {
            (ut1_jd.jd2(), ut1_jd.jd1(), false)
        };

        let adjusted_dut1 = adjust_dut1_for_leap_second(big, small, dut1_seconds)?;
        let small_corrected = small - adjusted_dut1 / SECONDS_PER_DAY_F64;

        let (utc_jd1, utc_jd2) = if big_first {
            (big, small_corrected)
        } else {
            (small_corrected, big)
        };
        Ok(UTC::from_julian_date(JulianDate::new(utc_jd1, utc_jd2)))
    }
}

/// Alternative UT1 to UTC conversion path that chains through TAI.
///
/// This trait provides a verification mechanism: the direct path (`ToUTCWithDUT1`)
/// and the TAI path should produce identical results. Any discrepancy indicates
/// a bug in one of the conversion implementations.
pub trait ToUTCViaTAI {
    /// Convert UT1 to UTC by chaining: UT1 → TAI → UTC.
    ///
    /// # Arguments
    ///
    /// * `dut1_seconds` - The UT1-UTC offset in seconds (from IERS Bulletin A)
    ///
    /// # Algorithm
    ///
    /// 1. Compute UT1-TAI offset from DUT1 and TAI-UTC for the date
    /// 2. Convert UT1 to TAI using the computed offset
    /// 3. Convert TAI to UTC using the leap second table
    ///
    /// This path uses the same TAI-UTC lookup as the direct conversion but
    /// exercises different code paths, making it useful for testing.
    fn to_utc_via_tai_with_dut1(&self, dut1_seconds: f64) -> TimeResult<UTC>;
}

impl ToUTCViaTAI for UT1 {
    fn to_utc_via_tai_with_dut1(&self, dut1_seconds: f64) -> TimeResult<UTC> {
        let ut1_jd = self.to_julian_date();
        let (year, month, day, day_fraction) = julian_to_calendar(ut1_jd.jd1(), ut1_jd.jd2())?;
        let tai_utc_seconds = get_tai_utc_offset(year, month, day, day_fraction);
        let ut1_tai_offset = dut1_seconds - tai_utc_seconds;

        let tai = self.to_tai_with_offset(ut1_tai_offset)?;

        tai.to_utc()
    }
}

#[cfg(test)]
mod tests {
    use super::super::{ToUT1, ToUTC};
    use super::*;
    use celestial_core::constants::J2000_JD;

    #[test]
    fn test_identity_conversions() {
        let ut1 = UT1::from_julian_date(JulianDate::new(J2000_JD, 0.5));
        let identity_ut1 = ut1.to_ut1().unwrap();
        assert_eq!(
            ut1.to_julian_date().jd1(),
            identity_ut1.to_julian_date().jd1()
        );
        assert_eq!(
            ut1.to_julian_date().jd2(),
            identity_ut1.to_julian_date().jd2()
        );

        let utc = UTC::from_julian_date(JulianDate::new(J2000_JD, 0.5));
        let identity_utc = utc.to_utc().unwrap();
        assert_eq!(
            utc.to_julian_date().jd1(),
            identity_utc.to_julian_date().jd1()
        );
        assert_eq!(
            utc.to_julian_date().jd2(),
            identity_utc.to_julian_date().jd2()
        );
    }

    #[test]
    fn test_dut1_offset_relationship() {
        let dut1_values = [-0.9, 0.0, 0.9];

        for dut1 in dut1_values {
            let utc = UTC::from_julian_date(JulianDate::new(J2000_JD, 0.0));
            let ut1 = utc.to_ut1_with_dut1(dut1).unwrap();

            let ut1_jd = ut1.to_julian_date().to_f64();
            let diff_seconds = (ut1_jd - J2000_JD) * SECONDS_PER_DAY_F64;

            assert!(
                diff_seconds > -1.0 && diff_seconds < 1.0,
                "DUT1={}: UT1-UTC difference should be within 1 second: {} seconds",
                dut1,
                diff_seconds
            );

            let ut1_reverse = UT1::from_julian_date(JulianDate::new(J2000_JD, 0.0));
            let utc_reverse = ut1_reverse.to_utc_with_dut1(dut1).unwrap();
            let utc_jd = utc_reverse.to_julian_date().to_f64();
            let reverse_diff = (J2000_JD - utc_jd) * SECONDS_PER_DAY_F64;

            assert!(
                (reverse_diff - dut1).abs() < 0.1,
                "DUT1={}: UTC should be behind UT1 by ~DUT1: {} seconds",
                dut1,
                reverse_diff
            );
        }

        let ut1_normal = UT1::from_julian_date(JulianDate::new(J2000_JD, 0.5));
        let utc_normal = ut1_normal.to_utc_with_dut1(0.3).unwrap();
        assert!(
            utc_normal.to_julian_date().jd1().abs() > utc_normal.to_julian_date().jd2().abs(),
            "Should preserve larger JD1 component"
        );

        let ut1_flipped = UT1::from_julian_date(JulianDate::new(0.1, J2000_JD));
        let utc_flipped = ut1_flipped.to_utc_with_dut1(0.3).unwrap();
        assert!(
            utc_flipped.to_julian_date().jd2().abs() > utc_flipped.to_julian_date().jd1().abs(),
            "Should preserve larger JD2 component"
        );
    }

    #[test]
    fn test_utc_ut1_round_trip_precision() {
        let tolerance = 1e-14; // ~1 picosecond

        let jd_splits = [
            (J2000_JD, 0.123456789),
            (J2000_JD, 0.5),
            (0.1, J2000_JD),
            (J2000_JD, 0.999999999),
            (J2000_JD, 0.25),
        ];

        let dut1_values = [-0.9, 0.0, 0.9];

        for (jd1, jd2) in jd_splits {
            for dut1 in dut1_values {
                let original_utc = UTC::from_julian_date(JulianDate::new(jd1, jd2));
                let ut1 = original_utc.to_ut1_with_dut1(dut1).unwrap();
                let round_trip_utc = ut1.to_utc_with_dut1(dut1).unwrap();

                let diff = (original_utc.to_julian_date().to_f64()
                    - round_trip_utc.to_julian_date().to_f64())
                .abs();
                assert!(
                    diff < tolerance,
                    "UTC->UT1->UTC round trip (jd1={}, jd2={}, dut1={}): {:.2e} days exceeds {:.0e}",
                    jd1,
                    jd2,
                    dut1,
                    diff,
                    tolerance
                );

                let original_ut1 = UT1::from_julian_date(JulianDate::new(jd1, jd2));
                let utc = original_ut1.to_utc_with_dut1(dut1).unwrap();
                let round_trip_ut1 = utc.to_ut1_with_dut1(dut1).unwrap();

                let diff_reverse = (original_ut1.to_julian_date().to_f64()
                    - round_trip_ut1.to_julian_date().to_f64())
                .abs();
                assert!(
                    diff_reverse < tolerance,
                    "UT1->UTC->UT1 round trip (jd1={}, jd2={}, dut1={}): {:.2e} days exceeds {:.0e}",
                    jd1,
                    jd2,
                    dut1,
                    diff_reverse,
                    tolerance
                );
            }
        }
    }

    #[test]
    fn test_leap_second_boundary_handling() {
        let leap_dates = [
            (2441499.5, "1972-07-01"),
            (2441683.5, "1973-01-01"),
            (2442048.5, "1974-01-01"),
        ];

        for (jd, description) in leap_dates {
            let ut1_at_leap = UT1::from_julian_date(JulianDate::new(jd, 0.0));
            let utc = ut1_at_leap.to_utc_with_dut1(0.0).unwrap();
            assert!(
                utc.to_julian_date().to_f64() > 0.0,
                "{}: Should produce valid UTC at leap second",
                description
            );

            let ut1_after_leap = UT1::from_julian_date(JulianDate::new(jd, 0.001));
            let utc_after = ut1_after_leap.to_utc_with_dut1(0.0).unwrap();
            assert!(
                utc_after.to_julian_date().to_f64() > jd,
                "{}: UTC should be after leap second start",
                description
            );

            let utc_direct = ut1_at_leap.to_utc_with_dut1(0.0).unwrap();
            let utc_via_tai = ut1_at_leap.to_utc_via_tai_with_dut1(0.0).unwrap();
            let diff = (utc_direct.to_julian_date().to_f64()
                - utc_via_tai.to_julian_date().to_f64())
            .abs();
            assert!(
                diff < 1e-14,
                "{}: Direct vs TAI-intermediate should match: {:.2e} days",
                description,
                diff
            );
        }
    }
}