Skip to main content

deep_time/dt/
tdb.rs

1use crate::{Dt, Real, Scale, TDB0_ATTOS, TT_TAI_OFFSET, sin};
2
3impl Dt {
4    /// DE440/LTE440-tuned compact analytical TT–TDB model
5    ///
6    /// Exact 13-term Fourier decomposition from LTE440 (Lu et al. 2025, Table 3)
7    /// + physical VSOP2013 annual term + tiny JPL secular corrections.
8    pub const fn tdb_minus_tt(seconds_since_j2000_tt: Real) -> Real {
9        // J2000.0 = 2000-01-01 12:00:00 TT → 100 Julian years = exactly 3_155_760_000 s
10        const J2000_SEC_PER_MILLENNIUM: Real = 31_557_600_000.0;
11
12        let t = seconds_since_j2000_tt / J2000_SEC_PER_MILLENNIUM; // centuries since J2000
13        let mut correction = f!(0.0);
14
15        // Physical annual term (VSOP2013 secular e(t) — replaces LTE440 term #1)
16        let g = f!(6283.075849991) * t + f!(6.240054195);
17        let e = f!(0.0167086342) - f!(0.0004203654) * t - f!(0.0000126734) * t * t
18            + f!(0.0000001444) * t * t * t
19            - f!(0.0000000002) * t * t * t * t
20            + f!(0.0000000003) * t * t * t * t * t;
21        let k = f!(0.09897232);
22        let varpi = f!(-0.00000257) - f!(0.05551247) * t;
23        correction += k * e * sin(g + varpi + f!(0.01671) * sin(g));
24
25        // Exact LTE440 Fourier terms #2–#13 (all amplitudes >1 µs from DE440 item #15)
26        const LTE440_TERMS: [(Real, Real, Real); 12] = [
27            (0.00012630813184, 77713.771468120, 5.18472464), // #2 D (lunar synodic)
28            (0.00001937467715, 5753.384884897, 1.33855843),  // #3 E–J (Earth–Jupiter)
29            (0.00001370088760, 12566.151699983, 3.07602294), // #4 2E (semi-annual)
30            (0.00000747520418, 5574.656149776, 3.32446352),  // #5 D–L
31            (0.00000424397312, 4320.34946237, 3.43186281),   // #6 J (Jupiter)
32            (0.00000376051430, 377.97977422, 0.92358639),    // #7 E–S (Earth–Saturn)
33            (0.00000293368121, 161002.466707021, 1.09317212), // #8 D+L
34            (0.00000267752983, 6208.659051973, 1.51225314),  // #9 E–U (Earth–Uranus)
35            (0.00000236687890, 71430.993657045, 5.21748801), // #10 E–D (Earth–Moon difference)
36            (0.00000185820098, 211.334300759, 2.56843762),   // #11 S (Saturn long-period)
37            (0.00000109742615, 3929.675646567, 4.67635157),  // #12 V–E (Venus–Earth)
38            (0.00000108850698, 7859.351293133, 2.99248981),  // #13 2V–2E
39        ];
40
41        let mut i = 0;
42        while i < 12 {
43            let (amp, freq, phase) = LTE440_TERMS[i];
44            correction += amp * sin(freq * t + phase);
45            i += 1;
46        }
47
48        // Tiny JPL wj mass adjustments + quadratic secular (<1 ns)
49        correction += f!(0.00000000065) * sin(f!(6069.776754) * t + f!(4.021194));
50        correction += f!(0.00000000033) * sin(f!(213.299095) * t + f!(5.543132));
51        correction += f!(-0.00000000196) * sin(f!(6208.294251) * t + f!(5.696701));
52        correction += f!(-0.00000000173) * sin(f!(74.781599) * t + f!(2.435900));
53        correction += f!(0.00000003638) * t * t; // quadratic secular
54
55        correction
56    }
57
58    /// Converts a TAI [`Dt`] to TDB.
59    pub const fn tai_to_tdb(tai: Self) -> Self {
60        let tt = tai.add(TT_TAI_OFFSET);
61        let correction = Self::tdb_minus_tt(tt.to_sec_f());
62        tt.add(Dt::from_sec_f(correction))
63    }
64
65    /// Converts a TDB [`Dt`] to TAI.
66    pub const fn tdb_to_tai(tdb: Self) -> Self {
67        // Linear-rate + constant initial guess (dominant part of the forward transformation)
68        let elapsed = Self::to_attos_since_tcg_tcb_epoch(tdb);
69        let linear_span = Self::mul_lb(elapsed); // LB * elapsed
70        let mut tt = tdb
71            .sub(Dt::from_attos(linear_span, Scale::TAI))
72            .sub(Dt::from_attos(TDB0_ATTOS, Scale::TAI));
73
74        // Fixed-point iteration: TT_{n+1} = TDB − P(TT_n)
75        let mut i = 0u32;
76        while i < 8 {
77            let p = Self::tdb_minus_tt(tt.to_sec_f());
78            let new_tt = tdb.sub(Dt::from_sec_f(p));
79
80            // Early exit when change is smaller than ~1 atto-second
81            let delta = new_tt.to_diff_raw(tt);
82            if delta.sec == 0 && delta.attos < 1 {
83                tt = new_tt;
84                break;
85            }
86
87            tt = new_tt;
88            i += 1;
89        }
90
91        tt.sub(TT_TAI_OFFSET)
92    }
93}