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
23#[derive(Debug, Clone, Copy, PartialEq)]
28pub struct JulianDate {
29 jd: f64,
31 delta_t: f64,
33}
34
35impl JulianDate {
36 #[cfg(feature = "chrono")]
50 #[cfg_attr(docsrs, doc(cfg(feature = "chrono")))]
51 pub fn from_datetime<Tz: TimeZone>(
52 datetime: &chrono::DateTime<Tz>,
53 delta_t: f64,
54 ) -> Result<Self> {
55 let utc_datetime = datetime.with_timezone(&chrono::Utc);
57 Self::from_utc(
58 utc_datetime.year(),
59 utc_datetime.month(),
60 utc_datetime.day(),
61 utc_datetime.hour(),
62 utc_datetime.minute(),
63 f64::from(utc_datetime.second()) + f64::from(utc_datetime.nanosecond()) / 1e9,
64 delta_t,
65 )
66 }
67
68 pub fn from_utc(
93 year: i32,
94 month: u32,
95 day: u32,
96 hour: u32,
97 minute: u32,
98 second: f64,
99 delta_t: f64,
100 ) -> Result<Self> {
101 if !(1..=12).contains(&month) {
103 return Err(Error::invalid_datetime("month must be between 1 and 12"));
104 }
105 if !(1..=31).contains(&day) {
106 return Err(Error::invalid_datetime("day must be between 1 and 31"));
107 }
108 if hour > 23 {
109 return Err(Error::invalid_datetime("hour must be between 0 and 23"));
110 }
111 if minute > 59 {
112 return Err(Error::invalid_datetime("minute must be between 0 and 59"));
113 }
114 if !(0.0..60.0).contains(&second) {
115 return Err(Error::invalid_datetime(
116 "second must be between 0 and 59.999...",
117 ));
118 }
119 if !delta_t.is_finite() {
120 return Err(Error::invalid_datetime("delta_t must be finite"));
121 }
122
123 if day > days_in_month(year, month, day)? {
124 return Err(Error::invalid_datetime("day is out of range for month"));
125 }
126
127 let jd = calculate_julian_date(year, month, day, hour, minute, second);
128 Ok(Self { jd, delta_t })
129 }
130
131 pub fn from_utc_simple(
147 year: i32,
148 month: u32,
149 day: u32,
150 hour: u32,
151 minute: u32,
152 second: f64,
153 ) -> Result<Self> {
154 Self::from_utc(year, month, day, hour, minute, second, 0.0)
155 }
156
157 #[must_use]
162 pub const fn julian_date(&self) -> f64 {
163 self.jd
164 }
165
166 #[must_use]
171 pub const fn delta_t(&self) -> f64 {
172 self.delta_t
173 }
174
175 #[must_use]
182 pub fn julian_ephemeris_day(&self) -> f64 {
183 self.jd + self.delta_t / SECONDS_PER_DAY
184 }
185
186 #[must_use]
193 pub fn julian_century(&self) -> f64 {
194 (self.jd - J2000_JDN) / DAYS_PER_CENTURY
195 }
196
197 #[must_use]
204 pub fn julian_ephemeris_century(&self) -> f64 {
205 (self.julian_ephemeris_day() - J2000_JDN) / DAYS_PER_CENTURY
206 }
207
208 #[must_use]
215 pub fn julian_ephemeris_millennium(&self) -> f64 {
216 self.julian_ephemeris_century() / 10.0
217 }
218
219 pub(crate) fn add_days(self, days: f64) -> Self {
221 Self {
222 jd: self.jd + days,
223 delta_t: self.delta_t,
224 }
225 }
226}
227
228fn calculate_julian_date(
233 year: i32,
234 month: u32,
235 day: u32,
236 hour: u32,
237 minute: u32,
238 second: f64,
239) -> f64 {
240 let mut y = year;
241 let mut m = i32::try_from(month).expect("month should be valid i32");
242
243 if m < 3 {
245 y -= 1;
246 m += 12;
247 }
248
249 let d = f64::from(day) + (f64::from(hour) + (f64::from(minute) + second / 60.0) / 60.0) / 24.0;
251
252 let mut jd =
254 floor(365.25 * (f64::from(y) + 4716.0)) + floor(30.6001 * f64::from(m + 1)) + d - 1524.5;
255
256 if jd >= 2_299_161.0 {
259 let a = floor(f64::from(y) / 100.0);
260 let b = 2.0 - a + floor(a / 4.0);
261 jd += b;
262 }
263
264 jd
265}
266
267const fn is_gregorian_date(year: i32, month: u32, day: u32) -> bool {
268 year > 1582 || (year == 1582 && (month > 10 || (month == 10 && day >= 15)))
269}
270
271const fn is_leap_year(year: i32, is_gregorian: bool) -> bool {
272 if is_gregorian {
273 (year % 4 == 0 && year % 100 != 0) || year % 400 == 0
274 } else {
275 year % 4 == 0
276 }
277}
278
279fn days_in_month(year: i32, month: u32, day: u32) -> Result<u32> {
280 if year == 1582 && month == 10 && (5..=14).contains(&day) {
281 return Err(Error::invalid_datetime(
282 "dates 1582-10-05 through 1582-10-14 do not exist in Gregorian calendar",
283 ));
284 }
285
286 let is_gregorian = is_gregorian_date(year, month, day);
287 let days = match month {
288 1 | 3 | 5 | 7 | 8 | 10 | 12 => 31,
289 4 | 6 | 9 | 11 => 30,
290 2 => {
291 if is_leap_year(year, is_gregorian) {
292 29
293 } else {
294 28
295 }
296 }
297 _ => unreachable!("month already validated"),
298 };
299 Ok(days)
300}
301
302pub struct DeltaT;
307
308impl DeltaT {
309 #[allow(clippy::too_many_lines)] pub fn estimate(decimal_year: f64) -> Result<f64> {
331 let year = decimal_year;
332
333 if !year.is_finite() {
334 return Err(Error::invalid_datetime("year must be finite"));
335 }
336
337 if year < -500.0 {
338 return Err(Error::invalid_datetime(
339 "ΔT estimates not available before year -500",
340 ));
341 }
342
343 let delta_t = if year < 500.0 {
344 let u = year / 100.0;
345 polynomial(
346 &[
347 10583.6,
348 -1014.41,
349 33.78311,
350 -5.952053,
351 -0.1798452,
352 0.022174192,
353 0.0090316521,
354 ],
355 u,
356 )
357 } else if year < 1600.0 {
358 let u = (year - 1000.0) / 100.0;
359 polynomial(
360 &[
361 1574.2,
362 -556.01,
363 71.23472,
364 0.319781,
365 -0.8503463,
366 -0.005050998,
367 0.0083572073,
368 ],
369 u,
370 )
371 } else if year < 1700.0 {
372 let t = year - 1600.0;
373 polynomial(&[120.0, -0.9808, -0.01532, 1.0 / 7129.0], t)
374 } else if year < 1800.0 {
375 let t = year - 1700.0;
376 polynomial(
377 &[8.83, 0.1603, -0.0059285, 0.00013336, -1.0 / 1_174_000.0],
378 t,
379 )
380 } else if year < 1860.0 {
381 let t = year - 1800.0;
382 polynomial(
383 &[
384 13.72,
385 -0.332447,
386 0.0068612,
387 0.0041116,
388 -0.00037436,
389 0.0000121272,
390 -0.0000001699,
391 0.000000000875,
392 ],
393 t,
394 )
395 } else if year < 1900.0 {
396 let t = year - 1860.0;
397 polynomial(
398 &[
399 7.62,
400 0.5737,
401 -0.251754,
402 0.01680668,
403 -0.0004473624,
404 1.0 / 233_174.0,
405 ],
406 t,
407 )
408 } else if year < 1920.0 {
409 let t = year - 1900.0;
410 polynomial(&[-2.79, 1.494119, -0.0598939, 0.0061966, -0.000197], t)
411 } else if year < 1941.0 {
412 let t = year - 1920.0;
413 polynomial(&[21.20, 0.84493, -0.076100, 0.0020936], t)
414 } else if year < 1961.0 {
415 let t = year - 1950.0;
416 polynomial(&[29.07, 0.407, -1.0 / 233.0, 1.0 / 2547.0], t)
417 } else if year < 1986.0 {
418 let t = year - 1975.0;
419 polynomial(&[45.45, 1.067, -1.0 / 260.0, -1.0 / 718.0], t)
420 } else if year < 2005.0 {
421 let t = year - 2000.0;
422 polynomial(
423 &[
424 63.86,
425 0.3345,
426 -0.060374,
427 0.0017275,
428 0.000651814,
429 0.00002373599,
430 ],
431 t,
432 )
433 } else if year < 2015.0 {
434 let t = year - 2005.0;
435 polynomial(&[64.69, 0.2930], t)
436 } else if year <= 3000.0 {
437 let t = year - 2015.0;
438 polynomial(&[67.62, 0.3645, 0.0039755], t)
439 } else {
440 return Err(Error::invalid_datetime(
441 "ΔT estimates not available beyond year 3000",
442 ));
443 };
444
445 Ok(delta_t)
446 }
447
448 pub fn estimate_from_date(year: i32, month: u32) -> Result<f64> {
465 if !(1..=12).contains(&month) {
466 return Err(Error::invalid_datetime("month must be between 1 and 12"));
467 }
468
469 let decimal_year = f64::from(year) + (f64::from(month) - 0.5) / 12.0;
470 Self::estimate(decimal_year)
471 }
472
473 #[cfg(feature = "chrono")]
502 #[cfg_attr(docsrs, doc(cfg(feature = "chrono")))]
503 #[allow(clippy::needless_pass_by_value)]
504 pub fn estimate_from_date_like<D: Datelike>(date: D) -> Result<f64> {
505 Self::estimate_from_date(date.year(), date.month())
506 }
507}
508
509#[cfg(test)]
510mod tests {
511 use super::*;
512
513 const EPSILON: f64 = 1e-10;
514
515 #[test]
516 fn test_julian_date_creation() {
517 let jd = JulianDate::from_utc(2000, 1, 1, 12, 0, 0.0, 0.0).unwrap();
518
519 assert!((jd.julian_date() - J2000_JDN).abs() < EPSILON);
521 assert_eq!(jd.delta_t(), 0.0);
522 }
523
524 #[test]
525 fn test_julian_date_invalid_day_validation() {
526 assert!(JulianDate::from_utc(2024, 2, 30, 0, 0, 0.0, 0.0).is_err());
527 assert!(JulianDate::from_utc(2024, 2, 29, 0, 0, 0.0, 0.0).is_ok());
528 assert!(JulianDate::from_utc(1900, 2, 29, 0, 0, 0.0, 0.0).is_err());
529 assert!(JulianDate::from_utc(1500, 2, 29, 0, 0, 0.0, 0.0).is_ok());
530 assert!(JulianDate::from_utc(1582, 10, 10, 0, 0, 0.0, 0.0).is_err());
531 assert!(JulianDate::from_utc(1582, 10, 4, 0, 0, 0.0, 0.0).is_ok());
532 assert!(JulianDate::from_utc(1582, 10, 15, 0, 0, 0.0, 0.0).is_ok());
533 }
534
535 #[test]
536 fn test_julian_date_validation() {
537 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());
544 }
546
547 #[test]
548 fn test_julian_centuries() {
549 let jd = JulianDate::from_utc(2000, 1, 1, 12, 0, 0.0, 0.0).unwrap();
550
551 assert!(jd.julian_century().abs() < EPSILON);
553 assert!(jd.julian_ephemeris_century().abs() < EPSILON);
554 assert!(jd.julian_ephemeris_millennium().abs() < EPSILON);
555 }
556
557 #[test]
558 fn test_julian_ephemeris_day() {
559 let delta_t = 69.0; let jd = JulianDate::from_utc(2023, 6, 21, 12, 0, 0.0, delta_t).unwrap();
561
562 let jde = jd.julian_ephemeris_day();
563 let expected = jd.julian_date() + delta_t / SECONDS_PER_DAY;
564
565 assert!((jde - expected).abs() < EPSILON);
566 }
567
568 #[test]
569 fn test_gregorian_calendar_correction() {
570 let julian_date = JulianDate::from_utc(1582, 10, 4, 12, 0, 0.0, 0.0).unwrap();
573 let gregorian_date = JulianDate::from_utc(1582, 10, 15, 12, 0, 0.0, 0.0).unwrap();
574
575 let diff = gregorian_date.julian_date() - julian_date.julian_date();
578 assert!(
579 (diff - 1.0).abs() < 1e-6,
580 "Expected 1 day difference in JD, got {diff}"
581 );
582
583 let pre_gregorian = JulianDate::from_utc(1582, 10, 1, 12, 0, 0.0, 0.0).unwrap();
586 let post_gregorian = JulianDate::from_utc(1583, 1, 1, 12, 0, 0.0, 0.0).unwrap();
587
588 assert!(pre_gregorian.julian_date() > 2_000_000.0);
590 assert!(post_gregorian.julian_date() > pre_gregorian.julian_date());
591 }
592
593 #[test]
594 fn test_delta_t_modern_estimates() {
595 let delta_t_2000 = DeltaT::estimate(2000.0).unwrap();
597 let delta_t_2020 = DeltaT::estimate(2020.0).unwrap();
598
599 assert!(delta_t_2000 > 60.0 && delta_t_2000 < 70.0);
600 assert!(delta_t_2020 > 65.0 && delta_t_2020 < 75.0);
601 assert!(delta_t_2020 > delta_t_2000); }
603
604 #[test]
605 fn test_delta_t_historical_estimates() {
606 let delta_t_1900 = DeltaT::estimate(1900.0).unwrap();
607 let delta_t_1950 = DeltaT::estimate(1950.0).unwrap();
608
609 assert!(delta_t_1900 < 0.0); assert!(delta_t_1950 > 25.0 && delta_t_1950 < 35.0);
611 }
612
613 #[test]
614 fn test_delta_t_boundary_conditions() {
615 assert!(DeltaT::estimate(-500.0).is_ok());
617 assert!(DeltaT::estimate(3000.0).is_ok());
618 assert!(DeltaT::estimate(-501.0).is_err());
619 assert!(DeltaT::estimate(3001.0).is_err()); }
621
622 #[test]
623 fn test_delta_t_from_date() {
624 let delta_t = DeltaT::estimate_from_date(2024, 6).unwrap();
625 let delta_t_decimal = DeltaT::estimate(2024.5 - 1.0 / 24.0).unwrap(); assert!((delta_t - delta_t_decimal).abs() < 0.01);
629
630 assert!(DeltaT::estimate_from_date(2024, 13).is_err());
632 assert!(DeltaT::estimate_from_date(2024, 0).is_err());
633 }
634
635 #[test]
636 #[cfg(feature = "chrono")]
637 fn test_delta_t_from_date_like() {
638 use chrono::{DateTime, FixedOffset, NaiveDate, Utc};
639
640 let datetime_fixed = "2024-06-15T12:00:00-07:00"
642 .parse::<DateTime<FixedOffset>>()
643 .unwrap();
644 let delta_t_fixed = DeltaT::estimate_from_date_like(datetime_fixed).unwrap();
645
646 let datetime_utc = "2024-06-15T19:00:00Z".parse::<DateTime<Utc>>().unwrap();
648 let delta_t_utc = DeltaT::estimate_from_date_like(datetime_utc).unwrap();
649
650 let naive_date = NaiveDate::from_ymd_opt(2024, 6, 15).unwrap();
652 let delta_t_naive_date = DeltaT::estimate_from_date_like(naive_date).unwrap();
653
654 let naive_datetime = naive_date.and_hms_opt(12, 0, 0).unwrap();
656 let delta_t_naive_datetime = DeltaT::estimate_from_date_like(naive_datetime).unwrap();
657
658 assert_eq!(delta_t_fixed, delta_t_utc);
660 assert_eq!(delta_t_fixed, delta_t_naive_date);
661 assert_eq!(delta_t_fixed, delta_t_naive_datetime);
662
663 let delta_t_date = DeltaT::estimate_from_date(2024, 6).unwrap();
665 assert_eq!(delta_t_fixed, delta_t_date);
666
667 assert!(delta_t_fixed > 60.0 && delta_t_fixed < 80.0);
669 }
670
671 #[test]
672 fn test_specific_julian_dates() {
673 let unix_epoch = JulianDate::from_utc(1970, 1, 1, 0, 0, 0.0, 0.0).unwrap();
677 assert!((unix_epoch.julian_date() - 2_440_587.5).abs() < 1e-6);
678
679 let y2k = JulianDate::from_utc(2000, 1, 1, 0, 0, 0.0, 0.0).unwrap();
681 assert!((y2k.julian_date() - 2_451_544.5).abs() < 1e-6);
682 }
683}