falak 1.0.0

//! Ephemeris computation — Julian date, sidereal time, planetary/lunar positions.
//!
//! Provides the time foundation for orbital mechanics: Julian dates,
//! Modified Julian dates, Greenwich Mean Sidereal Time, and simplified
//! planetary/lunar position models.

use tracing::instrument;

use crate::error::{FalakError, Result};

/// Julian date of the J2000.0 epoch (2000 January 1, 12:00 TT).
pub const J2000_JD: f64 = 2_451_545.0;

/// Julian date of the Unix epoch (1970 January 1, 00:00 UTC).
pub const UNIX_EPOCH_JD: f64 = 2_440_587.5;

/// Modified Julian Date offset (MJD = JD − 2400000.5).
pub const MJD_OFFSET: f64 = 2_400_000.5;

/// Seconds per Julian day.
pub const SECONDS_PER_DAY: f64 = 86_400.0;

/// Seconds per sidereal day (IERS 2010).
pub const SECONDS_PER_SIDEREAL_DAY: f64 = 86_164.090_5;

/// Earth's sidereal rotation rate (rad/s).
pub const EARTH_ROTATION_RATE: f64 = std::f64::consts::TAU / SECONDS_PER_SIDEREAL_DAY;

/// Julian century in days.
pub const DAYS_PER_JULIAN_CENTURY: f64 = 36_525.0;

// ── Julian Date conversions ───────────────────────────────────────────────

/// Convert a calendar date to Julian Date.
///
/// Uses the algorithm valid for dates after the Gregorian calendar reform (1582).
///
/// # Arguments
///
/// * `year` — Calendar year (e.g. 2024)
/// * `month` — Month (1–12)
/// * `day` — Day of month (1–31, may include fractional part for time)
///
/// # Errors
///
/// Returns [`FalakError::InvalidParameter`] if month is out of range.
#[must_use = "returns the computed Julian Date"]
#[instrument(level = "trace")]
pub fn calendar_to_jd(year: i32, month: u32, day: f64) -> Result<f64> {
    if !(1..=12).contains(&month) {
        return Err(FalakError::InvalidParameter(
            format!("month must be 1–12, got {month}").into(),
        ));
    }

    let (y, m) = if month <= 2 {
        (year as f64 - 1.0, month as f64 + 12.0)
    } else {
        (year as f64, month as f64)
    };

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

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

/// Convert Julian Date to calendar date.
///
/// Returns `(year, month, day)` where day may include a fractional part.
#[must_use]
pub fn jd_to_calendar(jd: f64) -> (i32, u32, f64) {
    let jd = jd + 0.5;
    let z = jd.floor();
    let f = jd - z;

    let a = if z < 2_299_161.0 {
        z
    } else {
        let alpha = ((z - 1_867_216.25) / 36_524.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 = b - d - (30.6001 * e).floor() + f;
    let month = if e < 14.0 { e - 1.0 } else { e - 13.0 };
    let year = if month > 2.0 { c - 4716.0 } else { c - 4715.0 };

    (year as i32, month as u32, day)
}

/// Convert Julian Date to Modified Julian Date.
#[must_use]
#[inline]
pub fn jd_to_mjd(jd: f64) -> f64 {
    jd - MJD_OFFSET
}

/// Convert Modified Julian Date to Julian Date.
#[must_use]
#[inline]
pub fn mjd_to_jd(mjd: f64) -> f64 {
    mjd + MJD_OFFSET
}

/// Convert Unix timestamp (seconds since 1970-01-01 00:00 UTC) to Julian Date.
#[must_use]
#[inline]
pub fn unix_to_jd(unix_seconds: f64) -> f64 {
    UNIX_EPOCH_JD + unix_seconds / SECONDS_PER_DAY
}

/// Convert Julian Date to Unix timestamp.
#[must_use]
#[inline]
pub fn jd_to_unix(jd: f64) -> f64 {
    (jd - UNIX_EPOCH_JD) * SECONDS_PER_DAY
}

/// Julian centuries since J2000.0.
#[must_use]
#[inline]
pub fn julian_centuries_since_j2000(jd: f64) -> f64 {
    (jd - J2000_JD) / DAYS_PER_JULIAN_CENTURY
}

// ── Sidereal Time ─────────────────────────────────────────────────────────

/// Compute Greenwich Mean Sidereal Time (GMST) in radians.
///
/// Uses the IAU 1982 expression for GMST at 0h UT1, then adds the
/// fractional UT1 day scaled by the sidereal/solar ratio.
///
/// # Arguments
///
/// * `jd_ut1` — Julian Date in UT1 time scale
#[must_use]
#[inline]
pub fn gmst(jd_ut1: f64) -> f64 {
    // Separate into 0h UT1 and fractional day
    let jd_0h = (jd_ut1 + 0.5).floor() - 0.5;
    let frac_day = jd_ut1 - jd_0h;

    // Julian centuries from J2000.0 to 0h UT1 (NOT the full JD)
    let t0 = (jd_0h - J2000_JD) / DAYS_PER_JULIAN_CENTURY;

    // GMST at 0h UT1 in seconds (IAU 1982)
    let gmst_0h_sec =
        24_110.548_41 + 8_640_184.812_866 * t0 + 0.093_104 * t0 * t0 - 6.2e-6 * t0 * t0 * t0;

    // Add fractional day scaled by sidereal/solar ratio
    let gmst_total_sec = gmst_0h_sec + frac_day * SECONDS_PER_DAY * 1.002_737_909_350_795;

    // Convert to radians and normalise to [0, 2π)
    let gmst_rad = (gmst_total_sec / SECONDS_PER_DAY) * std::f64::consts::TAU;
    gmst_rad.rem_euclid(std::f64::consts::TAU)
}

// ── Day of year ───────────────────────────────────────────────────────────

/// Compute the day of year (1–366) from calendar date.
///
/// # Errors
///
/// Returns [`FalakError::InvalidParameter`] if month or day is out of range.
#[must_use = "returns the computed day of year"]
#[instrument(level = "trace")]
pub fn day_of_year(year: i32, month: u32, day: u32) -> Result<u32> {
    if !(1..=12).contains(&month) {
        return Err(FalakError::InvalidParameter(
            format!("month must be 1–12, got {month}").into(),
        ));
    }
    let is_leap = (year % 4 == 0 && year % 100 != 0) || year % 400 == 0;
    let days_in_month: [u32; 12] = if is_leap {
        [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    } else {
        [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    };

    let max_day = days_in_month[month as usize - 1];
    if day == 0 || day > max_day {
        return Err(FalakError::InvalidParameter(
            format!("day must be 1–{max_day} for month {month}, got {day}").into(),
        ));
    }

    let mut doy = day;
    for &dim in &days_in_month[..(month as usize - 1)] {
        doy += dim;
    }

    Ok(doy)
}

// ── Simplified planetary positions (VSOP87 truncated) ─────────────────────

/// Ecliptic longitude and latitude of a planet (radians) plus distance (AU).
#[derive(Debug, Clone, PartialEq, serde::Serialize, serde::Deserialize)]
#[non_exhaustive]
pub struct PlanetaryPosition {
    /// Ecliptic longitude (radians).
    pub longitude: f64,
    /// Ecliptic latitude (radians).
    pub latitude: f64,
    /// Distance from the Sun (AU).
    pub distance: f64,
}

/// Planet identifier for simplified ephemeris.
#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
#[non_exhaustive]
pub enum Planet {
    /// Mercury.
    Mercury,
    /// Venus.
    Venus,
    /// Earth.
    Earth,
    /// Mars.
    Mars,
    /// Jupiter.
    Jupiter,
    /// Saturn.
    Saturn,
}

/// Compute a simplified planetary position (truncated VSOP87-like series).
///
/// Accuracy: ~1° in longitude for dates within ±50 years of J2000.
/// This is a low-order approximation suitable for visualization and
/// third-body perturbation estimates, not precision navigation.
///
/// # Arguments
///
/// * `planet` — Which planet to compute.
/// * `jd` — Julian Date.
#[must_use]
pub fn planetary_position(planet: Planet, jd: f64) -> PlanetaryPosition {
    let t = julian_centuries_since_j2000(jd);

    // Orbital elements at epoch + secular rates (simplified)
    // Format: (L0, L_rate, lon_peri0, lon_peri_rate, a_AU, e0, e_rate)
    // L = mean longitude (deg), ϖ = longitude of perihelion (deg)
    let (l0, l_rate, wp0, wp_rate, a, e0, e_rate) = match planet {
        Planet::Mercury => (
            252.251,
            149_472.675,
            77.456,
            0.160,
            0.387_098,
            0.205_630,
            0.000_02,
        ),
        Planet::Venus => (
            181.980, 58_517.816, 131.564, 0.009, 0.723_332, 0.006_773, -0.000_05,
        ),
        Planet::Earth => (
            100.464, 35_999.373, 102.937, 0.032, 1.000_000, 0.016_709, -0.000_04,
        ),
        Planet::Mars => (
            355.433, 19_140.299, 336.060, 0.443, 1.523_688, 0.093_405, 0.000_09,
        ),
        Planet::Jupiter => (
            34.351, 3_034.906, 14.331, 0.216, 5.202_561, 0.048_498, 0.000_16,
        ),
        Planet::Saturn => (
            50.077, 1_222.114, 93.057, 0.891, 9.554_909, 0.055_509, -0.000_35,
        ),
    };

    let mean_lon_deg = l0 + l_rate * t;
    let lon_peri_deg = wp0 + wp_rate * t;
    let ecc = e0 + e_rate * t;

    // Mean anomaly M = L - ϖ
    let m = (mean_lon_deg - lon_peri_deg)
        .to_radians()
        .rem_euclid(std::f64::consts::TAU);

    // Solve Kepler's equation
    let mut ea = m + ecc * m.sin();
    for _ in 0..10 {
        let f = ea - ecc * ea.sin() - m;
        let fp = 1.0 - ecc * ea.cos();
        if fp.abs() < 1e-30 || f.abs() < 1e-12 {
            break;
        }
        ea -= f / fp;
    }

    // True anomaly
    let factor = ((1.0 + ecc) / (1.0 - ecc)).sqrt();
    let nu = 2.0 * (factor * (ea / 2.0).tan()).atan();

    // Heliocentric distance
    let distance = a * (1.0 - ecc * ea.cos());

    // Ecliptic longitude = longitude of perihelion + true anomaly
    let longitude = (lon_peri_deg.to_radians() + nu).rem_euclid(std::f64::consts::TAU);

    PlanetaryPosition {
        longitude,
        latitude: 0.0, // simplified: ecliptic plane
        distance,
    }
}

/// Convert heliocentric ecliptic position to Cartesian `[x, y, z]` in AU.
///
/// Z is perpendicular to the ecliptic (zero for simplified positions).
#[must_use]
#[inline]
pub fn ecliptic_to_cartesian(pos: &PlanetaryPosition) -> [f64; 3] {
    let r = pos.distance;
    [
        r * pos.longitude.cos() * pos.latitude.cos(),
        r * pos.longitude.sin() * pos.latitude.cos(),
        r * pos.latitude.sin(),
    ]
}

// ── Simple lunar ephemeris ────────────────────────────────────────────────

/// Lunar position in geocentric ecliptic coordinates.
#[derive(Debug, Clone, PartialEq, serde::Serialize, serde::Deserialize)]
#[non_exhaustive]
pub struct LunarPosition {
    /// Ecliptic longitude (radians).
    pub longitude: f64,
    /// Ecliptic latitude (radians).
    pub latitude: f64,
    /// Geocentric distance (km).
    pub distance_km: f64,
}

/// Compute a simplified lunar position.
///
/// Uses a low-order series approximation. Accuracy: ~1° longitude, ~0.5° latitude,
/// ~1000 km distance. Suitable for tidal perturbation estimates and visualization.
///
/// # Arguments
///
/// * `jd` — Julian Date.
#[must_use]
pub fn lunar_position(jd: f64) -> LunarPosition {
    let t = julian_centuries_since_j2000(jd);

    // Fundamental arguments (degrees, then convert)
    // Mean longitude of the Moon
    let l_prime = 218.316_447_7 + 481_267.881_343_6 * t;
    // Mean anomaly of the Moon
    let m_moon = 134.963_396_4 + 477_198.867_505_5 * t;
    // Mean anomaly of the Sun
    let m_sun = 357.529_109_2 + 35_999.050_290_9 * t;
    // Mean elongation of the Moon
    let d = 297.850_192_1 + 445_267.111_403_4 * t;
    // Argument of latitude of the Moon
    let f = 93.272_095_0 + 483_202.017_523_3 * t;

    let m_moon_r = m_moon.to_radians();
    let m_sun_r = m_sun.to_radians();
    let d_r = d.to_radians();
    let f_r = f.to_radians();

    // Longitude perturbations (largest terms, degrees)
    let mut lon_pert = 6.289 * m_moon_r.sin();
    lon_pert += 1.274 * (2.0 * d_r - m_moon_r).sin();
    lon_pert += 0.658 * (2.0 * d_r).sin();
    lon_pert += 0.214 * (2.0 * m_moon_r).sin();
    lon_pert -= 0.186 * m_sun_r.sin();
    lon_pert -= 0.114 * (2.0 * f_r).sin();

    let longitude = (l_prime + lon_pert)
        .to_radians()
        .rem_euclid(std::f64::consts::TAU);

    // Latitude (largest terms)
    let mut lat_pert = 5.128 * f_r.sin();
    lat_pert += 0.281 * (m_moon_r + f_r).sin();
    lat_pert += 0.278 * (m_moon_r - f_r).sin();

    let latitude = lat_pert.to_radians();

    // Distance (km, mean + largest perturbations)
    let mut dist = 385_000.56;
    dist -= 20_905.36 * m_moon_r.cos();
    dist -= 3_699.11 * (2.0 * d_r - m_moon_r).cos();
    dist -= 2_955.97 * (2.0 * d_r).cos();

    LunarPosition {
        longitude,
        latitude,
        distance_km: dist,
    }
}

/// Convert lunar position to geocentric Cartesian `[x, y, z]` in metres.
#[must_use]
#[inline]
pub fn lunar_to_cartesian_metres(pos: &LunarPosition) -> [f64; 3] {
    let r = pos.distance_km * 1000.0;
    [
        r * pos.longitude.cos() * pos.latitude.cos(),
        r * pos.longitude.sin() * pos.latitude.cos(),
        r * pos.latitude.sin(),
    ]
}

// ── Rise / set / transit ─────────────────────────────────────────────────

/// Rise, transit, and set times for a celestial body.
#[derive(Debug, Clone, PartialEq, serde::Serialize, serde::Deserialize)]
#[non_exhaustive]
pub struct RiseTransitSet {
    /// Julian date of rise (body crosses above the horizon), or `None` if
    /// the body is circumpolar or never rises.
    pub rise: Option<f64>,
    /// Julian date of transit (body crosses the observer's meridian).
    pub transit: Option<f64>,
    /// Julian date of set (body crosses below the horizon), or `None` if
    /// the body is circumpolar or never rises.
    pub set: Option<f64>,
}

/// Compute rise, transit, and set times for a body at a given right ascension
/// and declination, as observed from a location on Earth.
///
/// Uses the algorithm from Meeus (1991) Chapter 15, which provides times
/// accurate to about 1 minute for slowly-moving objects (stars, Sun).
///
/// # Arguments
///
/// * `jd_0h` — Julian date at 0h UT on the day of interest
/// * `ra` — Apparent right ascension (radians, 0..2π)
/// * `dec` — Apparent declination (radians, −π/2..π/2)
/// * `observer_lat` — Observer geodetic latitude (radians, −π/2..π/2)
/// * `observer_lon` — Observer longitude (radians, east positive)
/// * `horizon_elev` — Horizon elevation angle (radians, typically −0.5667° for
///   standard atmospheric refraction, or 0.0 for geometric horizon)
///
/// # Errors
///
/// Returns [`FalakError::EphemerisError`] if the body is circumpolar (never
/// rises or never sets at this latitude).
#[must_use = "returns rise/transit/set times"]
#[instrument(level = "trace")]
pub fn rise_transit_set(
    jd_0h: f64,
    ra: f64,
    dec: f64,
    observer_lat: f64,
    observer_lon: f64,
    horizon_elev: f64,
) -> Result<RiseTransitSet> {
    // Hour angle at rise/set: cos(H₀) = (sin(h₀) - sin(φ)sin(δ)) / (cos(φ)cos(δ))
    let sin_h0 = horizon_elev.sin();
    let cos_lat = observer_lat.cos();
    let sin_lat = observer_lat.sin();
    let cos_dec = dec.cos();
    let sin_dec = dec.sin();

    let denom = cos_lat * cos_dec;

    if denom.abs() < 1e-15 {
        // Observer at pole or body at pole — special case
        return if sin_lat * sin_dec > 0.0 {
            // Same hemisphere → circumpolar (never sets)
            Ok(RiseTransitSet {
                rise: None,
                transit: Some(transit_time(jd_0h, ra, observer_lon)),
                set: None,
            })
        } else {
            // Opposite hemisphere → never rises
            Ok(RiseTransitSet {
                rise: None,
                transit: None,
                set: None,
            })
        };
    }

    let cos_h0 = (sin_h0 - sin_lat * sin_dec) / denom;

    if cos_h0 < -1.0 {
        // Body is circumpolar (always above horizon)
        return Ok(RiseTransitSet {
            rise: None,
            transit: Some(transit_time(jd_0h, ra, observer_lon)),
            set: None,
        });
    }
    if cos_h0 > 1.0 {
        // Body never rises
        return Ok(RiseTransitSet {
            rise: None,
            transit: None,
            set: None,
        });
    }

    let h0 = cos_h0.acos(); // hour angle at rise/set (radians)

    // GMST at 0h UT
    let gmst_0h = gmst(jd_0h);

    // Transit: when hour angle = 0 → local sidereal time = RA
    // m₀ = (RA - lon - GMST) / 2π  (fraction of day)
    let m0 = (ra - observer_lon - gmst_0h) / std::f64::consts::TAU;
    let m0 = m0.rem_euclid(1.0); // normalise to [0, 1)

    // Rise: m₁ = m₀ - H₀/(2π)
    let m1 = (m0 - h0 / std::f64::consts::TAU).rem_euclid(1.0);

    // Set: m₂ = m₀ + H₀/(2π)
    let m2 = (m0 + h0 / std::f64::consts::TAU).rem_euclid(1.0);

    Ok(RiseTransitSet {
        rise: Some(jd_0h + m1),
        transit: Some(jd_0h + m0),
        set: Some(jd_0h + m2),
    })
}

/// Compute transit time (meridian crossing) for a body.
fn transit_time(jd_0h: f64, ra: f64, observer_lon: f64) -> f64 {
    let gmst_0h = gmst(jd_0h);
    let m0 = (ra - observer_lon - gmst_0h) / std::f64::consts::TAU;
    jd_0h + m0.rem_euclid(1.0)
}

/// Standard atmospheric refraction correction for rise/set (radians).
///
/// The standard value is −0°34' = −0.5667°, which accounts for atmospheric
/// refraction at the horizon. Use this as `horizon_elev` in [`rise_transit_set`].
pub const STANDARD_REFRACTION: f64 = -0.009_890_199_5; // -0.5667° in radians

// ── Eclipse prediction ───────────────────────────────────────────────────

/// Eclipse illumination state of a satellite.
#[derive(Debug, Clone, Copy, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
#[non_exhaustive]
pub enum EclipseState {
    /// Fully illuminated by the Sun.
    Sunlit,
    /// In the penumbral shadow (partial illumination).
    Penumbra,
    /// In the umbral shadow (no direct sunlight).
    Umbra,
}

/// Result of an eclipse check.
#[derive(Debug, Clone, PartialEq, serde::Serialize, serde::Deserialize)]
#[non_exhaustive]
pub struct EclipseInfo {
    /// Current eclipse state.
    pub state: EclipseState,
    /// Shadow fraction: 0.0 = fully sunlit, 1.0 = full umbra.
    /// Values between 0 and 1 indicate penumbra.
    pub shadow_fraction: f64,
}

/// Check whether a satellite is in eclipse using a cylindrical shadow model.
///
/// The cylindrical model assumes the shadow is a cylinder of radius equal to
/// the occulting body's radius, cast in the anti-Sun direction. Simple and
/// fast; suitable for LEO/MEO eclipse estimation.
///
/// # Arguments
///
/// * `sat_pos` — Satellite position `[x, y, z]` (metres, ECI or body-centred inertial)
/// * `sun_pos` — Sun position `[x, y, z]` (metres, same frame)
/// * `body_radius` — Radius of the occulting body (metres)
///
/// # Returns
///
/// [`EclipseInfo`] with the eclipse state and shadow fraction.
///
/// **Limitation**: assumes parallel solar rays (Sun at infinity). Accurate for
/// LEO/MEO; use [`eclipse_conical`] for GEO or cislunar orbits.
#[must_use]
#[inline]
pub fn eclipse_cylindrical(sat_pos: [f64; 3], sun_pos: [f64; 3], body_radius: f64) -> EclipseInfo {
    // Unit vector from body centre toward the Sun
    let sun_mag =
        (sun_pos[0] * sun_pos[0] + sun_pos[1] * sun_pos[1] + sun_pos[2] * sun_pos[2]).sqrt();
    if sun_mag < 1e-10 {
        return EclipseInfo {
            state: EclipseState::Sunlit,
            shadow_fraction: 0.0,
        };
    }
    let sun_hat = [
        sun_pos[0] / sun_mag,
        sun_pos[1] / sun_mag,
        sun_pos[2] / sun_mag,
    ];

    // Project satellite position onto Sun direction
    let sat_dot_sun = sat_pos[0] * sun_hat[0] + sat_pos[1] * sun_hat[1] + sat_pos[2] * sun_hat[2];

    // If satellite is on the Sun side of the body, it's sunlit
    if sat_dot_sun > 0.0 {
        return EclipseInfo {
            state: EclipseState::Sunlit,
            shadow_fraction: 0.0,
        };
    }

    // Perpendicular distance from satellite to the Sun-body line
    let perp = [
        sat_pos[0] - sat_dot_sun * sun_hat[0],
        sat_pos[1] - sat_dot_sun * sun_hat[1],
        sat_pos[2] - sat_dot_sun * sun_hat[2],
    ];
    let perp_dist = (perp[0] * perp[0] + perp[1] * perp[1] + perp[2] * perp[2]).sqrt();

    if perp_dist < body_radius {
        EclipseInfo {
            state: EclipseState::Umbra,
            shadow_fraction: 1.0,
        }
    } else {
        EclipseInfo {
            state: EclipseState::Sunlit,
            shadow_fraction: 0.0,
        }
    }
}

/// Check whether a satellite is in eclipse using a conical shadow model.
///
/// The conical model accounts for the finite angular sizes of both the Sun
/// and the occulting body, producing accurate umbra/penumbra boundaries.
/// More accurate than cylindrical for high-altitude orbits (GEO, cislunar).
///
/// # Arguments
///
/// * `sat_pos` — Satellite position `[x, y, z]` (metres, body-centred inertial)
/// * `sun_pos` — Sun position `[x, y, z]` (metres, same frame)
/// * `body_radius` — Radius of the occulting body (metres)
/// * `sun_radius` — Radius of the Sun (metres, default 6.957e8)
#[must_use]
#[inline]
pub fn eclipse_conical(
    sat_pos: [f64; 3],
    sun_pos: [f64; 3],
    body_radius: f64,
    sun_radius: f64,
) -> EclipseInfo {
    let sat_mag =
        (sat_pos[0] * sat_pos[0] + sat_pos[1] * sat_pos[1] + sat_pos[2] * sat_pos[2]).sqrt();
    let sun_mag =
        (sun_pos[0] * sun_pos[0] + sun_pos[1] * sun_pos[1] + sun_pos[2] * sun_pos[2]).sqrt();

    if sat_mag < 1e-10 || sun_mag < 1e-10 {
        return EclipseInfo {
            state: EclipseState::Sunlit,
            shadow_fraction: 0.0,
        };
    }

    // Apparent angular radii as seen from the satellite
    let theta_body = (body_radius / sat_mag).asin(); // angular radius of occulting body
    let theta_sun = (sun_radius / sun_mag).asin(); // angular radius of Sun

    // Angle between satellite-to-centre and satellite-to-sun directions
    // sat_to_sun = sun_pos - sat_pos... no, we want body-sun angle as seen from sat.
    // The body is at origin. The angle is between -sat_pos (toward body) and
    // (sun_pos - sat_pos) (toward sun).
    let to_sun = [
        sun_pos[0] - sat_pos[0],
        sun_pos[1] - sat_pos[1],
        sun_pos[2] - sat_pos[2],
    ];
    let to_sun_mag = (to_sun[0] * to_sun[0] + to_sun[1] * to_sun[1] + to_sun[2] * to_sun[2]).sqrt();

    if to_sun_mag < 1e-10 {
        return EclipseInfo {
            state: EclipseState::Sunlit,
            shadow_fraction: 0.0,
        };
    }

    // Angle between "toward body centre" and "toward Sun" as seen from satellite
    let cos_sep = -(sat_pos[0] * to_sun[0] + sat_pos[1] * to_sun[1] + sat_pos[2] * to_sun[2])
        / (sat_mag * to_sun_mag);
    let sep = cos_sep.clamp(-1.0, 1.0).acos();

    if sep > theta_body + theta_sun {
        // No overlap — fully sunlit
        EclipseInfo {
            state: EclipseState::Sunlit,
            shadow_fraction: 0.0,
        }
    } else if sep < theta_body - theta_sun && theta_body > theta_sun {
        // Sun fully behind body — total eclipse (umbra)
        EclipseInfo {
            state: EclipseState::Umbra,
            shadow_fraction: 1.0,
        }
    } else if sep < theta_sun - theta_body && theta_sun > theta_body {
        // Body in front of Sun but smaller — annular eclipse
        // Shadow fraction is the area ratio
        let frac = (theta_body / theta_sun).powi(2);
        EclipseInfo {
            state: EclipseState::Penumbra,
            shadow_fraction: frac,
        }
    } else {
        // Partial overlap — penumbra
        // Approximate shadow fraction from overlap geometry
        let frac = overlap_fraction(sep, theta_body, theta_sun);
        EclipseInfo {
            state: EclipseState::Penumbra,
            shadow_fraction: frac,
        }
    }
}

/// Approximate fractional overlap area of two circles (discs on the sky).
/// Fraction of the Sun's disc (angular radius `r2`) obscured by the body
/// (angular radius `r1`) at angular separation `sep`.
///
/// This equals the fraction of sunlight blocked — a satellite in this
/// penumbral zone receives `(1 - overlap_fraction)` of full illumination.
fn overlap_fraction(sep: f64, r1: f64, r2: f64) -> f64 {
    if sep >= r1 + r2 {
        return 0.0;
    }
    if sep <= (r1 - r2).abs() {
        return r1.min(r2).powi(2) / r2.powi(2);
    }

    // Area of intersection of two circles
    let d = sep;
    let cos_a1 = (d * d + r1 * r1 - r2 * r2) / (2.0 * d * r1);
    let cos_a2 = (d * d + r2 * r2 - r1 * r1) / (2.0 * d * r2);
    let a1 = cos_a1.clamp(-1.0, 1.0).acos();
    let a2 = cos_a2.clamp(-1.0, 1.0).acos();

    let overlap = r1 * r1 * (a1 - a1.sin() * a1.cos()) + r2 * r2 * (a2 - a2.sin() * a2.cos());
    let sun_area = std::f64::consts::PI * r2 * r2;

    if sun_area > 0.0 {
        (overlap / sun_area).clamp(0.0, 1.0)
    } else {
        0.0
    }
}

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

    // ── Julian Date ──────────────────────────────────────────────────

    #[test]
    fn j2000_epoch() {
        // J2000.0 = 2000 Jan 1, 12:00 TT → JD 2451545.0
        let jd = calendar_to_jd(2000, 1, 1.5).unwrap();
        assert!(
            (jd - J2000_JD).abs() < 1e-6,
            "J2000: {} vs {}",
            jd,
            J2000_JD
        );
    }

    #[test]
    fn known_date_sputnik() {
        // Sputnik launch: 1957 Oct 4 → JD ~2436116.0
        let jd = calendar_to_jd(1957, 10, 4.0).unwrap();
        assert!((jd - 2_436_115.5).abs() < 0.5, "Sputnik JD: {}", jd);
    }

    #[test]
    fn calendar_roundtrip() {
        let jd = calendar_to_jd(2024, 7, 15.75).unwrap();
        let (y, m, d) = jd_to_calendar(jd);
        assert_eq!(y, 2024);
        assert_eq!(m, 7);
        assert!((d - 15.75).abs() < 1e-10, "day: {d}");
    }

    #[test]
    fn calendar_roundtrip_leap() {
        let jd = calendar_to_jd(2024, 2, 29.0).unwrap();
        let (y, m, d) = jd_to_calendar(jd);
        assert_eq!(y, 2024);
        assert_eq!(m, 2);
        assert!((d - 29.0).abs() < 1e-10);
    }

    #[test]
    fn calendar_invalid_month() {
        assert!(calendar_to_jd(2024, 0, 1.0).is_err());
        assert!(calendar_to_jd(2024, 13, 1.0).is_err());
    }

    // ── MJD ──────────────────────────────────────────────────────────

    #[test]
    fn mjd_roundtrip() {
        let jd = 2_460_000.5;
        let mjd = jd_to_mjd(jd);
        let jd2 = mjd_to_jd(mjd);
        assert!((jd - jd2).abs() < 1e-15);
    }

    // ── Unix ─────────────────────────────────────────────────────────

    #[test]
    fn unix_epoch() {
        let jd = unix_to_jd(0.0);
        assert!((jd - UNIX_EPOCH_JD).abs() < 1e-10, "unix epoch JD: {}", jd);
    }

    #[test]
    fn unix_roundtrip() {
        let ts = 1_700_000_000.0; // ~2023-11-14
        let jd = unix_to_jd(ts);
        let ts2 = jd_to_unix(jd);
        // Large values lose some f64 precision
        assert!((ts - ts2).abs() < 0.01, "unix roundtrip: {ts} vs {ts2}");
    }

    // ── GMST ─────────────────────────────────────────────────────────

    #[test]
    fn gmst_range() {
        // GMST should always be in [0, 2π)
        for day_offset in 0..365 {
            let jd = J2000_JD + day_offset as f64;
            let g = gmst(jd);
            assert!(
                (0.0..std::f64::consts::TAU).contains(&g),
                "GMST out of range at JD {jd}: {g}"
            );
        }
    }

    #[test]
    fn gmst_monotonic_over_day() {
        // GMST should increase over the course of a day (modulo 2π wrap)
        let jd_base = J2000_JD + 100.0; // arbitrary day
        let g1 = gmst(jd_base);
        let g2 = gmst(jd_base + 0.25); // 6 hours later
        // 6 hours ≈ π/2 radians of sidereal rotation
        let diff = (g2 - g1).rem_euclid(std::f64::consts::TAU);
        assert!(
            (diff - std::f64::consts::FRAC_PI_2).abs() < 0.05,
            "6h should ≈ π/2 rad: diff={diff}"
        );
    }

    #[test]
    fn gmst_j2000() {
        // At J2000.0 (2000 Jan 1.5 UT1), GMST = 18h 41m 50.55s = 280.4606°
        // IAU 1982: at 0h UT1 on Jan 1 2000, GMST = 24110.54841s = 6h 41m 50.5s
        // Plus 12h of sidereal rotation for the half-day → 280.46°
        let g = gmst(J2000_JD);
        let g_deg = g.to_degrees();
        assert!((g_deg - 280.46).abs() < 0.1, "GMST at J2000: {g_deg}°");
    }

    // ── Day of year ──────────────────────────────────────────────────

    #[test]
    fn doy_jan1() {
        assert_eq!(day_of_year(2024, 1, 1).unwrap(), 1);
    }

    #[test]
    fn doy_leap_dec31() {
        assert_eq!(day_of_year(2024, 12, 31).unwrap(), 366);
    }

    #[test]
    fn doy_non_leap_dec31() {
        assert_eq!(day_of_year(2023, 12, 31).unwrap(), 365);
    }

    #[test]
    fn doy_march1_leap() {
        // 2024 is leap: Jan(31) + Feb(29) + 1 = 61
        assert_eq!(day_of_year(2024, 3, 1).unwrap(), 61);
    }

    #[test]
    fn doy_invalid() {
        assert!(day_of_year(2024, 0, 1).is_err());
        assert!(day_of_year(2024, 1, 0).is_err());
    }

    #[test]
    fn doy_feb30_rejected() {
        assert!(day_of_year(2024, 2, 30).is_err());
        // But Feb 29 is valid in leap year
        assert!(day_of_year(2024, 2, 29).is_ok());
        // And rejected in non-leap year
        assert!(day_of_year(2023, 2, 29).is_err());
    }

    #[test]
    fn doy_apr31_rejected() {
        assert!(day_of_year(2024, 4, 31).is_err());
        assert!(day_of_year(2024, 4, 30).is_ok());
    }

    // ── Julian centuries ─────────────────────────────────────────────

    #[test]
    fn centuries_at_j2000() {
        assert!((julian_centuries_since_j2000(J2000_JD)).abs() < 1e-15);
    }

    #[test]
    fn centuries_one_century() {
        let jd = J2000_JD + DAYS_PER_JULIAN_CENTURY;
        assert!((julian_centuries_since_j2000(jd) - 1.0).abs() < 1e-12);
    }

    // ── Planetary positions ──────────────────────────────────────────

    #[test]
    fn earth_distance_1au() {
        let pos = planetary_position(Planet::Earth, J2000_JD);
        assert!(
            (pos.distance - 1.0).abs() < 0.02,
            "Earth distance: {} AU",
            pos.distance
        );
    }

    #[test]
    fn mercury_closer_than_earth() {
        let m = planetary_position(Planet::Mercury, J2000_JD);
        let e = planetary_position(Planet::Earth, J2000_JD);
        assert!(m.distance < e.distance, "Mercury should be closer to Sun");
    }

    #[test]
    fn jupiter_further_than_mars() {
        let j = planetary_position(Planet::Jupiter, J2000_JD);
        let m = planetary_position(Planet::Mars, J2000_JD);
        assert!(
            j.distance > m.distance,
            "Jupiter should be further than Mars"
        );
    }

    #[test]
    fn planet_longitude_range() {
        for planet in [
            Planet::Mercury,
            Planet::Venus,
            Planet::Earth,
            Planet::Mars,
            Planet::Jupiter,
            Planet::Saturn,
        ] {
            let pos = planetary_position(planet, J2000_JD + 500.0);
            assert!(
                (0.0..std::f64::consts::TAU).contains(&pos.longitude),
                "{planet:?} longitude out of range: {}",
                pos.longitude
            );
        }
    }

    #[test]
    fn ecliptic_to_cartesian_roundtrip() {
        let pos = planetary_position(Planet::Earth, J2000_JD);
        let cart = ecliptic_to_cartesian(&pos);
        let r = (cart[0] * cart[0] + cart[1] * cart[1] + cart[2] * cart[2]).sqrt();
        assert!(
            (r - pos.distance).abs() < 1e-10,
            "cartesian distance: {r} vs {}",
            pos.distance
        );
    }

    // ── Lunar position ───────────────────────────────────────────────

    #[test]
    fn lunar_distance_range() {
        // Moon distance varies ~356,500–406,700 km
        let pos = lunar_position(J2000_JD);
        assert!(
            pos.distance_km > 350_000.0 && pos.distance_km < 410_000.0,
            "lunar distance: {} km",
            pos.distance_km
        );
    }

    #[test]
    fn lunar_latitude_range() {
        // Lunar latitude should be within ±5.3°
        let pos = lunar_position(J2000_JD);
        assert!(
            pos.latitude.abs() < 6.0_f64.to_radians(),
            "lunar latitude: {}°",
            pos.latitude.to_degrees()
        );
    }

    #[test]
    fn lunar_cartesian_distance() {
        let pos = lunar_position(J2000_JD);
        let cart = lunar_to_cartesian_metres(&pos);
        let r = (cart[0] * cart[0] + cart[1] * cart[1] + cart[2] * cart[2]).sqrt();
        assert!(
            (r - pos.distance_km * 1000.0).abs() < 1.0,
            "cartesian: {r} vs {}",
            pos.distance_km * 1000.0
        );
    }

    #[test]
    fn lunar_position_varies() {
        let p1 = lunar_position(J2000_JD);
        let p2 = lunar_position(J2000_JD + 15.0); // half a lunar month
        // Longitude should differ significantly
        let diff = (p2.longitude - p1.longitude).rem_euclid(std::f64::consts::TAU);
        assert!(
            diff > 2.0,
            "lunar longitude should change over 15 days: {diff} rad"
        );
    }

    // ── Eclipse ──────────────────────────────────────────────────────

    const R_EARTH: f64 = 6_378_137.0;
    const R_SUN: f64 = 6.957e8;
    const AU: f64 = 1.495_978_707e11;

    #[test]
    fn eclipse_cylindrical_sunlit() {
        // Satellite on Sun side of Earth — clearly sunlit
        let sat = [R_EARTH + 400e3, 0.0, 0.0];
        let sun = [AU, 0.0, 0.0]; // Sun along +x
        let info = eclipse_cylindrical(sat, sun, R_EARTH);
        assert_eq!(info.state, EclipseState::Sunlit);
        assert!((info.shadow_fraction - 0.0).abs() < 1e-10);
    }

    #[test]
    fn eclipse_cylindrical_umbra() {
        // Satellite directly behind Earth (opposite Sun)
        let sat = [-(R_EARTH + 400e3), 0.0, 0.0];
        let sun = [AU, 0.0, 0.0];
        let info = eclipse_cylindrical(sat, sun, R_EARTH);
        assert_eq!(info.state, EclipseState::Umbra);
        assert!((info.shadow_fraction - 1.0).abs() < 1e-10);
    }

    #[test]
    fn eclipse_cylindrical_edge() {
        // Satellite behind Earth but far off-axis — should be sunlit
        let sat = [-1e7, R_EARTH * 2.0, 0.0]; // way off to the side
        let sun = [AU, 0.0, 0.0];
        let info = eclipse_cylindrical(sat, sun, R_EARTH);
        assert_eq!(info.state, EclipseState::Sunlit);
    }

    #[test]
    fn eclipse_conical_sunlit() {
        let sat = [R_EARTH + 400e3, 0.0, 0.0];
        let sun = [AU, 0.0, 0.0];
        let info = eclipse_conical(sat, sun, R_EARTH, R_SUN);
        assert_eq!(info.state, EclipseState::Sunlit);
    }

    #[test]
    fn eclipse_conical_umbra() {
        // Satellite directly behind Earth, close enough for umbra
        let sat = [-(R_EARTH + 400e3), 0.0, 0.0];
        let sun = [AU, 0.0, 0.0];
        let info = eclipse_conical(sat, sun, R_EARTH, R_SUN);
        // At LEO altitude, Earth's angular size > Sun's → full eclipse
        assert_eq!(info.state, EclipseState::Umbra);
        assert!((info.shadow_fraction - 1.0).abs() < 1e-10);
    }

    #[test]
    fn eclipse_conical_geo_shadow() {
        // At GEO (~42164 km), Earth's shadow is narrower — test that
        // a satellite directly behind Earth is still in shadow
        let r_geo = 42_164e3;
        let sat = [-r_geo, 0.0, 0.0];
        let sun = [AU, 0.0, 0.0];
        let info = eclipse_conical(sat, sun, R_EARTH, R_SUN);
        // At GEO, Earth (angular radius ~8.7°) still bigger than Sun (~0.27°)
        assert!(
            info.state == EclipseState::Umbra || info.state == EclipseState::Penumbra,
            "GEO satellite behind Earth should be eclipsed: {:?}",
            info.state
        );
        // GEO behind Earth → umbra (Earth angular radius ~8.7° >> Sun ~0.27°)
        assert_eq!(info.state, EclipseState::Umbra);
    }

    #[test]
    fn eclipse_cylindrical_vs_conical_agreement() {
        // Both models should agree on clearly-eclipsed and clearly-sunlit cases
        let sun = [AU, 0.0, 0.0];

        // Clearly sunlit (Sun-side)
        let sat_sunlit = [R_EARTH + 400e3, 0.0, 0.0];
        let cyl = eclipse_cylindrical(sat_sunlit, sun, R_EARTH);
        let con = eclipse_conical(sat_sunlit, sun, R_EARTH, R_SUN);
        assert_eq!(cyl.state, EclipseState::Sunlit);
        assert_eq!(con.state, EclipseState::Sunlit);

        // Clearly in shadow (directly behind Earth)
        let sat_shadow = [-(R_EARTH + 400e3), 0.0, 0.0];
        let cyl = eclipse_cylindrical(sat_shadow, sun, R_EARTH);
        let con = eclipse_conical(sat_shadow, sun, R_EARTH, R_SUN);
        assert_eq!(cyl.state, EclipseState::Umbra);
        assert_eq!(con.state, EclipseState::Umbra);
    }

    #[test]
    fn eclipse_penumbra_region() {
        // Place satellite at the edge of Earth's shadow — should be penumbra
        // in conical model, but binary in cylindrical model.
        let sun = [AU, 0.0, 0.0];
        // Satellite behind Earth, offset by slightly more than R_EARTH
        let offset = R_EARTH * 1.001; // just outside cylindrical shadow
        let sat = [-(R_EARTH + 400e3), offset, 0.0];

        let cyl = eclipse_cylindrical(sat, sun, R_EARTH);
        // Cylindrical: just outside shadow → sunlit
        assert_eq!(cyl.state, EclipseState::Sunlit);

        // Conical may still detect penumbra (or sunlit depending on geometry)
        let con = eclipse_conical(sat, sun, R_EARTH, R_SUN);
        assert!(
            con.state == EclipseState::Sunlit || con.state == EclipseState::Penumbra,
            "edge case should be sunlit or penumbra: {:?}",
            con.state
        );
    }

    // ── Rise/set/transit ─────────────────────────────────────────────

    #[test]
    fn rise_transit_set_sun_at_equinox() {
        // Vernal equinox: Sun at RA=0, Dec=0.
        // Observer at Greenwich (lat=51.5°N, lon=0°).
        // Sun should transit around noon (~12h UT).
        let jd_0h = calendar_to_jd(2024, 3, 20.0).unwrap();
        let lat = 51.5_f64.to_radians();
        let lon = 0.0;
        let ra = 0.0;
        let dec = 0.0;

        let rts = rise_transit_set(jd_0h, ra, dec, lat, lon, STANDARD_REFRACTION).unwrap();

        assert!(rts.rise.is_some(), "Sun should rise at 51.5°N");
        assert!(rts.transit.is_some(), "Sun should transit");
        assert!(rts.set.is_some(), "Sun should set at 51.5°N");

        // Transit should be near 0.5 days (noon)
        let transit_frac = rts.transit.unwrap() - jd_0h;
        assert!(
            (transit_frac - 0.5).abs() < 0.1,
            "transit at {transit_frac:.3} days, expected ~0.5"
        );

        // Rise before transit, set after transit
        assert!(rts.rise.unwrap() < rts.transit.unwrap());
        assert!(rts.set.unwrap() > rts.transit.unwrap());
    }

    #[test]
    fn rise_transit_set_circumpolar() {
        // Star at Dec=+80° seen from 70°N latitude → circumpolar (never sets)
        let jd_0h = calendar_to_jd(2024, 6, 21.0).unwrap();
        let lat = 70.0_f64.to_radians();
        let dec = 80.0_f64.to_radians();
        let ra = 1.0;

        let rts = rise_transit_set(jd_0h, ra, dec, lat, 0.0, 0.0).unwrap();

        // Circumpolar: rise and set are None, transit exists
        assert!(rts.rise.is_none(), "circumpolar body should not rise");
        assert!(rts.set.is_none(), "circumpolar body should not set");
        assert!(rts.transit.is_some(), "circumpolar body should transit");
    }

    #[test]
    fn rise_transit_set_never_rises() {
        // Star at Dec=−80° seen from 70°N latitude → never above horizon
        let jd_0h = calendar_to_jd(2024, 6, 21.0).unwrap();
        let lat = 70.0_f64.to_radians();
        let dec = -80.0_f64.to_radians();
        let ra = 1.0;

        let rts = rise_transit_set(jd_0h, ra, dec, lat, 0.0, 0.0).unwrap();

        assert!(rts.rise.is_none());
        assert!(rts.set.is_none());
        assert!(rts.transit.is_none());
    }

    #[test]
    fn rise_transit_set_equatorial_observer() {
        // Observer at equator: all bodies should rise and set (except pole stars)
        let jd_0h = calendar_to_jd(2024, 1, 15.0).unwrap();
        let lat = 0.0;
        let dec = 30.0_f64.to_radians(); // Dec=30°

        let rts = rise_transit_set(jd_0h, 3.0, dec, lat, 0.0, 0.0).unwrap();
        assert!(rts.rise.is_some());
        assert!(rts.set.is_some());
        assert!(rts.transit.is_some());

        // Day length should be ~12 hours at equator
        let mut day_frac = rts.set.unwrap() - rts.rise.unwrap();
        if day_frac < 0.0 {
            day_frac += 1.0; // set wrapped past midnight
        }
        assert!(
            day_frac > 0.3 && day_frac < 0.7,
            "day fraction={day_frac:.3}, expected ~0.5"
        );
    }

    #[test]
    fn standard_refraction_value() {
        // Verify the constant is approximately -0.5667°
        let deg = STANDARD_REFRACTION.to_degrees();
        assert!(
            (deg + 0.5667).abs() < 0.001,
            "standard refraction = {deg}°, expected -0.5667°"
        );
    }
}