1use super::ToTT;
63use crate::{
64 constants::FAIRHD,
65 julian::JulianDate,
66 scales::{TDB, TT},
67};
68
69use crate::{TimeError, TimeResult};
70use celestial_core::constants::{
71 DAYS_PER_JULIAN_MILLENNIUM, DEG_TO_RAD, J2000_JD, SECONDS_PER_DAY_F64, TWOPI,
72};
73use celestial_core::math::fmod;
74use celestial_core::Location;
75
76fn greenwich_location() -> Location {
81 Location::from_degrees(51.477928, 0.0, 46.0).expect("Greenwich coordinates should be valid")
82}
83
84const DEFAULT_UT1_OFFSET_SECONDS: f64 = 0.0;
85
86pub fn compute_tdb_tt_offset(
101 date_jd: &JulianDate,
102 ut1_fraction: f64,
103 location: &Location,
104) -> TimeResult<f64> {
105 let (u, v) = location.to_geocentric_km()?;
106
107 let dtr = calculate_tdb_tt_difference(
108 date_jd.jd1(),
109 date_jd.jd2(),
110 ut1_fraction,
111 location.longitude,
112 u,
113 v,
114 );
115
116 Ok(dtr)
117}
118
119pub trait ToTDB {
131 fn to_tdb_greenwich(&self) -> TimeResult<TDB>;
133 fn to_tdb_with_location(&self, location: &Location) -> TimeResult<TDB>;
135 fn to_tdb_with_location_and_ut1_offset(
137 &self,
138 location: &Location,
139 ut1_offset_seconds: f64,
140 ) -> TimeResult<TDB>;
141 fn to_tdb_with_offset(&self, dtr_seconds: f64) -> TimeResult<TDB>;
143}
144
145impl ToTDB for TDB {
146 fn to_tdb_greenwich(&self) -> TimeResult<TDB> {
147 Ok(*self)
148 }
149
150 fn to_tdb_with_location(&self, _location: &Location) -> TimeResult<TDB> {
151 Ok(*self)
152 }
153
154 fn to_tdb_with_location_and_ut1_offset(
155 &self,
156 _location: &Location,
157 _ut1_offset_seconds: f64,
158 ) -> TimeResult<TDB> {
159 Ok(*self)
160 }
161
162 fn to_tdb_with_offset(&self, _dtr_seconds: f64) -> TimeResult<TDB> {
163 Ok(*self)
164 }
165}
166
167impl ToTDB for TT {
168 fn to_tdb_greenwich(&self) -> TimeResult<TDB> {
169 let location = greenwich_location();
170 self.to_tdb_with_location_and_ut1_offset(&location, DEFAULT_UT1_OFFSET_SECONDS)
171 }
172
173 fn to_tdb_with_location(&self, location: &Location) -> TimeResult<TDB> {
174 self.to_tdb_with_location_and_ut1_offset(location, DEFAULT_UT1_OFFSET_SECONDS)
175 }
176
177 fn to_tdb_with_location_and_ut1_offset(
178 &self,
179 location: &Location,
180 ut1_offset_seconds: f64,
181 ) -> TimeResult<TDB> {
182 let tt_jd = self.to_julian_date();
183
184 let tt_f64 = tt_jd.to_f64();
185 let ut1_fraction = ((tt_f64 - libm::trunc(tt_f64)) * SECONDS_PER_DAY_F64
186 + ut1_offset_seconds)
187 / SECONDS_PER_DAY_F64;
188 let ut1_fraction = ut1_fraction - libm::floor(ut1_fraction);
189
190 let dtr = compute_tdb_tt_offset(&tt_jd, ut1_fraction, location)?;
191 self.to_tdb_with_offset(dtr)
192 }
193
194 fn to_tdb_with_offset(&self, dtr_seconds: f64) -> TimeResult<TDB> {
195 let tt_jd = self.to_julian_date();
196
197 let dtr_days = dtr_seconds / celestial_core::constants::SECONDS_PER_DAY_F64;
198
199 let (tdb_jd1, tdb_jd2) = if tt_jd.jd1().abs() > tt_jd.jd2().abs() {
200 (tt_jd.jd1(), tt_jd.jd2() + dtr_days)
201 } else {
202 (tt_jd.jd1() + dtr_days, tt_jd.jd2())
203 };
204
205 let tdb_jd = JulianDate::new(tdb_jd1, tdb_jd2);
206 Ok(TDB::from_julian_date(tdb_jd))
207 }
208}
209
210impl ToTT for TDB {
211 fn to_tt(&self) -> TimeResult<TT> {
212 Err(TimeError::ConversionError(
213 "TDB→TT conversion requires observer location. \
214 Use to_tt_greenwich() for Greenwich or to_tt_with_location() for other locations."
215 .to_string(),
216 ))
217 }
218}
219
220pub trait ToTTFromTDB {
238 fn to_tt_greenwich(&self) -> TimeResult<TT>;
240 fn to_tt_with_location(&self, location: &Location) -> TimeResult<TT>;
242 fn to_tt_with_location_and_ut1_offset(
244 &self,
245 location: &Location,
246 ut1_offset_seconds: f64,
247 ) -> TimeResult<TT>;
248 fn to_tt_with_offset(&self, dtr_seconds: f64) -> TimeResult<TT>;
250}
251
252impl ToTTFromTDB for TDB {
253 fn to_tt_greenwich(&self) -> TimeResult<TT> {
254 let location = greenwich_location();
255 self.to_tt_with_location_and_ut1_offset(&location, DEFAULT_UT1_OFFSET_SECONDS)
256 }
257
258 fn to_tt_with_location(&self, location: &Location) -> TimeResult<TT> {
259 self.to_tt_with_location_and_ut1_offset(location, DEFAULT_UT1_OFFSET_SECONDS)
260 }
261
262 fn to_tt_with_location_and_ut1_offset(
263 &self,
264 location: &Location,
265 ut1_offset_seconds: f64,
266 ) -> TimeResult<TT> {
267 let tdb_jd = self.to_julian_date();
268
269 let tdb_f64 = tdb_jd.to_f64();
270 let ut1_fraction_approx = ((tdb_f64 - libm::trunc(tdb_f64)) * SECONDS_PER_DAY_F64
271 + ut1_offset_seconds)
272 / SECONDS_PER_DAY_F64;
273 let ut1_fraction_approx = ut1_fraction_approx - libm::floor(ut1_fraction_approx);
274
275 let dtr_approx = compute_tdb_tt_offset(&tdb_jd, ut1_fraction_approx, location)?;
276
277 let tt_approx = self.to_tt_with_offset(dtr_approx)?;
278 let tt_jd_approx = tt_approx.to_julian_date();
279
280 let tt_jd_approx_f64 = tt_jd_approx.jd1() + tt_jd_approx.jd2();
281 let ut1_fraction = ((tt_jd_approx_f64 - libm::trunc(tt_jd_approx_f64))
282 * SECONDS_PER_DAY_F64
283 + ut1_offset_seconds)
284 / SECONDS_PER_DAY_F64;
285 let ut1_fraction = ut1_fraction - libm::floor(ut1_fraction);
286
287 let dtr_refined = compute_tdb_tt_offset(&tt_jd_approx, ut1_fraction, location)?;
288
289 let tt_refined = self.to_tt_with_offset(dtr_refined)?;
290 let tt_jd_refined = tt_refined.to_julian_date();
291
292 let tt_jd_refined_f64 = tt_jd_refined.jd1() + tt_jd_refined.jd2();
293 let ut1_fraction_final = ((tt_jd_refined_f64 - libm::trunc(tt_jd_refined_f64))
294 * SECONDS_PER_DAY_F64
295 + ut1_offset_seconds)
296 / SECONDS_PER_DAY_F64;
297 let ut1_fraction_final = ut1_fraction_final - libm::floor(ut1_fraction_final);
298
299 let dtr_final = compute_tdb_tt_offset(&tt_jd_refined, ut1_fraction_final, location)?;
300
301 self.to_tt_with_offset(dtr_final)
302 }
303
304 fn to_tt_with_offset(&self, dtr_seconds: f64) -> TimeResult<TT> {
305 let tdb_jd = self.to_julian_date();
306
307 let dtr_days = dtr_seconds / celestial_core::constants::SECONDS_PER_DAY_F64;
308
309 let (tt_jd1, tt_jd2) = if tdb_jd.jd1().abs() > tdb_jd.jd2().abs() {
310 (tdb_jd.jd1(), tdb_jd.jd2() - dtr_days)
311 } else {
312 (tdb_jd.jd1() - dtr_days, tdb_jd.jd2())
313 };
314
315 let tt_jd = JulianDate::new(tt_jd1, tt_jd2);
316 Ok(TT::from_julian_date(tt_jd))
317 }
318}
319
320fn calculate_tdb_tt_difference(date1: f64, date2: f64, ut: f64, elong: f64, u: f64, v: f64) -> f64 {
338 let t = ((date1 - J2000_JD) + date2) / DAYS_PER_JULIAN_MILLENNIUM;
339
340 let tsol = fmod(ut, 1.0) * TWOPI + elong;
341
342 let w = t / 3600.0;
343
344 let elsun = fmod(280.46645683 + 1296027711.03429 * w, 360.0) * DEG_TO_RAD;
345
346 let emsun = fmod(357.52910918 + 1295965810.481 * w, 360.0) * DEG_TO_RAD;
347
348 let d = fmod(297.85019547 + 16029616012.090 * w, 360.0) * DEG_TO_RAD;
349
350 let elj = fmod(34.35151874 + 109306899.89453 * w, 360.0) * DEG_TO_RAD;
351
352 let els = fmod(50.07744430 + 44046398.47038 * w, 360.0) * DEG_TO_RAD;
353
354 let wt = 0.00029e-10 * u * libm::sin(tsol + elsun - els)
355 + 0.00100e-10 * u * libm::sin(tsol - 2.0 * emsun)
356 + 0.00133e-10 * u * libm::sin(tsol - d)
357 + 0.00133e-10 * u * libm::sin(tsol + elsun - elj)
358 - 0.00229e-10 * u * libm::sin(tsol + 2.0 * elsun + emsun)
359 - 0.02200e-10 * v * libm::cos(elsun + emsun)
360 + 0.05312e-10 * u * libm::sin(tsol - emsun)
361 - 0.13677e-10 * u * libm::sin(tsol + 2.0 * elsun)
362 - 1.31840e-10 * v * libm::cos(elsun)
363 + 3.17679e-10 * u * libm::sin(tsol);
364
365 let mut w0 = 0.0;
366 for j in (0..474).rev() {
367 w0 += FAIRHD[j][0] * libm::sin(FAIRHD[j][1] * t + FAIRHD[j][2]);
368 }
369
370 let mut w1 = 0.0;
371 for j in (474..679).rev() {
372 w1 += FAIRHD[j][0] * libm::sin(FAIRHD[j][1] * t + FAIRHD[j][2]);
373 }
374
375 let mut w2 = 0.0;
376 for j in (679..764).rev() {
377 if FAIRHD[j][0] != 0.0 {
378 w2 += FAIRHD[j][0] * libm::sin(FAIRHD[j][1] * t + FAIRHD[j][2]);
379 }
380 }
381
382 let mut w3 = 0.0;
383 for j in (764..784).rev() {
384 if FAIRHD[j][0] != 0.0 {
385 w3 += FAIRHD[j][0] * libm::sin(FAIRHD[j][1] * t + FAIRHD[j][2]);
386 }
387 }
388
389 let mut w4 = 0.0;
390 for j in (784..787).rev() {
391 if FAIRHD[j][0] != 0.0 {
392 w4 += FAIRHD[j][0] * libm::sin(FAIRHD[j][1] * t + FAIRHD[j][2]);
393 }
394 }
395
396 let wf = t * (t * (t * (t * w4 + w3) + w2) + w1) + w0;
397
398 let wj = 0.00065e-6 * libm::sin(6069.776754 * t + 4.021194)
399 + 0.00033e-6 * libm::sin(213.299095 * t + 5.543132)
400 - 0.00196e-6 * libm::sin(6208.294251 * t + 5.696701)
401 - 0.00173e-6 * libm::sin(74.781599 * t + 2.435900)
402 + 0.03638e-6 * t * t;
403
404 wt + wf + wj
405}
406
407#[cfg(test)]
408mod tests {
409 use super::*;
410 use celestial_core::constants::{DAYS_PER_JULIAN_CENTURY, J2000_JD};
411
412 #[test]
413 fn test_identity_conversions() {
414 let tdb = TDB::from_julian_date(JulianDate::new(J2000_JD, 0.999999999999999));
415 let location = Location::from_degrees(45.0, 90.0, 100.0).unwrap();
416
417 let via_greenwich = tdb.to_tdb_greenwich().unwrap();
418 assert_eq!(
419 tdb.to_julian_date().jd1(),
420 via_greenwich.to_julian_date().jd1(),
421 "TDB→TDB via Greenwich should preserve JD1"
422 );
423 assert_eq!(
424 tdb.to_julian_date().jd2(),
425 via_greenwich.to_julian_date().jd2(),
426 "TDB→TDB via Greenwich should preserve JD2"
427 );
428
429 let via_location = tdb.to_tdb_with_location(&location).unwrap();
430 assert_eq!(
431 tdb.to_julian_date().jd1(),
432 via_location.to_julian_date().jd1(),
433 "TDB→TDB via location should preserve JD1"
434 );
435 assert_eq!(
436 tdb.to_julian_date().jd2(),
437 via_location.to_julian_date().jd2(),
438 "TDB→TDB via location should preserve JD2"
439 );
440
441 let via_ut1_offset = tdb
442 .to_tdb_with_location_and_ut1_offset(&location, 0.3)
443 .unwrap();
444 assert_eq!(
445 tdb.to_julian_date().jd1(),
446 via_ut1_offset.to_julian_date().jd1(),
447 "TDB→TDB via UT1 offset should preserve JD1"
448 );
449 assert_eq!(
450 tdb.to_julian_date().jd2(),
451 via_ut1_offset.to_julian_date().jd2(),
452 "TDB→TDB via UT1 offset should preserve JD2"
453 );
454
455 let via_offset = tdb.to_tdb_with_offset(0.001).unwrap();
456 assert_eq!(
457 tdb.to_julian_date().jd1(),
458 via_offset.to_julian_date().jd1(),
459 "TDB→TDB via offset should preserve JD1"
460 );
461 assert_eq!(
462 tdb.to_julian_date().jd2(),
463 via_offset.to_julian_date().jd2(),
464 "TDB→TDB via offset should preserve JD2"
465 );
466 }
467
468 #[test]
469 fn test_tt_tdb_offset_verification() {
470 let test_cases = [
471 (J2000_JD, "J2000.0", greenwich_location()),
472 (
473 J2000_JD - DAYS_PER_JULIAN_CENTURY,
474 "1900",
475 greenwich_location(),
476 ),
477 (J2000_JD + 18262.5, "2050", greenwich_location()),
478 (
479 J2000_JD + DAYS_PER_JULIAN_CENTURY,
480 "2100",
481 greenwich_location(),
482 ),
483 (
484 J2000_JD,
485 "J2000 Tokyo",
486 Location::from_degrees(35.6762, 139.6503, 40.0).unwrap(),
487 ),
488 (
489 J2000_JD,
490 "J2000 Sydney",
491 Location::from_degrees(-33.8688, 151.2093, 58.0).unwrap(),
492 ),
493 ];
494
495 for (jd, description, location) in test_cases {
496 let tt = TT::from_julian_date(JulianDate::new(jd, 0.0));
497 let tdb = tt.to_tdb_with_location(&location).unwrap();
498
499 let diff_seconds = (tdb.to_julian_date().to_f64() - tt.to_julian_date().to_f64())
500 * SECONDS_PER_DAY_F64;
501
502 assert!(
503 diff_seconds.abs() < 0.002,
504 "{}: TT→TDB offset should be < 2ms, got {:.6} seconds",
505 description,
506 diff_seconds
507 );
508
509 let tdb = TDB::from_julian_date(JulianDate::new(jd, 0.0));
510 let tt = tdb.to_tt_with_location(&location).unwrap();
511
512 let diff_seconds = (tdb.to_julian_date().to_f64() - tt.to_julian_date().to_f64())
513 * SECONDS_PER_DAY_F64;
514
515 assert!(
516 diff_seconds.abs() < 0.002,
517 "{}: TDB→TT offset should be < 2ms, got {:.6} seconds",
518 description,
519 diff_seconds
520 );
521 }
522 }
523
524 #[test]
525 fn test_tt_tdb_round_trip_precision() {
526 const TOLERANCE_DAYS: f64 = 1e-11;
527
528 let test_jd2_values = [0.0, 0.5, 0.123456789012345, -0.123456789012345];
529
530 for jd2 in test_jd2_values {
531 let original_tt = TT::from_julian_date(JulianDate::new(J2000_JD, jd2));
532 let tdb = original_tt.to_tdb_greenwich().unwrap();
533 let round_trip_tt = tdb.to_tt_greenwich().unwrap();
534
535 let jd1_diff =
536 (original_tt.to_julian_date().jd1() - round_trip_tt.to_julian_date().jd1()).abs();
537 let jd2_diff =
538 (original_tt.to_julian_date().jd2() - round_trip_tt.to_julian_date().jd2()).abs();
539
540 assert!(
541 jd1_diff < TOLERANCE_DAYS,
542 "TT→TDB→TT JD1 diff {} exceeds tolerance for jd2={}",
543 jd1_diff,
544 jd2
545 );
546 assert!(
547 jd2_diff < TOLERANCE_DAYS,
548 "TT→TDB→TT JD2 diff {} exceeds tolerance for jd2={}",
549 jd2_diff,
550 jd2
551 );
552
553 let original_tdb = TDB::from_julian_date(JulianDate::new(J2000_JD, jd2));
554 let tt = original_tdb.to_tt_greenwich().unwrap();
555 let round_trip_tdb = tt.to_tdb_greenwich().unwrap();
556
557 let jd1_diff =
558 (original_tdb.to_julian_date().jd1() - round_trip_tdb.to_julian_date().jd1()).abs();
559 let jd2_diff =
560 (original_tdb.to_julian_date().jd2() - round_trip_tdb.to_julian_date().jd2()).abs();
561
562 assert!(
563 jd1_diff < TOLERANCE_DAYS,
564 "TDB→TT→TDB JD1 diff {} exceeds tolerance for jd2={}",
565 jd1_diff,
566 jd2
567 );
568 assert!(
569 jd2_diff < TOLERANCE_DAYS,
570 "TDB→TT→TDB JD2 diff {} exceeds tolerance for jd2={}",
571 jd2_diff,
572 jd2
573 );
574 }
575
576 let alt_tt = TT::from_julian_date(JulianDate::new(0.5, J2000_JD));
577 let alt_tdb = alt_tt.to_tdb_greenwich().unwrap();
578 let alt_round_trip = alt_tdb.to_tt_greenwich().unwrap();
579
580 let jd1_diff =
581 (alt_tt.to_julian_date().jd1() - alt_round_trip.to_julian_date().jd1()).abs();
582 let jd2_diff =
583 (alt_tt.to_julian_date().jd2() - alt_round_trip.to_julian_date().jd2()).abs();
584
585 assert!(
586 jd1_diff < TOLERANCE_DAYS,
587 "Alternate split TT→TDB→TT JD1 diff {} exceeds tolerance",
588 jd1_diff
589 );
590 assert!(
591 jd2_diff < TOLERANCE_DAYS,
592 "Alternate split TT→TDB→TT JD2 diff {} exceeds tolerance",
593 jd2_diff
594 );
595 }
596
597 #[test]
598 fn test_api_equivalence() {
599 let tt = TT::from_julian_date(JulianDate::new(J2000_JD, 0.123456));
600
601 let tdb_greenwich = tt.to_tdb_greenwich().unwrap();
602 let tdb_explicit = tt.to_tdb_with_location(&greenwich_location()).unwrap();
603
604 assert_eq!(
605 tdb_greenwich.to_julian_date().jd1(),
606 tdb_explicit.to_julian_date().jd1(),
607 "TT: to_tdb_greenwich() should match to_tdb_with_location(greenwich) JD1"
608 );
609 assert_eq!(
610 tdb_greenwich.to_julian_date().jd2(),
611 tdb_explicit.to_julian_date().jd2(),
612 "TT: to_tdb_greenwich() should match to_tdb_with_location(greenwich) JD2"
613 );
614
615 let tdb = TDB::from_julian_date(JulianDate::new(J2000_JD, 0.987654));
616
617 let tt_greenwich = tdb.to_tt_greenwich().unwrap();
618 let tt_explicit = tdb.to_tt_with_location(&greenwich_location()).unwrap();
619
620 assert_eq!(
621 tt_greenwich.to_julian_date().jd1(),
622 tt_explicit.to_julian_date().jd1(),
623 "TDB: to_tt_greenwich() should match to_tt_with_location(greenwich) JD1"
624 );
625 assert_eq!(
626 tt_greenwich.to_julian_date().jd2(),
627 tt_explicit.to_julian_date().jd2(),
628 "TDB: to_tt_greenwich() should match to_tt_with_location(greenwich) JD2"
629 );
630 }
631
632 #[test]
633 fn test_tdb_to_tt_requires_location() {
634 let tdb = TDB::from_julian_date(JulianDate::new(J2000_JD, 0.0));
635 let result = tdb.to_tt();
636
637 assert!(result.is_err(), "TDB.to_tt() should return error");
638
639 match result {
640 Err(TimeError::ConversionError(msg)) => {
641 assert!(
642 msg.contains("location"),
643 "Error message should mention location requirement: {}",
644 msg
645 );
646 assert!(
647 msg.contains("greenwich") || msg.contains("Greenwich"),
648 "Error message should mention Greenwich option: {}",
649 msg
650 );
651 }
652 _ => panic!("Expected ConversionError, got {:?}", result),
653 }
654 }
655}