#![allow(clippy::unreadable_literal)]
#![allow(clippy::many_single_char_names)]
use crate::math::{floor, polynomial};
use crate::{Error, Result};
#[cfg(feature = "chrono")]
use chrono::{Datelike, TimeZone, Timelike};
const SECONDS_PER_DAY: f64 = 86_400.0;
const J2000_JDN: f64 = 2_451_545.0;
const DAYS_PER_CENTURY: f64 = 36_525.0;
pub(crate) fn validate_utc_components(
year: i32,
month: u32,
day: u32,
hour: u32,
minute: u32,
second: f64,
) -> Result<()> {
if !(1..=12).contains(&month) {
return Err(Error::InvalidDateTime {
message: "month must be between 1 and 12",
});
}
if !(1..=31).contains(&day) {
return Err(Error::InvalidDateTime {
message: "day must be between 1 and 31",
});
}
if hour > 23 {
return Err(Error::InvalidDateTime {
message: "hour must be between 0 and 23",
});
}
if minute > 59 {
return Err(Error::InvalidDateTime {
message: "minute must be between 0 and 59",
});
}
if !(0.0..60.0).contains(&second) {
return Err(Error::InvalidDateTime {
message: "second must be between 0 and 59.999...",
});
}
if year == 1582 && month == 10 && (5..=14).contains(&day) {
return Err(Error::InvalidDateTime {
message: "dates 1582-10-05 through 1582-10-14 do not exist in Gregorian calendar",
});
}
if day > days_in_month(year, month) {
return Err(Error::InvalidDateTime {
message: "day is out of range for month",
});
}
Ok(())
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct JulianDate {
jd: f64,
delta_t: f64,
}
impl JulianDate {
#[cfg(feature = "chrono")]
#[cfg_attr(docsrs, doc(cfg(feature = "chrono")))]
pub fn from_datetime<Tz: TimeZone>(
datetime: &chrono::DateTime<Tz>,
delta_t: f64,
) -> Result<Self> {
let utc_datetime = datetime.with_timezone(&chrono::Utc);
Self::from_utc(
utc_datetime.year(),
utc_datetime.month(),
utc_datetime.day(),
utc_datetime.hour(),
utc_datetime.minute(),
f64::from(utc_datetime.second()) + f64::from(utc_datetime.nanosecond()) / 1e9,
delta_t,
)
}
pub fn from_utc(
year: i32,
month: u32,
day: u32,
hour: u32,
minute: u32,
second: f64,
delta_t: f64,
) -> Result<Self> {
validate_utc_components(year, month, day, hour, minute, second)?;
if !delta_t.is_finite() {
return Err(Error::InvalidDateTime {
message: "delta_t must be finite",
});
}
let jd = calculate_julian_date(year, month, day, hour, minute, second);
Ok(Self { jd, delta_t })
}
pub fn from_utc_simple(
year: i32,
month: u32,
day: u32,
hour: u32,
minute: u32,
second: f64,
) -> Result<Self> {
Self::from_utc(year, month, day, hour, minute, second, 0.0)
}
#[must_use]
pub const fn julian_date(&self) -> f64 {
self.jd
}
#[must_use]
pub const fn delta_t(&self) -> f64 {
self.delta_t
}
#[must_use]
pub fn julian_ephemeris_day(&self) -> f64 {
self.jd + self.delta_t / SECONDS_PER_DAY
}
#[must_use]
pub fn julian_century(&self) -> f64 {
(self.jd - J2000_JDN) / DAYS_PER_CENTURY
}
#[must_use]
pub fn julian_ephemeris_century(&self) -> f64 {
(self.julian_ephemeris_day() - J2000_JDN) / DAYS_PER_CENTURY
}
#[must_use]
pub fn julian_ephemeris_millennium(&self) -> f64 {
self.julian_ephemeris_century() / 10.0
}
pub(crate) fn add_days(self, days: f64) -> Self {
Self {
jd: self.jd + days,
delta_t: self.delta_t,
}
}
}
fn calculate_julian_date(
year: i32,
month: u32,
day: u32,
hour: u32,
minute: u32,
second: f64,
) -> f64 {
let mut y = year;
let mut m = i32::try_from(month).expect("month should be valid i32");
if m < 3 {
y -= 1;
m += 12;
}
let d = f64::from(day) + (f64::from(hour) + (f64::from(minute) + second / 60.0) / 60.0) / 24.0;
let mut jd =
floor(365.25 * (f64::from(y) + 4716.0)) + floor(30.6001 * f64::from(m + 1)) + d - 1524.5;
if jd >= 2_299_161.0 {
let a = floor(f64::from(y) / 100.0);
let b = 2.0 - a + floor(a / 4.0);
jd += b;
}
jd
}
fn days_in_month(year: i32, month: u32) -> u32 {
let is_leap_year = if year > 1582 || (year == 1582 && month > 10) {
(year % 4 == 0 && year % 100 != 0) || year % 400 == 0
} else {
year % 4 == 0
};
match month {
1 | 3 | 5 | 7 | 8 | 10 | 12 => 31,
4 | 6 | 9 | 11 => 30,
2 if is_leap_year => 29,
2 => 28,
_ => unreachable!("month already validated"),
}
}
pub struct DeltaT;
impl DeltaT {
#[allow(clippy::too_many_lines)] pub fn estimate(decimal_year: f64) -> Result<f64> {
let year = decimal_year;
if !year.is_finite() {
return Err(Error::InvalidDateTime {
message: "year must be finite",
});
}
if year < -500.0 {
return Err(Error::InvalidDateTime {
message: "ΔT estimates not available before year -500",
});
}
let delta_t = if year < 500.0 {
let u = year / 100.0;
polynomial(
&[
10583.6,
-1014.41,
33.78311,
-5.952053,
-0.1798452,
0.022174192,
0.0090316521,
],
u,
)
} else if year < 1600.0 {
let u = (year - 1000.0) / 100.0;
polynomial(
&[
1574.2,
-556.01,
71.23472,
0.319781,
-0.8503463,
-0.005050998,
0.0083572073,
],
u,
)
} else if year < 1700.0 {
let t = year - 1600.0;
polynomial(&[120.0, -0.9808, -0.01532, 1.0 / 7129.0], t)
} else if year < 1800.0 {
let t = year - 1700.0;
polynomial(
&[8.83, 0.1603, -0.0059285, 0.00013336, -1.0 / 1_174_000.0],
t,
)
} else if year < 1860.0 {
let t = year - 1800.0;
polynomial(
&[
13.72,
-0.332447,
0.0068612,
0.0041116,
-0.00037436,
0.0000121272,
-0.0000001699,
0.000000000875,
],
t,
)
} else if year < 1900.0 {
let t = year - 1860.0;
polynomial(
&[
7.62,
0.5737,
-0.251754,
0.01680668,
-0.0004473624,
1.0 / 233_174.0,
],
t,
)
} else if year < 1920.0 {
let t = year - 1900.0;
polynomial(&[-2.79, 1.494119, -0.0598939, 0.0061966, -0.000197], t)
} else if year < 1941.0 {
let t = year - 1920.0;
polynomial(&[21.20, 0.84493, -0.076100, 0.0020936], t)
} else if year < 1961.0 {
let t = year - 1950.0;
polynomial(&[29.07, 0.407, -1.0 / 233.0, 1.0 / 2547.0], t)
} else if year < 1986.0 {
let t = year - 1975.0;
polynomial(&[45.45, 1.067, -1.0 / 260.0, -1.0 / 718.0], t)
} else if year < 2005.0 {
let t = year - 2000.0;
polynomial(
&[
63.86,
0.3345,
-0.060374,
0.0017275,
0.000651814,
0.00002373599,
],
t,
)
} else if year < 2015.0 {
let t = year - 2005.0;
polynomial(&[64.69, 0.2930], t)
} else if year <= 3000.0 {
let t = year - 2015.0;
polynomial(&[67.62, 0.3645, 0.0039755], t)
} else {
return Err(Error::InvalidDateTime {
message: "ΔT estimates not available beyond year 3000",
});
};
Ok(delta_t)
}
pub fn estimate_from_date(year: i32, month: u32) -> Result<f64> {
if !(1..=12).contains(&month) {
return Err(Error::InvalidDateTime {
message: "month must be between 1 and 12",
});
}
let decimal_year = f64::from(year) + (f64::from(month) - 0.5) / 12.0;
Self::estimate(decimal_year)
}
#[cfg(feature = "chrono")]
#[cfg_attr(docsrs, doc(cfg(feature = "chrono")))]
#[allow(clippy::needless_pass_by_value)]
pub fn estimate_from_date_like<D: Datelike>(date: D) -> Result<f64> {
Self::estimate_from_date(date.year(), date.month())
}
}
#[cfg(test)]
mod tests {
use super::*;
const EPSILON: f64 = 1e-10;
#[test]
fn test_julian_date_creation() {
let jd = JulianDate::from_utc(2000, 1, 1, 12, 0, 0.0, 0.0).unwrap();
assert!((jd.julian_date() - J2000_JDN).abs() < EPSILON);
assert_eq!(jd.delta_t(), 0.0);
}
#[test]
fn test_julian_date_invalid_day_validation() {
assert!(JulianDate::from_utc(2024, 2, 30, 0, 0, 0.0, 0.0).is_err());
assert!(JulianDate::from_utc(2024, 2, 29, 0, 0, 0.0, 0.0).is_ok());
assert!(JulianDate::from_utc(1900, 2, 29, 0, 0, 0.0, 0.0).is_err());
assert!(JulianDate::from_utc(1500, 2, 29, 0, 0, 0.0, 0.0).is_ok());
assert!(JulianDate::from_utc(1582, 10, 10, 0, 0, 0.0, 0.0).is_err());
assert!(JulianDate::from_utc(1582, 10, 4, 0, 0, 0.0, 0.0).is_ok());
assert!(JulianDate::from_utc(1582, 10, 15, 0, 0, 0.0, 0.0).is_ok());
}
#[test]
fn test_julian_date_validation() {
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());
}
#[test]
fn test_julian_centuries() {
let jd = JulianDate::from_utc(2000, 1, 1, 12, 0, 0.0, 0.0).unwrap();
assert!(jd.julian_century().abs() < EPSILON);
assert!(jd.julian_ephemeris_century().abs() < EPSILON);
assert!(jd.julian_ephemeris_millennium().abs() < EPSILON);
}
#[test]
fn test_julian_ephemeris_day() {
let delta_t = 69.0; let jd = JulianDate::from_utc(2023, 6, 21, 12, 0, 0.0, delta_t).unwrap();
let jde = jd.julian_ephemeris_day();
let expected = jd.julian_date() + delta_t / SECONDS_PER_DAY;
assert!((jde - expected).abs() < EPSILON);
}
#[test]
fn test_gregorian_calendar_correction() {
let julian_date = JulianDate::from_utc(1582, 10, 4, 12, 0, 0.0, 0.0).unwrap();
let gregorian_date = JulianDate::from_utc(1582, 10, 15, 12, 0, 0.0, 0.0).unwrap();
let diff = gregorian_date.julian_date() - julian_date.julian_date();
assert!(
(diff - 1.0).abs() < 1e-6,
"Expected 1 day difference in JD, got {diff}"
);
let pre_gregorian = JulianDate::from_utc(1582, 10, 1, 12, 0, 0.0, 0.0).unwrap();
let post_gregorian = JulianDate::from_utc(1583, 1, 1, 12, 0, 0.0, 0.0).unwrap();
assert!(pre_gregorian.julian_date() > 2_000_000.0);
assert!(post_gregorian.julian_date() > pre_gregorian.julian_date());
}
#[test]
fn test_delta_t_modern_estimates() {
let delta_t_2000 = DeltaT::estimate(2000.0).unwrap();
let delta_t_2020 = DeltaT::estimate(2020.0).unwrap();
assert!(delta_t_2000 > 60.0 && delta_t_2000 < 70.0);
assert!(delta_t_2020 > 65.0 && delta_t_2020 < 75.0);
assert!(delta_t_2020 > delta_t_2000); }
#[test]
fn test_delta_t_historical_estimates() {
let delta_t_1900 = DeltaT::estimate(1900.0).unwrap();
let delta_t_1950 = DeltaT::estimate(1950.0).unwrap();
assert!(delta_t_1900 < 0.0); assert!(delta_t_1950 > 25.0 && delta_t_1950 < 35.0);
}
#[test]
fn test_delta_t_boundary_conditions() {
assert!(DeltaT::estimate(-500.0).is_ok());
assert!(DeltaT::estimate(3000.0).is_ok());
assert!(DeltaT::estimate(-501.0).is_err());
assert!(DeltaT::estimate(3001.0).is_err()); }
#[test]
fn test_delta_t_from_date() {
let delta_t = DeltaT::estimate_from_date(2024, 6).unwrap();
let delta_t_decimal = DeltaT::estimate(2024.5 - 1.0 / 24.0).unwrap();
assert!((delta_t - delta_t_decimal).abs() < 0.01);
assert!(DeltaT::estimate_from_date(2024, 13).is_err());
assert!(DeltaT::estimate_from_date(2024, 0).is_err());
}
#[test]
#[cfg(feature = "chrono")]
fn test_delta_t_from_date_like() {
use chrono::{DateTime, FixedOffset, NaiveDate, Utc};
let datetime_fixed = "2024-06-15T12:00:00-07:00"
.parse::<DateTime<FixedOffset>>()
.unwrap();
let delta_t_fixed = DeltaT::estimate_from_date_like(datetime_fixed).unwrap();
let datetime_utc = "2024-06-15T19:00:00Z".parse::<DateTime<Utc>>().unwrap();
let delta_t_utc = DeltaT::estimate_from_date_like(datetime_utc).unwrap();
let naive_date = NaiveDate::from_ymd_opt(2024, 6, 15).unwrap();
let delta_t_naive_date = DeltaT::estimate_from_date_like(naive_date).unwrap();
let naive_datetime = naive_date.and_hms_opt(12, 0, 0).unwrap();
let delta_t_naive_datetime = DeltaT::estimate_from_date_like(naive_datetime).unwrap();
assert_eq!(delta_t_fixed, delta_t_utc);
assert_eq!(delta_t_fixed, delta_t_naive_date);
assert_eq!(delta_t_fixed, delta_t_naive_datetime);
let delta_t_date = DeltaT::estimate_from_date(2024, 6).unwrap();
assert_eq!(delta_t_fixed, delta_t_date);
assert!(delta_t_fixed > 60.0 && delta_t_fixed < 80.0);
}
#[test]
fn test_specific_julian_dates() {
let unix_epoch = JulianDate::from_utc(1970, 1, 1, 0, 0, 0.0, 0.0).unwrap();
assert!((unix_epoch.julian_date() - 2_440_587.5).abs() < 1e-6);
let y2k = JulianDate::from_utc(2000, 1, 1, 0, 0, 0.0, 0.0).unwrap();
assert!((y2k.julian_date() - 2_451_544.5).abs() < 1e-6);
}
}