use super::model::{Instant, InstantRepr, JulianDateSplit, TimeModelError, TimeScale};
use super::scales::julian_day_number;
use crate::astro::constants::time::SECONDS_PER_DAY_I64;
use crate::constants::{J2000_JD, SECONDS_PER_DAY};
pub const MJD_JD_OFFSET: f64 = 2_400_000.5;
pub const J2000_JULIAN_DAY_NUMBER: i64 = 2_451_545;
pub const J2000_NOON_OFFSET_S: i64 = 43_200;
#[must_use]
pub fn split_julian_date(
year: i32,
month: i32,
day: i32,
hour: i32,
minute: i32,
second: f64,
) -> (f64, f64) {
let jd_whole = julian_day_number(year, month, day) as f64 - 0.5;
let day_seconds = hour as f64 * 3600.0 + minute as f64 * 60.0 + second;
let fraction = day_seconds / SECONDS_PER_DAY;
(jd_whole, fraction)
}
impl Instant {
pub fn from_utc_civil(
year: i32,
month: i32,
day: i32,
hour: i32,
minute: i32,
second: f64,
) -> Result<Self, TimeModelError> {
let (jd_whole, fraction) = split_julian_date(year, month, day, hour, minute, second);
Ok(Self::from_julian_date(
TimeScale::Utc,
JulianDateSplit::new(jd_whole, fraction)?,
))
}
}
#[must_use]
pub fn j2000_seconds(year: i32, month: i32, day: i32, hour: i32, minute: i32, second: f64) -> f64 {
if !second.is_finite() {
return f64::NAN;
}
let whole_second = second.trunc();
let fraction = second - whole_second;
let noon_seconds_from_j2000 = (julian_day_number(year, month, day) - J2000_JD as i64)
.saturating_mul(SECONDS_PER_DAY as i64);
let day_seconds = (i64::from(hour) * 3600)
.saturating_add(i64::from(minute) * 60)
.saturating_add(whole_second as i64);
(noon_seconds_from_j2000
.saturating_sub(43_200)
.saturating_add(day_seconds)) as f64
+ fraction
}
#[must_use]
pub fn second_of_day(hour: i32, minute: i32, second: f64) -> f64 {
hour as f64 * 3600.0 + minute as f64 * 60.0 + second
}
#[must_use]
pub fn day_of_year(year: i32, month: i32, day: i32, hour: i32, minute: i32, second: f64) -> f64 {
let integer_day_of_year =
(julian_day_number(year, month, day) - julian_day_number(year, 1, 1) + 1) as f64;
let sod = second_of_day(hour, minute, second);
integer_day_of_year + sod / SECONDS_PER_DAY
}
#[must_use]
pub fn day_of_year_int(year: i32, month: i32, day: i32) -> i64 {
julian_day_number(year, month, day) - julian_day_number(year, 1, 1) + 1
}
#[must_use]
pub fn civil_from_julian_day_number(jdn: i64) -> (i64, i64, i64) {
let l = jdn + 68_569;
let n = 4 * l / 146_097;
let l = l - (146_097 * n + 3) / 4;
let i = 4_000 * (l + 1) / 1_461_001;
let l = l - 1_461 * i / 4 + 31;
let j = 80 * l / 2_447;
let day = l - 2_447 * j / 80;
let l = j / 11;
let month = j + 2 - 12 * l;
let year = 100 * (n - 49) + i + l;
(year, month, day)
}
#[must_use]
pub fn civil_from_j2000_seconds(seconds: i64) -> (i64, i64, i64, i64, i64, i64) {
let from_midnight = seconds + J2000_NOON_OFFSET_S;
let day_index = from_midnight.div_euclid(SECONDS_PER_DAY_I64);
let second_of_day = from_midnight.rem_euclid(SECONDS_PER_DAY_I64);
let (year, month, day) = civil_from_julian_day_number(day_index + J2000_JULIAN_DAY_NUMBER);
let hour = second_of_day / 3_600;
let minute = (second_of_day % 3_600) / 60;
let second = second_of_day % 60;
(year, month, day, hour, minute, second)
}
#[must_use]
pub fn mjd_from_jd(jd: f64) -> f64 {
jd - MJD_JD_OFFSET
}
#[must_use]
pub fn j2000_seconds_from_split(jd_whole: f64, fraction: f64) -> f64 {
(jd_whole - J2000_JD) * SECONDS_PER_DAY + fraction * SECONDS_PER_DAY
}
#[must_use]
pub fn seconds_between_splits(
later_whole: f64,
later_fraction: f64,
earlier_whole: f64,
earlier_fraction: f64,
) -> f64 {
((later_whole - earlier_whole) + (later_fraction - earlier_fraction)) * SECONDS_PER_DAY
}
#[must_use]
pub fn split_julian_date_from_j2000_seconds(seconds: i64) -> (f64, f64) {
let days = seconds.div_euclid(SECONDS_PER_DAY_I64);
let rem_s = seconds.rem_euclid(SECONDS_PER_DAY_I64);
(J2000_JD + days as f64, rem_s as f64 / SECONDS_PER_DAY)
}
#[must_use]
pub fn civil_from_split_julian_date(
jd_whole: f64,
fraction: f64,
) -> (i64, i64, i64, i64, i64, f64) {
debug_assert!(
(jd_whole.fract().abs() - 0.5).abs() < 1e-6,
"civil_from_split_julian_date expects a *.5 civil-midnight jd_whole, not a *.0 noon boundary"
);
let jdn = (jd_whole + 0.5).round() as i64;
let (year, month, day) = civil_from_julian_day_number(jdn);
let seconds_of_day = fraction * SECONDS_PER_DAY;
let hour = (seconds_of_day / 3600.0).floor() as i64;
let minute = ((seconds_of_day - hour as f64 * 3600.0) / 60.0).floor() as i64;
let second = seconds_of_day - hour as f64 * 3600.0 - minute as f64 * 60.0;
(year, month, day, hour, minute, second)
}
#[must_use]
pub fn split_julian_date_add_seconds(jd_whole: f64, fraction: f64, seconds: f64) -> (f64, f64) {
let mut whole = jd_whole;
let mut fraction = fraction + seconds / SECONDS_PER_DAY;
let carry = fraction.floor();
whole += carry;
fraction -= carry;
(whole, fraction)
}
#[must_use]
pub fn julian_date_from_instant(epoch: Instant) -> f64 {
match epoch.repr {
InstantRepr::JulianDate(split) => split.jd_whole + split.fraction,
InstantRepr::Nanos(nanos) => (nanos as f64) / 1.0e9 / SECONDS_PER_DAY + J2000_JD,
}
}
#[must_use]
pub fn fractional_day_of_year_from_instant(epoch: Instant) -> f64 {
let jd = julian_date_from_instant(epoch);
let jd_midnight = jd + 0.5;
let day_floor = jd_midnight.floor();
let day_fraction = jd_midnight - day_floor;
let (year, month, day) = civil_from_julian_day_number(day_floor as i64);
let doy_integer = day_of_year_int(year as i32, month as i32, day as i32);
doy_integer as f64 + day_fraction
}
#[must_use]
pub fn second_of_day_from_instant(epoch: Instant) -> f64 {
match epoch.repr {
InstantRepr::JulianDate(jd) => {
let from_midnight = jd.jd_whole + 0.5 + jd.fraction;
let day_fraction = from_midnight - from_midnight.floor();
day_fraction * SECONDS_PER_DAY
}
InstantRepr::Nanos(nanos) => {
let ns_per_day: i128 = 86_400 * 1_000_000_000;
let mut rem = nanos % ns_per_day;
if rem < 0 {
rem += ns_per_day;
}
rem as f64 / 1.0e9
}
}
}
#[must_use]
pub const fn is_leap_year(year: i64) -> bool {
(year % 4 == 0 && year % 100 != 0) || year % 400 == 0
}
#[must_use]
pub const fn days_in_month(year: i64, month: i64) -> i64 {
match month {
1 | 3 | 5 | 7 | 8 | 10 | 12 => 31,
4 | 6 | 9 | 11 => 30,
2 if is_leap_year(year) => 29,
2 => 28,
_ => 0,
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn split_jd_j2000_noon_is_exact() {
let (whole, frac) = split_julian_date(2000, 1, 1, 12, 0, 0.0);
assert_eq!(whole, 2_451_544.5);
assert_eq!(frac, 0.5);
assert_eq!(whole + frac, J2000_JD);
}
#[test]
fn j2000_seconds_epoch_is_zero() {
assert_eq!(j2000_seconds(2000, 1, 1, 12, 0, 0.0), 0.0);
assert_eq!(j2000_seconds(2000, 1, 2, 12, 0, 0.0), 86_400.0);
assert_eq!(j2000_seconds(2000, 1, 1, 12, 0, 0.25), 0.25);
}
#[test]
fn second_of_day_is_clock_arithmetic() {
assert_eq!(second_of_day(0, 0, 0.0), 0.0);
assert_eq!(second_of_day(1, 2, 3.5), 3723.5);
}
#[test]
fn day_of_year_jan1_midnight_is_one() {
assert_eq!(day_of_year(2021, 1, 1, 0, 0, 0.0), 1.0);
assert_eq!(day_of_year(2021, 1, 1, 12, 0, 0.0), 1.5);
assert_eq!(day_of_year(2020, 3, 1, 0, 0, 0.0), 61.0);
assert_eq!(day_of_year(2021, 3, 1, 0, 0, 0.0), 60.0);
}
#[test]
fn day_of_year_int_matches_fractional_and_cumulative_table() {
fn cum_table(year: i64, month: i64, day: i64) -> i64 {
const CUM: [i64; 12] = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334];
let leap = is_leap_year(year);
let mut doy = CUM[(month - 1) as usize] + day;
if leap && month > 2 {
doy += 1;
}
doy
}
let cases: [(i64, i64, i64); 7] = [
(2000, 1, 1),
(2000, 3, 1),
(2001, 3, 1),
(2020, 12, 31),
(1999, 6, 15),
(2100, 3, 1),
(2400, 2, 29),
];
for (y, m, d) in cases {
assert_eq!(
day_of_year_int(y as i32, m as i32, d as i32),
cum_table(y, m, d)
);
assert_eq!(
day_of_year_int(y as i32, m as i32, d as i32) as f64,
day_of_year(y as i32, m as i32, d as i32, 0, 0, 0.0)
);
}
}
#[test]
fn civil_from_jdn_inverts_julian_day_number() {
let cases: [(i32, i64, i64); 7] = [
(2000, 1, 1),
(1980, 1, 6),
(2006, 1, 1),
(1970, 1, 1),
(2024, 2, 29),
(1, 1, 1),
(2099, 12, 31),
];
for (y, m, d) in cases {
let jdn = julian_day_number(y, m as i32, d as i32);
assert_eq!(civil_from_julian_day_number(jdn), (i64::from(y), m, d));
}
for jdn in (2_400_000..2_500_000).step_by(37) {
let a = jdn + 32_044;
let b = (4 * a + 3) / 146_097;
let c = a - (146_097 * b) / 4;
let dd = (4 * c + 3) / 1461;
let e = c - (1461 * dd) / 4;
let mm = (5 * e + 2) / 153;
let day = e - (153 * mm + 2) / 5 + 1;
let month = mm + 3 - 12 * (mm / 10);
let year = 100 * b + dd - 4800 + mm / 10;
assert_eq!(civil_from_julian_day_number(jdn), (year, month, day));
}
}
#[test]
fn civil_from_j2000_seconds_inverts_j2000_seconds() {
assert_eq!(civil_from_j2000_seconds(0), (2000, 1, 1, 12, 0, 0));
assert_eq!(civil_from_j2000_seconds(86_400), (2000, 1, 2, 12, 0, 0));
assert_eq!(
civil_from_j2000_seconds(646_315_200),
(2020, 6, 25, 0, 0, 0)
);
let cases: [(i32, i64, i64, i64, i64, i64); 4] = [
(2000, 1, 1, 12, 0, 0),
(1999, 12, 31, 23, 59, 59),
(2024, 2, 29, 6, 30, 15),
(2030, 7, 4, 0, 0, 1),
];
for (y, m, d, h, mi, s) in cases {
let secs = j2000_seconds(y, m as i32, d as i32, h as i32, mi as i32, s as f64) as i64;
assert_eq!(
civil_from_j2000_seconds(secs),
(i64::from(y), m, d, h, mi, s)
);
}
}
#[test]
fn mjd_and_leap_and_month_helpers() {
assert_eq!(mjd_from_jd(J2000_JD), 51_544.5);
assert!(is_leap_year(2000) && is_leap_year(2024) && !is_leap_year(2100));
assert!(!is_leap_year(2023));
assert_eq!(days_in_month(2024, 2), 29);
assert_eq!(days_in_month(2023, 2), 28);
assert_eq!(days_in_month(2024, 4), 30);
assert_eq!(days_in_month(2024, 13), 0);
}
#[test]
fn split_from_j2000_seconds_inverts_field_form() {
assert_eq!(split_julian_date_from_j2000_seconds(0), (J2000_JD, 0.0));
assert_eq!(
split_julian_date_from_j2000_seconds(86_400),
(J2000_JD + 1.0, 0.0)
);
let (whole, frac) = split_julian_date_from_j2000_seconds(646_315_200);
assert_eq!(j2000_seconds_from_split(whole, frac), 646_315_200.0);
}
#[test]
fn civil_from_split_round_trips_split_julian_date() {
let cases: [(i32, i32, i32, i32, i32, f64); 4] = [
(2000, 1, 1, 12, 0, 0.0),
(2020, 6, 25, 0, 0, 0.0),
(2024, 2, 29, 6, 30, 15.0),
(1999, 12, 31, 23, 59, 59.0),
];
for (y, mo, d, h, mi, s) in cases {
let (whole, frac) = split_julian_date(y, mo, d, h, mi, s);
let civil = civil_from_split_julian_date(whole, frac);
assert_eq!(
civil,
(
i64::from(y),
i64::from(mo),
i64::from(d),
i64::from(h),
i64::from(mi),
s
)
);
}
}
#[test]
fn split_add_seconds_carries_across_midnight() {
let (w, f) = split_julian_date_add_seconds(2_451_544.5, 0.25, 6.0 * 3600.0);
assert_eq!((w, f), (2_451_544.5, 0.5));
let (w, f) = split_julian_date_add_seconds(2_451_544.5, 0.5, 18.0 * 3600.0);
assert_eq!(w, 2_451_545.5);
assert!((f - 0.25).abs() < 1e-12);
}
#[test]
fn instant_julian_date_and_doy_and_sod() {
use super::super::model::{JulianDateSplit, TimeScale};
let (whole, frac) = split_julian_date(2020, 6, 25, 6, 0, 0.0);
let epoch = Instant::from_julian_date(
TimeScale::Gpst,
JulianDateSplit::new(whole, frac).expect("valid split"),
);
assert_eq!(julian_date_from_instant(epoch), whole + frac);
assert!((fractional_day_of_year_from_instant(epoch) - (177.0 + 0.25)).abs() < 1e-9);
assert!((second_of_day_from_instant(epoch) - 6.0 * 3600.0).abs() < 1e-6);
}
#[test]
fn from_utc_civil_matches_open_coded_path_and_accessors() {
use super::super::model::{JulianDateSplit, TimeScale};
let cases: [(i32, i32, i32, i32, i32, f64); 4] = [
(2000, 1, 1, 12, 0, 0.0),
(2020, 6, 25, 6, 0, 0.0),
(2024, 2, 29, 6, 30, 15.0),
(1999, 12, 31, 23, 59, 59.0),
];
for (y, mo, d, h, mi, s) in cases {
let instant = Instant::from_utc_civil(y, mo, d, h, mi, s).expect("valid civil instant");
let (whole, frac) = split_julian_date(y, mo, d, h, mi, s);
let expected = Instant::from_julian_date(
TimeScale::Utc,
JulianDateSplit::new(whole, frac).expect("valid split"),
);
assert_eq!(instant, expected);
assert_eq!(instant.scale, TimeScale::Utc);
assert_eq!(
instant.julian_date(),
Some(JulianDateSplit {
jd_whole: whole,
fraction: frac
})
);
assert_eq!(julian_date_from_instant(instant), whole + frac);
let civil = civil_from_split_julian_date(whole, frac);
assert_eq!(
civil,
(
i64::from(y),
i64::from(mo),
i64::from(d),
i64::from(h),
i64::from(mi),
s
)
);
}
}
#[test]
fn from_utc_civil_rejects_out_of_day_clock_field() {
assert!(Instant::from_utc_civil(2020, 6, 25, 25, 0, 0.0).is_err());
}
#[test]
fn j2000_seconds_from_split_matches_field_form() {
let (whole, frac) = split_julian_date(2020, 6, 25, 0, 0, 0.0);
assert_eq!(
j2000_seconds_from_split(whole, frac),
j2000_seconds(2020, 6, 25, 0, 0, 0.0)
);
}
}