//! Astronomical math for the almanac. Self-contained, no external API calls.
//!
//! All formulas trace to Jean Meeus, *Astronomical Algorithms* (2nd ed.):
//!   ch. 7  julian day
//!   ch. 22 mean obliquity of the ecliptic
//!   ch. 25 sun apparent position
//!   ch. 28 equation of time
//!   ch. 47 moon position and phase (table 47.A perturbation series, truncated)
//!
//! Conventions:
//!   * Time variable `t` is always Julian centuries since J2000.0
//!     (JD 2451545.0 = 2000-01-01 12:00 TT). One Julian century = 36525 days.
//!     We treat UTC as TT; the offset is ~70 sec, well below render precision.
//!   * Angle <-> time conversions: 1 deg of arc = 4 min of time, 15 deg = 1 hr
//!     (earth rotates 360 deg per 24*60 min).
//!   * Latitude is north-positive, longitude east-positive (Raleigh's lon is
//!     therefore negative).

use chrono::{DateTime, Datelike, Duration, TimeZone, Timelike, Utc};
use chrono_tz::Tz;

/// Latitude of Raleigh, NC (zone 7a target).
const LAT: f64 = 35.78;
/// Longitude of Raleigh, NC. East-positive, so this is negative.
const LON: f64 = -78.64;
/// Reference epoch for every "centuries since" calculation.
const J2000_JD: f64 = 2451545.0;
/// Days in one Julian century (the unit of `t` everywhere below).
const JULIAN_CENTURY: f64 = 36525.0;
/// Mean length of the lunar synodic month (new moon to new moon), in days.
const SYNODIC_MONTH: f64 = 29.53058867;

/// Julian Day for a UTC instant. Meeus formula 7.1, Gregorian calendar only
/// (fine for any modern date). Constants explained inline below.
fn julian_day(dt_utc: DateTime<Utc>) -> f64 {
    let mut y = dt_utc.year();
    let mut m = dt_utc.month() as i32;
    let d = dt_utc.day() as f64
        + (dt_utc.hour() as f64
            + (dt_utc.minute() as f64 + dt_utc.second() as f64 / 60.0) / 60.0)
            / 24.0;
    // Meeus's calendar trick: Jan/Feb are treated as months 13/14 of the
    // previous year so a single polynomial covers the whole shifted year.
    if m <= 2 {
        y -= 1;
        m += 12;
    }
    // `b` is the Gregorian leap-day correction (number of "missing" leap days
    // in the Gregorian rules vs. the pure Julian 365.25 average).
    let a = (y as f64 / 100.0).floor() as i32;
    let b = 2 - a + (a as f64 / 4.0).floor() as i32;
    // 365.25  = avg days/year (Julian)
    // 30.6001 = avg days/month over the shifted Mar->Feb calendar
    // 4716    = year offset that keeps the result non-negative for all dates
    // -1524.5 = anchors JD 0 to noon on -4712 Jan 1 (Julian proleptic), the
    //           astronomical zero point.
    (365.25 * (y as f64 + 4716.0)).floor()
        + (30.6001 * (m as f64 + 1.0)).floor()
        + d
        + b as f64
        - 1524.5
}

/// (age in days through the synodic month, illumination 0..1).
/// Truncated form of Meeus chapter 47 with the table 47.A perturbation series.
fn moon_state(date: DateTime<Tz>) -> (f64, f64) {
    let dt_utc = date.with_timezone(&Utc);
    let t = (julian_day(dt_utc) - J2000_JD) / JULIAN_CENTURY;
    // Meeus 47.2 - 47.5: the four fundamental angles of the lunar problem
    // (all degrees). Each polynomial is "value at J2000.0 + linear drift per
    // Julian century"; higher-order terms are below our render precision.
    //   d  = mean elongation of moon from sun
    //   ms = sun's mean anomaly
    //   mm = moon's mean anomaly
    //   f  = argument of latitude (moon's distance from ascending node)
    let d = (297.8501921_f64 + 445267.1114034 * t).rem_euclid(360.0);
    let ms = (357.5291092_f64 + 35999.0502909 * t).rem_euclid(360.0);
    let mm = (134.9633964_f64 + 477198.8675055 * t).rem_euclid(360.0);
    let f = (93.2720950_f64 + 483202.0175233 * t).rem_euclid(360.0);
    let dr = d.to_radians();
    let msr = ms.to_radians();
    let mmr = mm.to_radians();
    let fr = f.to_radians();
    // Meeus table 47.A: longitude perturbations of the moon (degrees). First
    // 13 terms are above the arcminute level; the rest of the table is too
    // small to matter for percent-illumination output. Coefficients come from
    // the table as-published.
    let dl_moon = 6.288774 * mmr.sin()
        + 1.274027 * (2.0 * dr - mmr).sin()
        + 0.658314 * (2.0 * dr).sin()
        + 0.213618 * (2.0 * mmr).sin()
        - 0.185116 * msr.sin()
        - 0.114332 * (2.0 * fr).sin()
        + 0.058793 * (2.0 * dr - 2.0 * mmr).sin()
        + 0.057066 * (2.0 * dr - msr - mmr).sin()
        + 0.053322 * (2.0 * dr + mmr).sin()
        + 0.045758 * (2.0 * dr - msr).sin()
        - 0.040923 * (msr - mmr).sin()
        - 0.034720 * dr.sin()
        - 0.030383 * (msr + mmr).sin();
    // Sun's longitude correction (the equation of centre, kept to 3 terms).
    let dl_sun = 1.914602 * msr.sin()
        + 0.019993 * (2.0 * msr).sin()
        + 0.000289 * (3.0 * msr).sin();
    // Elongation: how far around the cycle the moon is from the sun (degrees).
    let elong = (d + dl_moon - dl_sun).rem_euclid(360.0);
    // Map elongation linearly onto days through the mean synodic month.
    let age = elong / 360.0 * SYNODIC_MONTH;
    // Half-angle identity: fraction of disk illuminated, exact for a sphere.
    let illum = (1.0 - elong.to_radians().cos()) / 2.0;
    (age, illum)
}

pub fn moon_phase(date: DateTime<Tz>) -> f64 {
    moon_state(date).0
}

pub fn moon_illumination(date: DateTime<Tz>) -> f64 {
    moon_state(date).1
}

/// Name for a given lunar age in days. The four "moment" phases (new, first
/// quarter, full, last quarter) get a thin ~1.85-day window around the exact
/// moment (0.92 days each side); the four crescent/gibbous phases fill the
/// rest. Boundaries are placed so the cycle divides into eight pieces with
/// the moments centred on 0, 7.38 (cycle/4), 14.77 (cycle/2), 22.15 (3*cycle/4).
pub fn moon_name(phase: f64) -> &'static str {
    if phase < 1.85 {
        "new moon"
    } else if phase < 7.38 {
        "waxing crescent"
    } else if phase < 9.23 {
        "first quarter"
    } else if phase < 14.77 {
        "waxing gibbous"
    } else if phase < 16.61 {
        "full moon"
    } else if phase < 22.15 {
        "waning gibbous"
    } else if phase < 23.99 {
        "last quarter"
    } else if phase < 27.68 {
        "waning crescent"
    } else {
        "new moon"
    }
}

// Standard sun-event altitudes, in degrees of the sun's centre below the
// horizon. Sunrise/sunset uses -0.8333 deg = -50' total, which is the sum of
// average atmospheric refraction at the horizon (~34') and the sun's apparent
// semi-diameter (~16'); the upper limb appears tangent to the horizon at this
// geometric altitude. Twilight thresholds are the usual conventions.
const SUN_GEOMETRIC_ALT: f64 = -0.8333;
const CIVIL_DUSK_ALT: f64 = -6.0;
const NAUTICAL_DUSK_ALT: f64 = -12.0;
const ASTRO_DUSK_ALT: f64 = -18.0;

#[derive(Clone, Copy, Debug)]
pub struct SunTimes {
    pub sunrise: f64,
    pub sunset: f64,
    pub civil_dusk: f64,
    pub nautical_dusk: f64,
    pub astronomical_dusk: f64,
    pub daylight_hours: f64,
}

/// Apparent geocentric position of the sun, plus equation of time.
/// Returns (right ascension deg, declination rad, equation of time minutes).
/// Implements Meeus chapter 25 with low-accuracy nutation (chapter 22 / 25.8).
fn sun_apparent(jde: f64) -> (f64, f64, f64) {
    // `t` = Julian centuries since J2000.0, the unit every polynomial below
    // expects. All polynomial coefficients come straight from Meeus.
    let t = (jde - J2000_JD) / JULIAN_CENTURY;
    // 25.2: sun's mean longitude L0 (deg). Linear term ~36000 deg/century =
    // one full cycle/year. Quadratic term is tiny secular drift.
    let l0 = (280.46646 + 36000.76983 * t + 0.0003032 * t * t).rem_euclid(360.0);
    // 25.3: sun's mean anomaly M (deg). Same order of magnitude as L0; differs
    // because mean longitude tracks "where the sun would be on a circle" while
    // anomaly tracks position relative to perihelion.
    let m = (357.52911 + 35999.05029 * t - 0.0001537 * t * t).rem_euclid(360.0);
    let mr = m.to_radians();
    // Meeus 25, "equation of the centre" series. Three-term truncation gives
    // sub-arcsecond accuracy. Result is a correction to L0 to get the sun's
    // true (geometric) longitude.
    let c = (1.914602 - 0.004817 * t - 0.000014 * t * t) * mr.sin()
        + (0.019993 - 0.000101 * t) * (2.0 * mr).sin()
        + 0.000289 * (3.0 * mr).sin();
    let true_long = l0 + c;
    // Omega = longitude of ascending node of the moon's mean orbit (deg).
    // Required by the low-accuracy nutation corrections in 25.8 and 22.
    let omega_r = (125.04 - 1934.136 * t).to_radians();
    // 25.8: apparent longitude (corrects for aberration -0.00569 and the dominant
    // nutation term -0.00478*sin(Omega), both in degrees).
    let lambda = true_long - 0.00569 - 0.00478 * omega_r.sin();
    // 22.2: mean obliquity of the ecliptic (deg). Constant term =
    // 23 deg 26' 21.448" = 23.439291 deg; secular terms ~0.013 deg/century.
    let eps0 = 23.439291 - 0.0130042 * t - 1.64e-7 * t * t + 5.04e-7 * t * t * t;
    // Apparent obliquity = mean + dominant nutation term in obliquity.
    let eps = eps0 + 0.00256 * omega_r.cos();
    let lr = lambda.to_radians();
    let er = eps.to_radians();
    // 25.6 / 25.7: rotate ecliptic (lambda, beta=0) into equatorial coordinates.
    let alpha = (er.cos() * lr.sin())
        .atan2(lr.cos())
        .to_degrees()
        .rem_euclid(360.0);
    let decl = (er.sin() * lr.sin()).asin();
    // 28.1 (simplified): equation of time = L0 - aberration - alpha (deg),
    // wrapped into (-180, 180], then 4 min per degree gives minutes of time.
    // The 0.0057183 deg = 20.4965" is the standard aberration constant for
    // light's annual deflection.
    let mut diff = (l0 - 0.0057183 - alpha) % 360.0;
    if diff > 180.0 {
        diff -= 360.0;
    }
    if diff < -180.0 {
        diff += 360.0;
    }
    (alpha, decl, 4.0 * diff)
}

/// Half-day length, in hours: how long before/after solar noon the sun is at
/// altitude `alt_rad` for declination `decl`. Returns None if the sun never
/// reaches that altitude on this day (polar day or polar night). Standard
/// "sunrise equation": cos(H) = (sin(h) - sin(phi)*sin(delta)) / (cos(phi)*cos(delta)).
/// Convert from degrees to hours by dividing by 15 (earth turns 15 deg/hr).
fn hour_angle_at_alt(decl: f64, alt_rad: f64) -> Option<f64> {
    let lat_r = LAT.to_radians();
    let cos_h = (alt_rad.sin() - lat_r.sin() * decl.sin()) / (lat_r.cos() * decl.cos());
    if !(-1.0..=1.0).contains(&cos_h) {
        return None;
    }
    Some(cos_h.acos().to_degrees() / 15.0)
}

/// Local-tz fractional hour of the rise (`is_rise=true`) or set event when the
/// sun is at altitude `alt_deg`. Returns `None` if the sun never crosses that
/// altitude on the given date (polar day/night). Two-pass algorithm: first
/// solve using the sun's position at noon UTC, then re-solve at the predicted
/// event time so the answer accounts for declination drift across the day.
/// One refinement is enough for sub-second accuracy at temperate latitudes.
fn sun_event_hours(
    jd_noon: f64,
    utc_midnight: DateTime<Utc>,
    tz: Tz,
    alt_deg: f64,
    is_rise: bool,
) -> Option<f64> {
    let alt_r = alt_deg.to_radians();

    // Pass 1: estimate event using sun position at local noon.
    let (_, decl0, eot0) = sun_apparent(jd_noon);
    let h0 = hour_angle_at_alt(decl0, alt_r)?;
    // Apparent solar noon (UTC, fractional hours):
    //   12.0    = noon UTC at the prime meridian, ignoring everything else
    //   -LON/15 = longitude correction (15 deg = 1 hr; west is negative LON,
    //             so this term shifts solar noon later in UTC)
    //   -eot/60 = equation of time, converted from minutes to hours; this is
    //             how far apparent (sundial) noon is from mean (clock) noon
    let solar_noon_0 = 12.0 - LON / 15.0 - eot0 / 60.0;
    let evt0 = if is_rise { solar_noon_0 - h0 } else { solar_noon_0 + h0 };

    // Pass 2: recompute sun position at the predicted event time. `jd_noon`
    // is at 12:00 UTC, so `jd_noon - 0.5` is 00:00 UTC the same date; adding
    // `evt0 / 24` of a day lands at the event itself.
    let (_, decl1, eot1) = sun_apparent(jd_noon - 0.5 + evt0 / 24.0);
    let h1 = hour_angle_at_alt(decl1, alt_r)?;
    let solar_noon_1 = 12.0 - LON / 15.0 - eot1 / 60.0;
    let evt1 = if is_rise { solar_noon_1 - h1 } else { solar_noon_1 + h1 };

    // Convert UTC fractional hours back to a local-tz fractional hour.
    let event_utc = utc_midnight + Duration::nanoseconds((evt1 * 3600.0 * 1e9) as i64);
    let local = event_utc.with_timezone(&tz);
    Some(local.hour() as f64 + local.minute() as f64 / 60.0 + local.second() as f64 / 3600.0)
}

/// Sunrise, sunset, and three twilight ends (civil/nautical/astronomical dusk)
/// for the given local date. Uses Meeus chapter 25 sun position with one
/// refinement pass; matches USNO to a few seconds. Hours are local fractional
/// hours within the calendar date of `local_date`.
pub fn sun_times(local_date: DateTime<Tz>) -> SunTimes {
    let tz = local_date.timezone();
    let local_noon = tz
        .with_ymd_and_hms(
            local_date.year(),
            local_date.month(),
            local_date.day(),
            12,
            0,
            0,
        )
        .unwrap();
    let noon_utc = local_noon.with_timezone(&Utc);
    let utc_midnight = Utc
        .with_ymd_and_hms(noon_utc.year(), noon_utc.month(), noon_utc.day(), 0, 0, 0)
        .unwrap();
    let jd_noon = julian_day(noon_utc);
    let evt = |alt, is_rise| {
        sun_event_hours(jd_noon, utc_midnight, tz, alt, is_rise).unwrap_or(0.0)
    };
    let sunrise = evt(SUN_GEOMETRIC_ALT, true);
    let sunset = evt(SUN_GEOMETRIC_ALT, false);
    SunTimes {
        sunrise,
        sunset,
        civil_dusk: evt(CIVIL_DUSK_ALT, false),
        nautical_dusk: evt(NAUTICAL_DUSK_ALT, false),
        astronomical_dusk: evt(ASTRO_DUSK_ALT, false),
        daylight_hours: sunset - sunrise,
    }
}

pub fn format_hm(hours: f64) -> String {
    let h = hours.trunc() as i64;
    let m = ((hours - h as f64) * 60.0).round() as i64;
    format!("{h}h {m}m")
}

pub fn format_clock(hours: f64) -> String {
    let mut h = hours.trunc() as i64;
    let mut m = ((hours - h as f64) * 60.0).round() as i64;
    if m == 60 {
        h += 1;
        m = 0;
    }
    let suffix = if h >= 12 { "pm" } else { "am" };
    let display = if h > 12 {
        h - 12
    } else if h == 0 {
        12
    } else {
        h
    };
    format!("{display}:{m:02} {suffix}")
}

pub fn sky_data_lines(now: DateTime<Tz>) -> Vec<String> {
    let phase = moon_phase(now);
    let name = moon_name(phase);
    let illum = (moon_illumination(now) * 100.0).round() as i64;
    let st = sun_times(now);
    let yesterday = now - Duration::days(1);
    let gained = (st.daylight_hours - sun_times(yesterday).daylight_hours) * 60.0;
    let sign = if gained > 0.0 { "+" } else { "" };
    vec![
        format!("<strong>{name}</strong>, {illum}% lit"),
        format!(
            "sunrise <strong>{}</strong> \u{00b7} sunset <strong>{}</strong>",
            format_clock(st.sunrise),
            format_clock(st.sunset)
        ),
        format!(
            "<strong>{}</strong> of daylight ({sign}{:.1} minutes from yesterday)",
            format_hm(st.daylight_hours),
            gained
        ),
        format!(
            "civil dusk <strong>{}</strong> \u{00b7} sailor's dark <strong>{}</strong> \u{00b7} true dark <strong>{}</strong>",
            format_clock(st.civil_dusk),
            format_clock(st.nautical_dusk),
            format_clock(st.astronomical_dusk)
        ),
    ]
}

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

    fn fixed_now() -> DateTime<Tz> {
        New_York.with_ymd_and_hms(2026, 5, 6, 19, 0, 0).unwrap()
    }

    #[test]
    fn moon_matches_python() {
        let now = fixed_now();
        let phase = moon_phase(now);
        let illum = moon_illumination(now);
        assert!((phase - 19.474334776875626).abs() < 1e-9, "phase={phase}");
        assert!((illum - 0.7693359060289093).abs() < 1e-9, "illum={illum}");
        assert_eq!(moon_name(phase), "waning gibbous");
    }

    #[test]
    fn sun_matches_usno() {
        // 2026-05-06 USNO/Naval Observatory values for lat 35.78, lon -78.64 EDT:
        //   sunrise 06:17, sunset 20:06, civil dusk 20:34, nautical 21:07, astro 21:43.
        // Tolerance is 30 sec since both USNO and Meeus round display values.
        let now = fixed_now();
        let st = sun_times(now);
        let approx = |actual: f64, h: i32, m: i32, label: &str| {
            let expected = h as f64 + m as f64 / 60.0;
            let diff = (actual - expected).abs();
            assert!(diff < 1.0 / 60.0, "{label}: got {actual}, expected ~{h}:{m:02}");
        };
        approx(st.sunrise, 6, 17, "sunrise");
        approx(st.sunset, 20, 6, "sunset");
        approx(st.civil_dusk, 20, 34, "civil dusk");
        approx(st.nautical_dusk, 21, 7, "nautical dusk");
        approx(st.astronomical_dusk, 21, 43, "astro dusk");
        assert!(
            (st.daylight_hours - 13.81).abs() < 0.02,
            "daylight={}",
            st.daylight_hours
        );
    }
}