nsb 0.1.0

// SPDX-License-Identifier: AGPL-3.0-only
// Copyright (C) 2026 Vallés Puig, Ramon
use super::domain::AirglowNightPhase;
use crate::units::angular::Degrees;
use chrono::{Datelike, Timelike};
use siderust::bodies::Sun as SunBody;
use siderust::coordinates::centers::Geodetic;
use siderust::coordinates::frames::ECEF;
use siderust::event::altitude::{AltitudeEventsExt, SearchOpts};
use siderust::qtty::Days;
use siderust::time::{Interval as TimePeriod, ModifiedJulianDate, TT};
use tempoch::{Time, MJD, UTC};

pub(crate) const ASTRONOMICAL_TWILIGHT: Degrees = Degrees::new(-18.0);
const INITIAL_NIGHT_SEARCH_RADIUS: Days = Days::new(2.0);
const MAX_NIGHT_SEARCH_RADIUS: Days = Days::new(200.0);
const NIGHT_SEARCH_EXPANSION_FACTOR: f64 = 4.0;

#[derive(Debug, Clone, Copy)]
pub(crate) struct AstronomicalNightPeriod {
    pub(crate) period: TimePeriod<ModifiedJulianDate>,
    pub(crate) phase_bounded: bool,
}

#[derive(Debug, Clone, Copy)]
pub(crate) struct AirglowPhasePeriod {
    pub(crate) period: TimePeriod<ModifiedJulianDate>,
    #[cfg_attr(not(debug_assertions), allow(dead_code))]
    pub(crate) phase: AirglowNightPhase,
}

/// PALACE month and one-hour local-mean-solar-time bin for the planning proxy.
///
/// PALACE's climatology was fitted at Cerro Paranal. For arbitrary-location
/// planning, NSB transfers that climatology by local mean solar clock: the
/// observer's local solar month/hour selects the equivalent PALACE bin. This is
/// an explicit spatial-transfer assumption, not evidence that PALACE was fitted
/// at the observer longitude. Bins 1..=12 cover 18:00..06:00; astronomical-night
/// samples outside that interval use the nearest endpoint bin.
pub(crate) fn palace_climatology_coordinates(
    time: Time<UTC>,
    location: Geodetic<ECEF>,
) -> Option<(u32, usize)> {
    let dt = observer_local_mean_solar_datetime(time, location)?;
    let hour = dt.hour();
    let bin = match hour {
        18..=23 => (hour - 17) as usize,
        0..=5 => (hour + 7) as usize,
        6..=11 => 12,
        12..=17 => 1,
        _ => unreachable!("chrono hour is always 0..=23"),
    };
    Some((dt.month(), bin))
}

/// Return every local-mean-solar hour boundary strictly inside a TT window.
///
/// PALACE selects a discrete climatology row from calendar month and local
/// mean solar hour. Threshold scans must therefore treat these boundaries as
/// potential discontinuities independently of the legacy astronomical-night
/// third boundaries.
pub(crate) fn palace_climatology_boundaries_for_window(
    window: TimePeriod<ModifiedJulianDate>,
    location: Geodetic<ECEF>,
) -> Vec<ModifiedJulianDate> {
    if window.start >= window.end {
        return Vec::new();
    }

    let start_utc = Time::<TT>::from(window.start).to::<UTC>();
    let end_utc = Time::<TT>::from(window.end).to::<UTC>();
    let Some(start_local) = observer_local_mean_solar_datetime(start_utc, location) else {
        return Vec::new();
    };
    let Some(end_local) = observer_local_mean_solar_datetime(end_utc, location) else {
        return Vec::new();
    };
    let Some(mut boundary_local) = start_local
        .with_minute(0)
        .and_then(|dt| dt.with_second(0))
        .and_then(|dt| dt.with_nanosecond(0))
    else {
        return Vec::new();
    };
    if boundary_local <= start_local {
        boundary_local += chrono::Duration::hours(1);
    }

    let offset_seconds = local_mean_solar_offset_seconds(location);
    let mut boundaries = Vec::new();
    while boundary_local < end_local {
        let boundary_utc =
            Time::<UTC>::from_chrono(boundary_local - chrono::Duration::seconds(offset_seconds));
        boundaries.push(utc_time_to_tt_mjd(boundary_utc));
        boundary_local += chrono::Duration::hours(1);
    }
    boundaries
}

/// Site-aware Airglow night phase based on astronomical-night thirds.
///
/// Divide the full astronomical-night interval (`alt_sun < -18°`) into thirds
/// for the retained phase identity. PALACE evaluation itself uses local mean
/// solar hour bins, and threshold searches split those discontinuities separately.
///
/// The search expands adaptively so high-latitude winter nights are not
/// mistaken for missing Airglow merely because the first local window is clipped.
/// If the final expanded search still sees continuous night without both bounding
/// crossings, [`AirglowNightPhase::FullNight`] is used explicitly.
pub(crate) fn night_phase(time: Time<UTC>, location: Geodetic<ECEF>) -> Option<AirglowNightPhase> {
    let time_tt = utc_time_to_tt_mjd(time);
    let night = astronomical_night_containing(time_tt, location)?;
    night_phase_from_night(time_tt, &night)
}

pub(crate) fn is_astronomical_twilight(threshold: Degrees) -> bool {
    (threshold - ASTRONOMICAL_TWILIGHT).abs() <= Degrees::new(f64::EPSILON)
}

pub(crate) fn astronomical_nights_for_window(
    window: TimePeriod<ModifiedJulianDate>,
    location: Geodetic<ECEF>,
) -> Vec<AstronomicalNightPeriod> {
    if window.start >= window.end {
        return Vec::new();
    }

    let mut radius = INITIAL_NIGHT_SEARCH_RADIUS;
    loop {
        let search_window = expand_window(window, radius);
        let nights: Vec<_> =
            sun_below_threshold_periods(search_window, location, ASTRONOMICAL_TWILIGHT)
                .into_iter()
                .filter(|night| night.end > window.start && night.start < window.end)
                .map(|night| AstronomicalNightPeriod {
                    phase_bounded: night.start > search_window.start
                        && night.end < search_window.end,
                    period: night,
                })
                .collect();

        if nights.iter().all(|night| night.phase_bounded) || radius >= MAX_NIGHT_SEARCH_RADIUS {
            return nights;
        }
        radius = (radius * NIGHT_SEARCH_EXPANSION_FACTOR).min(MAX_NIGHT_SEARCH_RADIUS);
    }
}

pub(crate) fn clipped_night_periods(
    nights: &[AstronomicalNightPeriod],
    window: TimePeriod<ModifiedJulianDate>,
) -> Vec<TimePeriod<ModifiedJulianDate>> {
    nights
        .iter()
        .filter_map(|night| intersect_period(night.period, window))
        .collect()
}

pub(crate) fn airglow_phase_periods_for_window(
    nights: &[AstronomicalNightPeriod],
    window: TimePeriod<ModifiedJulianDate>,
) -> Vec<AirglowPhasePeriod> {
    let mut periods = Vec::new();
    for night in nights {
        if !night.phase_bounded {
            if let Some(period) = intersect_period(night.period, window) {
                periods.push(AirglowPhasePeriod {
                    period,
                    phase: AirglowNightPhase::FullNight,
                });
            }
            continue;
        }

        let duration = night.period.length();
        if !duration.is_finite() || duration <= Days::new(0.0) {
            continue;
        }

        let first = night.period.start + duration / 3.0;
        let second = night.period.start + duration * 2.0 / 3.0;
        for (period, phase) in [
            (
                TimePeriod::new(night.period.start, first),
                AirglowNightPhase::FirstThird,
            ),
            (
                TimePeriod::new(first, second),
                AirglowNightPhase::MiddleThird,
            ),
            (
                TimePeriod::new(second, night.period.end),
                AirglowNightPhase::LastThird,
            ),
        ] {
            if let Some(period) = intersect_period(period, window) {
                periods.push(AirglowPhasePeriod { period, phase });
            }
        }
    }
    periods
}

pub(crate) fn night_phase_from_nights(
    time_tt: ModifiedJulianDate,
    nights: &[AstronomicalNightPeriod],
) -> Option<AirglowNightPhase> {
    nights
        .iter()
        .find(|night| night.period.start < time_tt && time_tt < night.period.end)
        .and_then(|night| night_phase_from_night(time_tt, night))
}

#[cfg_attr(not(debug_assertions), allow(dead_code))]
pub(crate) fn night_phase_from_phase_periods(
    time_tt: ModifiedJulianDate,
    phases: &[AirglowPhasePeriod],
) -> Option<AirglowNightPhase> {
    phases
        .iter()
        .find(|phase| phase.period.start < time_tt && time_tt < phase.period.end)
        .map(|phase| phase.phase)
}

fn night_phase_from_night(
    time_tt: ModifiedJulianDate,
    night: &AstronomicalNightPeriod,
) -> Option<AirglowNightPhase> {
    if !(night.period.start < time_tt && time_tt < night.period.end) {
        return None;
    }

    if !night.phase_bounded {
        return Some(AirglowNightPhase::FullNight);
    }

    let duration = night.period.length();
    if !duration.is_finite() || duration <= Days::new(0.0) {
        return None;
    }

    let first = night.period.start + duration / 3.0;
    let second = night.period.start + duration * 2.0 / 3.0;

    Some(if time_tt < first {
        AirglowNightPhase::FirstThird
    } else if time_tt < second {
        AirglowNightPhase::MiddleThird
    } else {
        AirglowNightPhase::LastThird
    })
}

fn astronomical_night_containing(
    time_tt: ModifiedJulianDate,
    location: Geodetic<ECEF>,
) -> Option<AstronomicalNightPeriod> {
    let mut radius = INITIAL_NIGHT_SEARCH_RADIUS;

    loop {
        let search_window = TimePeriod::new(time_tt - radius, time_tt + radius);

        let night = sun_below_threshold_periods(search_window, location, ASTRONOMICAL_TWILIGHT)
            .into_iter()
            .find(|night| night.start < time_tt && time_tt < night.end)?;

        if night.start > search_window.start && night.end < search_window.end {
            return Some(AstronomicalNightPeriod {
                period: night,
                phase_bounded: true,
            });
        }

        if radius >= MAX_NIGHT_SEARCH_RADIUS {
            return Some(AstronomicalNightPeriod {
                period: night,
                phase_bounded: false,
            });
        }

        radius = (radius * NIGHT_SEARCH_EXPANSION_FACTOR).min(MAX_NIGHT_SEARCH_RADIUS);
    }
}

/// Delegate solar-event discovery, exact validation, and fallback to Siderust.
pub(crate) fn sun_below_threshold_periods(
    window: TimePeriod<ModifiedJulianDate>,
    location: Geodetic<ECEF>,
    threshold: Degrees,
) -> Vec<TimePeriod<ModifiedJulianDate>> {
    SunBody.below_threshold(&location, window, threshold, SearchOpts::default())
}

fn expand_window(
    window: TimePeriod<ModifiedJulianDate>,
    radius: Days,
) -> TimePeriod<ModifiedJulianDate> {
    TimePeriod::new(window.start - radius, window.end + radius)
}

fn intersect_period(
    lhs: TimePeriod<ModifiedJulianDate>,
    rhs: TimePeriod<ModifiedJulianDate>,
) -> Option<TimePeriod<ModifiedJulianDate>> {
    let start = lhs.start.max(rhs.start);
    let end = lhs.end.min(rhs.end);
    (start < end).then(|| TimePeriod::new(start, end))
}

fn utc_time_to_tt_mjd(time: Time<UTC>) -> ModifiedJulianDate {
    ModifiedJulianDate::from(time.to::<TT>().to::<MJD>())
}

fn local_mean_solar_offset_seconds(location: Geodetic<ECEF>) -> i64 {
    (location.lon.value() / 15.0 * 3600.0).round() as i64
}

fn observer_local_mean_solar_datetime(
    time: Time<UTC>,
    location: Geodetic<ECEF>,
) -> Option<chrono::DateTime<chrono::Utc>> {
    let dt = time.to_chrono()?;
    Some(dt + chrono::Duration::seconds(local_mean_solar_offset_seconds(location)))
}

#[cfg(test)]
pub(crate) mod test_support {
    use super::*;

    pub(crate) fn night_phase(
        time: Time<UTC>,
        location: Geodetic<ECEF>,
    ) -> Option<AirglowNightPhase> {
        super::night_phase(time, location)
    }

    pub(crate) fn astronomical_night(
        time: Time<UTC>,
        location: Geodetic<ECEF>,
    ) -> Option<TimePeriod<ModifiedJulianDate>> {
        astronomical_night_containing(utc_time_to_tt_mjd(time), location).map(|night| night.period)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::units::angular::Degrees;
    use chrono::{TimeZone, Utc};
    use siderust::coordinates::centers::Geodetic;
    use siderust::coordinates::frames::ECEF;
    use siderust::qtty::Meters;

    fn utc(year: i32, month: u32, day: u32) -> Time<UTC> {
        Time::<UTC>::from_chrono(
            Utc.with_ymd_and_hms(year, month, day, 12, 0, 0)
                .single()
                .unwrap(),
        )
    }

    fn equator() -> Geodetic<ECEF> {
        Geodetic::new_raw(Degrees::new(0.0), Degrees::new(0.0), Meters::new(0.0))
    }

    fn cta_south() -> Geodetic<ECEF> {
        Geodetic::new_raw(
            Degrees::new(-70.3147),
            Degrees::new(-24.6834),
            Meters::new(2_147.0),
        )
    }

    fn period(start: (i32, u32, u32), days: i64) -> TimePeriod<ModifiedJulianDate> {
        let start = utc(start.0, start.1, start.2);
        let end =
            Time::<UTC>::from_chrono(start.to_chrono().unwrap() + chrono::Duration::days(days));
        TimePeriod::new(utc_time_to_tt_mjd(start), utc_time_to_tt_mjd(end))
    }

    fn assert_solar_periods_match_reference(
        location: Geodetic<ECEF>,
        window: TimePeriod<ModifiedJulianDate>,
        threshold: Degrees,
    ) {
        let reference =
            SunBody.below_threshold(&location, window, threshold, SearchOpts::default());
        let actual = sun_below_threshold_periods(window, location, threshold);
        assert_eq!(
            actual.len(),
            reference.len(),
            "actual={actual:?} reference={reference:?}"
        );
        for (actual, reference) in actual.iter().zip(reference) {
            let start_error_seconds =
                (actual.start.raw().value() - reference.start.raw().value()).abs() * 86_400.0;
            let end_error_seconds =
                (actual.end.raw().value() - reference.end.raw().value()).abs() * 86_400.0;
            assert!(
                start_error_seconds <= 1.0,
                "solar start error {start_error_seconds}s"
            );
            assert!(
                end_error_seconds <= 1.0,
                "solar end error {end_error_seconds}s"
            );
        }
    }

    #[test]
    fn palace_planning_transfer_uses_observer_local_mean_solar_clock() {
        let time = Time::<UTC>::from_chrono(
            Utc.with_ymd_and_hms(2026, 9, 15, 22, 0, 0)
                .single()
                .unwrap(),
        );
        let greenwich = Geodetic::new_raw(Degrees::new(0.0), Degrees::new(0.0), Meters::new(0.0));
        let east_30 = Geodetic::new_raw(Degrees::new(30.0), Degrees::new(0.0), Meters::new(0.0));

        assert_eq!(
            palace_climatology_coordinates(time, greenwich),
            Some((9, 5))
        );
        assert_eq!(palace_climatology_coordinates(time, east_30), Some((9, 7)));
    }

    #[test]
    fn bounded_night_phase_uses_current_third_boundaries() {
        let night = AstronomicalNightPeriod {
            period: TimePeriod::new(ModifiedJulianDate::new(0.0), ModifiedJulianDate::new(3.0)),
            phase_bounded: true,
        };

        assert_eq!(
            night_phase_from_nights(ModifiedJulianDate::new(0.0), &[night]),
            None
        );

        assert_eq!(
            night_phase_from_night(ModifiedJulianDate::new(0.5), &night),
            Some(AirglowNightPhase::FirstThird)
        );
        assert_eq!(
            night_phase_from_night(ModifiedJulianDate::new(1.0), &night),
            Some(AirglowNightPhase::MiddleThird)
        );
        assert_eq!(
            night_phase_from_night(ModifiedJulianDate::new(1.5), &night),
            Some(AirglowNightPhase::MiddleThird)
        );
        assert_eq!(
            night_phase_from_night(ModifiedJulianDate::new(2.0), &night),
            Some(AirglowNightPhase::LastThird)
        );
        assert_eq!(
            night_phase_from_night(ModifiedJulianDate::new(2.5), &night),
            Some(AirglowNightPhase::LastThird)
        );
        assert_eq!(
            night_phase_from_nights(ModifiedJulianDate::new(3.0), &[night]),
            None
        );

        for time in [0.5, 1.0, 1.5, 2.0, 2.5] {
            let time = ModifiedJulianDate::new(time);
            assert_eq!(
                night_phase_from_nights(time, &[night]),
                night_phase_from_night(time, &night),
                "the threshold hot path must preserve canonical phase semantics"
            );
        }
    }

    #[test]
    fn precomputed_strict_intervals_do_not_define_third_boundary_semantics() {
        let night = AstronomicalNightPeriod {
            period: TimePeriod::new(ModifiedJulianDate::new(0.0), ModifiedJulianDate::new(3.0)),
            phase_bounded: true,
        };
        let phases = airglow_phase_periods_for_window(&[night], night.period);

        for (time, canonical) in [
            (1.0, AirglowNightPhase::MiddleThird),
            (2.0, AirglowNightPhase::LastThird),
        ] {
            let time = ModifiedJulianDate::new(time);
            assert_eq!(night_phase_from_nights(time, &[night]), Some(canonical));
            assert_eq!(night_phase_from_phase_periods(time, &phases), None);
        }
    }

    #[test]
    fn unbounded_astronomical_night_uses_explicit_full_night_phase() {
        let night = AstronomicalNightPeriod {
            period: TimePeriod::new(ModifiedJulianDate::new(0.0), ModifiedJulianDate::new(3.0)),
            phase_bounded: false,
        };
        assert_eq!(
            night_phase_from_night(ModifiedJulianDate::new(1.5), &night),
            Some(AirglowNightPhase::FullNight)
        );
    }

    #[test]
    fn delegated_siderust_solar_events_cover_sites_thresholds_and_year_boundary() {
        for (location, window, threshold) in [
            (cta_south(), period((2026, 1, 1), 31), Degrees::new(-18.0)),
            (cta_south(), period((2026, 6, 1), 31), Degrees::new(0.0)),
            (equator(), period((2026, 3, 15), 20), Degrees::new(-12.0)),
            (equator(), period((2026, 12, 30), 5), Degrees::new(0.0)),
            (
                Geodetic::new_raw(
                    Degrees::new(18.9553),
                    Degrees::new(69.6492),
                    Meters::new(0.0),
                ),
                period((2026, 5, 10), 40),
                Degrees::new(-18.0),
            ),
        ] {
            assert_solar_periods_match_reference(location, window, threshold);
        }
    }
}