1#![allow(clippy::unreadable_literal)]
7#![allow(clippy::many_single_char_names)]
8
9use crate::math::{floor, polynomial};
10use crate::{Error, Result};
11#[cfg(feature = "chrono")]
12use chrono::{Datelike, TimeZone, Timelike};
13
14const SECONDS_PER_DAY: f64 = 86_400.0;
16
17const J2000_JDN: f64 = 2_451_545.0;
19
20const DAYS_PER_CENTURY: f64 = 36_525.0;
22
23pub(crate) fn validate_utc_components(
25 year: i32,
26 month: u32,
27 day: u32,
28 hour: u32,
29 minute: u32,
30 second: f64,
31) -> Result<()> {
32 if !(1..=12).contains(&month) {
33 return Err(Error::InvalidDateTime {
34 message: "month must be between 1 and 12",
35 });
36 }
37 if !(1..=31).contains(&day) {
38 return Err(Error::InvalidDateTime {
39 message: "day must be between 1 and 31",
40 });
41 }
42 if hour > 23 {
43 return Err(Error::InvalidDateTime {
44 message: "hour must be between 0 and 23",
45 });
46 }
47 if minute > 59 {
48 return Err(Error::InvalidDateTime {
49 message: "minute must be between 0 and 59",
50 });
51 }
52 if !(0.0..60.0).contains(&second) {
53 return Err(Error::InvalidDateTime {
54 message: "second must be between 0 and 59.999...",
55 });
56 }
57 if year == 1582 && month == 10 && (5..=14).contains(&day) {
58 return Err(Error::InvalidDateTime {
59 message: "dates 1582-10-05 through 1582-10-14 do not exist in Gregorian calendar",
60 });
61 }
62 if day > days_in_month(year, month) {
63 return Err(Error::InvalidDateTime {
64 message: "day is out of range for month",
65 });
66 }
67
68 Ok(())
69}
70
71#[derive(Debug, Clone, Copy, PartialEq)]
76pub struct JulianDate {
77 jd: f64,
79 delta_t: f64,
81}
82
83impl JulianDate {
84 #[cfg(feature = "chrono")]
98 #[cfg_attr(docsrs, doc(cfg(feature = "chrono")))]
99 pub fn from_datetime<Tz: TimeZone>(
100 datetime: &chrono::DateTime<Tz>,
101 delta_t: f64,
102 ) -> Result<Self> {
103 let utc_datetime = datetime.with_timezone(&chrono::Utc);
105 Self::from_utc(
106 utc_datetime.year(),
107 utc_datetime.month(),
108 utc_datetime.day(),
109 utc_datetime.hour(),
110 utc_datetime.minute(),
111 f64::from(utc_datetime.second()) + f64::from(utc_datetime.nanosecond()) / 1e9,
112 delta_t,
113 )
114 }
115
116 pub fn from_utc(
141 year: i32,
142 month: u32,
143 day: u32,
144 hour: u32,
145 minute: u32,
146 second: f64,
147 delta_t: f64,
148 ) -> Result<Self> {
149 validate_utc_components(year, month, day, hour, minute, second)?;
150 if !delta_t.is_finite() {
151 return Err(Error::InvalidDateTime {
152 message: "delta_t must be finite",
153 });
154 }
155
156 let jd = calculate_julian_date(year, month, day, hour, minute, second);
157 Ok(Self { jd, delta_t })
158 }
159
160 pub fn from_utc_simple(
176 year: i32,
177 month: u32,
178 day: u32,
179 hour: u32,
180 minute: u32,
181 second: f64,
182 ) -> Result<Self> {
183 Self::from_utc(year, month, day, hour, minute, second, 0.0)
184 }
185
186 #[must_use]
191 pub const fn julian_date(&self) -> f64 {
192 self.jd
193 }
194
195 #[must_use]
200 pub const fn delta_t(&self) -> f64 {
201 self.delta_t
202 }
203
204 #[must_use]
211 pub fn julian_ephemeris_day(&self) -> f64 {
212 self.jd + self.delta_t / SECONDS_PER_DAY
213 }
214
215 #[must_use]
222 pub fn julian_century(&self) -> f64 {
223 (self.jd - J2000_JDN) / DAYS_PER_CENTURY
224 }
225
226 #[must_use]
233 pub fn julian_ephemeris_century(&self) -> f64 {
234 (self.julian_ephemeris_day() - J2000_JDN) / DAYS_PER_CENTURY
235 }
236
237 #[must_use]
244 pub fn julian_ephemeris_millennium(&self) -> f64 {
245 self.julian_ephemeris_century() / 10.0
246 }
247
248 pub(crate) fn add_days(self, days: f64) -> Self {
250 Self {
251 jd: self.jd + days,
252 delta_t: self.delta_t,
253 }
254 }
255}
256
257fn calculate_julian_date(
262 year: i32,
263 month: u32,
264 day: u32,
265 hour: u32,
266 minute: u32,
267 second: f64,
268) -> f64 {
269 let mut y = year;
270 let mut m = i32::try_from(month).expect("month should be valid i32");
271
272 if m < 3 {
274 y -= 1;
275 m += 12;
276 }
277
278 let d = f64::from(day) + (f64::from(hour) + (f64::from(minute) + second / 60.0) / 60.0) / 24.0;
280
281 let mut jd =
283 floor(365.25 * (f64::from(y) + 4716.0)) + floor(30.6001 * f64::from(m + 1)) + d - 1524.5;
284
285 if jd >= 2_299_161.0 {
288 let a = floor(f64::from(y) / 100.0);
289 let b = 2.0 - a + floor(a / 4.0);
290 jd += b;
291 }
292
293 jd
294}
295
296fn days_in_month(year: i32, month: u32) -> u32 {
297 let is_leap_year = if year > 1582 || (year == 1582 && month > 10) {
298 (year % 4 == 0 && year % 100 != 0) || year % 400 == 0
299 } else {
300 year % 4 == 0
301 };
302
303 match month {
304 1 | 3 | 5 | 7 | 8 | 10 | 12 => 31,
305 4 | 6 | 9 | 11 => 30,
306 2 if is_leap_year => 29,
307 2 => 28,
308 _ => unreachable!("month already validated"),
309 }
310}
311
312pub struct DeltaT;
317
318impl DeltaT {
319 #[allow(clippy::too_many_lines)] pub fn estimate(decimal_year: f64) -> Result<f64> {
341 let year = decimal_year;
342
343 if !year.is_finite() {
344 return Err(Error::InvalidDateTime {
345 message: "year must be finite",
346 });
347 }
348
349 if year < -500.0 {
350 return Err(Error::InvalidDateTime {
351 message: "ΔT estimates not available before year -500",
352 });
353 }
354
355 let delta_t = if year < 500.0 {
356 let u = year / 100.0;
357 polynomial(
358 &[
359 10583.6,
360 -1014.41,
361 33.78311,
362 -5.952053,
363 -0.1798452,
364 0.022174192,
365 0.0090316521,
366 ],
367 u,
368 )
369 } else if year < 1600.0 {
370 let u = (year - 1000.0) / 100.0;
371 polynomial(
372 &[
373 1574.2,
374 -556.01,
375 71.23472,
376 0.319781,
377 -0.8503463,
378 -0.005050998,
379 0.0083572073,
380 ],
381 u,
382 )
383 } else if year < 1700.0 {
384 let t = year - 1600.0;
385 polynomial(&[120.0, -0.9808, -0.01532, 1.0 / 7129.0], t)
386 } else if year < 1800.0 {
387 let t = year - 1700.0;
388 polynomial(
389 &[8.83, 0.1603, -0.0059285, 0.00013336, -1.0 / 1_174_000.0],
390 t,
391 )
392 } else if year < 1860.0 {
393 let t = year - 1800.0;
394 polynomial(
395 &[
396 13.72,
397 -0.332447,
398 0.0068612,
399 0.0041116,
400 -0.00037436,
401 0.0000121272,
402 -0.0000001699,
403 0.000000000875,
404 ],
405 t,
406 )
407 } else if year < 1900.0 {
408 let t = year - 1860.0;
409 polynomial(
410 &[
411 7.62,
412 0.5737,
413 -0.251754,
414 0.01680668,
415 -0.0004473624,
416 1.0 / 233_174.0,
417 ],
418 t,
419 )
420 } else if year < 1920.0 {
421 let t = year - 1900.0;
422 polynomial(&[-2.79, 1.494119, -0.0598939, 0.0061966, -0.000197], t)
423 } else if year < 1941.0 {
424 let t = year - 1920.0;
425 polynomial(&[21.20, 0.84493, -0.076100, 0.0020936], t)
426 } else if year < 1961.0 {
427 let t = year - 1950.0;
428 polynomial(&[29.07, 0.407, -1.0 / 233.0, 1.0 / 2547.0], t)
429 } else if year < 1986.0 {
430 let t = year - 1975.0;
431 polynomial(&[45.45, 1.067, -1.0 / 260.0, -1.0 / 718.0], t)
432 } else if year < 2005.0 {
433 let t = year - 2000.0;
434 polynomial(
435 &[
436 63.86,
437 0.3345,
438 -0.060374,
439 0.0017275,
440 0.000651814,
441 0.00002373599,
442 ],
443 t,
444 )
445 } else if year < 2015.0 {
446 let t = year - 2005.0;
447 polynomial(&[64.69, 0.2930], t)
448 } else if year <= 3000.0 {
449 let t = year - 2015.0;
450 polynomial(&[67.62, 0.3645, 0.0039755], t)
451 } else {
452 return Err(Error::InvalidDateTime {
453 message: "ΔT estimates not available beyond year 3000",
454 });
455 };
456
457 Ok(delta_t)
458 }
459
460 pub fn estimate_from_date(year: i32, month: u32) -> Result<f64> {
477 if !(1..=12).contains(&month) {
478 return Err(Error::InvalidDateTime {
479 message: "month must be between 1 and 12",
480 });
481 }
482
483 let decimal_year = f64::from(year) + (f64::from(month) - 0.5) / 12.0;
484 Self::estimate(decimal_year)
485 }
486
487 #[cfg(feature = "chrono")]
516 #[cfg_attr(docsrs, doc(cfg(feature = "chrono")))]
517 #[allow(clippy::needless_pass_by_value)]
518 pub fn estimate_from_date_like<D: Datelike>(date: D) -> Result<f64> {
519 Self::estimate_from_date(date.year(), date.month())
520 }
521}
522
523#[cfg(test)]
524mod tests {
525 use super::*;
526
527 const EPSILON: f64 = 1e-10;
528
529 #[test]
530 fn test_julian_date_creation() {
531 let jd = JulianDate::from_utc(2000, 1, 1, 12, 0, 0.0, 0.0).unwrap();
532
533 assert!((jd.julian_date() - J2000_JDN).abs() < EPSILON);
535 assert_eq!(jd.delta_t(), 0.0);
536 }
537
538 #[test]
539 fn test_julian_date_invalid_day_validation() {
540 assert!(JulianDate::from_utc(2024, 2, 30, 0, 0, 0.0, 0.0).is_err());
541 assert!(JulianDate::from_utc(2024, 2, 29, 0, 0, 0.0, 0.0).is_ok());
542 assert!(JulianDate::from_utc(1900, 2, 29, 0, 0, 0.0, 0.0).is_err());
543 assert!(JulianDate::from_utc(1500, 2, 29, 0, 0, 0.0, 0.0).is_ok());
544 assert!(JulianDate::from_utc(1582, 10, 10, 0, 0, 0.0, 0.0).is_err());
545 assert!(JulianDate::from_utc(1582, 10, 4, 0, 0, 0.0, 0.0).is_ok());
546 assert!(JulianDate::from_utc(1582, 10, 15, 0, 0, 0.0, 0.0).is_ok());
547 }
548
549 #[test]
550 fn test_julian_date_validation() {
551 assert!(JulianDate::from_utc(2024, 13, 1, 0, 0, 0.0, 0.0).is_err()); assert!(JulianDate::from_utc(2024, 1, 32, 0, 0, 0.0, 0.0).is_err()); assert!(JulianDate::from_utc(2024, 1, 1, 24, 0, 0.0, 0.0).is_err()); assert!(JulianDate::from_utc(2024, 1, 1, 0, 60, 0.0, 0.0).is_err()); assert!(JulianDate::from_utc(2024, 1, 1, 0, 0, 60.0, 0.0).is_err()); assert!(JulianDate::from_utc(2024, 1, 1, 0, 0, 0.0, f64::NAN).is_err()); assert!(JulianDate::from_utc(2024, 1, 1, 0, 0, 0.0, f64::INFINITY).is_err());
558 }
560
561 #[test]
562 fn test_julian_centuries() {
563 let jd = JulianDate::from_utc(2000, 1, 1, 12, 0, 0.0, 0.0).unwrap();
564
565 assert!(jd.julian_century().abs() < EPSILON);
567 assert!(jd.julian_ephemeris_century().abs() < EPSILON);
568 assert!(jd.julian_ephemeris_millennium().abs() < EPSILON);
569 }
570
571 #[test]
572 fn test_julian_ephemeris_day() {
573 let delta_t = 69.0; let jd = JulianDate::from_utc(2023, 6, 21, 12, 0, 0.0, delta_t).unwrap();
575
576 let jde = jd.julian_ephemeris_day();
577 let expected = jd.julian_date() + delta_t / SECONDS_PER_DAY;
578
579 assert!((jde - expected).abs() < EPSILON);
580 }
581
582 #[test]
583 fn test_gregorian_calendar_correction() {
584 let julian_date = JulianDate::from_utc(1582, 10, 4, 12, 0, 0.0, 0.0).unwrap();
587 let gregorian_date = JulianDate::from_utc(1582, 10, 15, 12, 0, 0.0, 0.0).unwrap();
588
589 let diff = gregorian_date.julian_date() - julian_date.julian_date();
592 assert!(
593 (diff - 1.0).abs() < 1e-6,
594 "Expected 1 day difference in JD, got {diff}"
595 );
596
597 let pre_gregorian = JulianDate::from_utc(1582, 10, 1, 12, 0, 0.0, 0.0).unwrap();
600 let post_gregorian = JulianDate::from_utc(1583, 1, 1, 12, 0, 0.0, 0.0).unwrap();
601
602 assert!(pre_gregorian.julian_date() > 2_000_000.0);
604 assert!(post_gregorian.julian_date() > pre_gregorian.julian_date());
605 }
606
607 #[test]
608 fn test_delta_t_modern_estimates() {
609 let delta_t_2000 = DeltaT::estimate(2000.0).unwrap();
611 let delta_t_2020 = DeltaT::estimate(2020.0).unwrap();
612
613 assert!(delta_t_2000 > 60.0 && delta_t_2000 < 70.0);
614 assert!(delta_t_2020 > 65.0 && delta_t_2020 < 75.0);
615 assert!(delta_t_2020 > delta_t_2000); }
617
618 #[test]
619 fn test_delta_t_historical_estimates() {
620 let delta_t_1900 = DeltaT::estimate(1900.0).unwrap();
621 let delta_t_1950 = DeltaT::estimate(1950.0).unwrap();
622
623 assert!(delta_t_1900 < 0.0); assert!(delta_t_1950 > 25.0 && delta_t_1950 < 35.0);
625 }
626
627 #[test]
628 fn test_delta_t_boundary_conditions() {
629 assert!(DeltaT::estimate(-500.0).is_ok());
631 assert!(DeltaT::estimate(3000.0).is_ok());
632 assert!(DeltaT::estimate(-501.0).is_err());
633 assert!(DeltaT::estimate(3001.0).is_err()); }
635
636 #[test]
637 fn test_delta_t_from_date() {
638 let delta_t = DeltaT::estimate_from_date(2024, 6).unwrap();
639 let delta_t_decimal = DeltaT::estimate(2024.5 - 1.0 / 24.0).unwrap(); assert!((delta_t - delta_t_decimal).abs() < 0.01);
643
644 assert!(DeltaT::estimate_from_date(2024, 13).is_err());
646 assert!(DeltaT::estimate_from_date(2024, 0).is_err());
647 }
648
649 #[test]
650 #[cfg(feature = "chrono")]
651 fn test_delta_t_from_date_like() {
652 use chrono::{DateTime, FixedOffset, NaiveDate, Utc};
653
654 let datetime_fixed = "2024-06-15T12:00:00-07:00"
656 .parse::<DateTime<FixedOffset>>()
657 .unwrap();
658 let delta_t_fixed = DeltaT::estimate_from_date_like(datetime_fixed).unwrap();
659
660 let datetime_utc = "2024-06-15T19:00:00Z".parse::<DateTime<Utc>>().unwrap();
662 let delta_t_utc = DeltaT::estimate_from_date_like(datetime_utc).unwrap();
663
664 let naive_date = NaiveDate::from_ymd_opt(2024, 6, 15).unwrap();
666 let delta_t_naive_date = DeltaT::estimate_from_date_like(naive_date).unwrap();
667
668 let naive_datetime = naive_date.and_hms_opt(12, 0, 0).unwrap();
670 let delta_t_naive_datetime = DeltaT::estimate_from_date_like(naive_datetime).unwrap();
671
672 assert_eq!(delta_t_fixed, delta_t_utc);
674 assert_eq!(delta_t_fixed, delta_t_naive_date);
675 assert_eq!(delta_t_fixed, delta_t_naive_datetime);
676
677 let delta_t_date = DeltaT::estimate_from_date(2024, 6).unwrap();
679 assert_eq!(delta_t_fixed, delta_t_date);
680
681 assert!(delta_t_fixed > 60.0 && delta_t_fixed < 80.0);
683 }
684
685 #[test]
686 fn test_specific_julian_dates() {
687 let unix_epoch = JulianDate::from_utc(1970, 1, 1, 0, 0, 0.0, 0.0).unwrap();
691 assert!((unix_epoch.julian_date() - 2_440_587.5).abs() < 1e-6);
692
693 let y2k = JulianDate::from_utc(2000, 1, 1, 0, 0, 0.0, 0.0).unwrap();
695 assert!((y2k.julian_date() - 2_451_544.5).abs() < 1e-6);
696 }
697}