open-standards/map/map-rs/crates/map-orbital/src/julian.rs

//! Julian date utilities for orbital mechanics.
//!
//! Provides conversions between calendar dates, Julian dates, and
//! Greenwich Mean Sidereal Time (GMST).

use chrono::{DateTime, Datelike, NaiveDate, NaiveDateTime, NaiveTime, TimeZone, Timelike, Utc};

use crate::constants::TWOPI;

/// Convert a `DateTime<Utc>` to Julian Date.
///
/// Uses the standard algorithm from Meeus, "Astronomical Algorithms".
pub fn datetime_to_jd(dt: &DateTime<Utc>) -> f64 {
    let year = dt.year();
    let month = dt.month() as i32;
    let day = dt.day() as f64;

    let hour = dt.hour() as f64;
    let minute = dt.minute() as f64;
    let second = dt.second() as f64 + dt.nanosecond() as f64 / 1_000_000_000.0;

    let frac = (hour + minute / 60.0 + second / 3600.0) / 24.0;

    let (y, m) = if month <= 2 {
        (year - 1, month + 12)
    } else {
        (year, month)
    };

    let a = (y as f64 / 100.0).floor();
    let b = 2.0 - a + (a / 4.0).floor();

    (365.25 * (y as f64 + 4716.0)).floor()
        + (30.6001 * (m as f64 + 1.0)).floor()
        + day
        + frac
        + b
        - 1524.5
}

/// Convert Julian Date back to `DateTime<Utc>`.
pub fn jd_to_datetime(jd: f64) -> DateTime<Utc> {
    let jd_shifted = jd + 0.5;
    let z = jd_shifted.floor();
    let f = jd_shifted - z;

    let a = if z < 2299161.0 {
        z
    } else {
        let alpha = ((z - 1867216.25) / 36524.25).floor();
        z + 1.0 + alpha - (alpha / 4.0).floor()
    };

    let b = a + 1524.0;
    let c = ((b - 122.1) / 365.25).floor();
    let d = (365.25 * c).floor();
    let e = ((b - d) / 30.6001).floor();

    let day_frac = b - d - (30.6001 * e).floor() + f;
    let day = day_frac.floor() as u32;
    let day_remainder = day_frac - day as f64;

    let month = if e < 14.0 {
        (e - 1.0) as u32
    } else {
        (e - 13.0) as u32
    };

    let year = if month > 2 {
        (c - 4716.0) as i32
    } else {
        (c - 4715.0) as i32
    };

    let total_seconds = (day_remainder * 86400.0 + 0.5).floor(); // round to nearest second
    let hour = (total_seconds / 3600.0).floor() as u32;
    let minute = ((total_seconds - hour as f64 * 3600.0) / 60.0).floor() as u32;
    let second_frac = total_seconds - hour as f64 * 3600.0 - minute as f64 * 60.0;
    let second = second_frac.floor() as u32;
    let nano = ((second_frac - second as f64) * 1_000_000_000.0).round() as u32;

    let date = NaiveDate::from_ymd_opt(year, month, day).unwrap_or(NaiveDate::from_ymd_opt(year, 1, 1).unwrap());
    let time = NaiveTime::from_hms_nano_opt(
        hour.min(23),
        minute.min(59),
        second.min(59),
        nano.min(999_999_999),
    )
    .unwrap();

    Utc.from_utc_datetime(&NaiveDateTime::new(date, time))
}

/// Convert TLE epoch (2-digit year + fractional day of year) to Julian Date.
pub fn epoch_to_jd(epoch_year: u16, epoch_day: f64) -> f64 {
    let year = epoch_year as i32;
    // January 1 of that year at 00:00 UTC
    let jan1 = NaiveDate::from_ymd_opt(year, 1, 1).unwrap();
    let dt = Utc.from_utc_datetime(
        &jan1.and_hms_opt(0, 0, 0).unwrap(),
    );
    let jd_jan1 = datetime_to_jd(&dt);
    // epoch_day is 1-based: day 1.0 = Jan 1 00:00
    jd_jan1 + epoch_day - 1.0
}

/// Convert TLE epoch to `DateTime<Utc>`.
pub fn epoch_to_datetime(epoch_year: u16, epoch_day: f64) -> DateTime<Utc> {
    jd_to_datetime(epoch_to_jd(epoch_year, epoch_day))
}

/// Greenwich Mean Sidereal Time in radians at the given Julian Date.
///
/// IAU 1982 expression from the Explanatory Supplement to the Astronomical Almanac.
pub fn gmst(jd: f64) -> f64 {
    // Julian centuries from J2000.0
    let t_ut1 = (jd - 2451545.0) / 36525.0;

    // GMST in seconds of time
    let gmst_sec = 67310.54841
        + (876600.0 * 3600.0 + 8640184.812866) * t_ut1
        + 0.093104 * t_ut1 * t_ut1
        - 6.2e-6 * t_ut1 * t_ut1 * t_ut1;

    // Convert to radians: seconds -> fraction of day -> radians
    // 86400 seconds in a day, one day = 2*PI radians of rotation
    let mut gmst_rad = (gmst_sec / 86400.0).rem_euclid(1.0) * TWOPI;

    if gmst_rad < 0.0 {
        gmst_rad += TWOPI;
    }

    gmst_rad
}

/// Greenwich Mean Sidereal Time from a `DateTime<Utc>`.
pub fn gmst_from_datetime(dt: &DateTime<Utc>) -> f64 {
    gmst(datetime_to_jd(dt))
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_j2000_jd() {
        // J2000.0 = January 1, 2000, 12:00 TT = JD 2451545.0
        let dt = Utc.with_ymd_and_hms(2000, 1, 1, 12, 0, 0).unwrap();
        let jd = datetime_to_jd(&dt);
        assert!((jd - 2451545.0).abs() < 1e-6, "J2000 JD: {jd}");
    }

    #[test]
    fn test_jd_roundtrip() {
        let dt = Utc.with_ymd_and_hms(2024, 6, 15, 14, 30, 0).unwrap();
        let jd = datetime_to_jd(&dt);
        let dt2 = jd_to_datetime(jd);
        assert_eq!(dt.year(), dt2.year());
        assert_eq!(dt.month(), dt2.month());
        assert_eq!(dt.day(), dt2.day());
        assert_eq!(dt.hour(), dt2.hour());
        // Allow 1-second tolerance due to floating-point JD precision
        let diff = (dt.timestamp() - dt2.timestamp()).abs();
        assert!(diff <= 1, "roundtrip time diff: {diff}s");
    }

    #[test]
    fn test_epoch_to_jd() {
        // 2024 day 1.5 = January 1, 2024 at 12:00 UTC
        let jd = epoch_to_jd(2024, 1.5);
        let dt = jd_to_datetime(jd);
        assert_eq!(dt.year(), 2024);
        assert_eq!(dt.month(), 1);
        assert_eq!(dt.day(), 1);
        assert_eq!(dt.hour(), 12);
    }

    #[test]
    fn test_gmst_j2000() {
        // At J2000.0, GMST should be approximately 280.46 degrees = 4.894 rad
        let gmst_val = gmst(2451545.0);
        let gmst_deg = gmst_val * 180.0 / std::f64::consts::PI;
        // GMST at J2000.0 is approximately 280.46 degrees
        assert!(
            (gmst_deg - 280.46).abs() < 1.0,
            "GMST at J2000.0: {gmst_deg} deg"
        );
    }
}