use crate::lighting::solar_position::NEPTUNEAXIALTILTDEG;
use crate::temporal::epoch::J2000EPOCH;

pub const AXIALTILTDEG: f64 = NEPTUNEAXIALTILTDEG;
pub const NEPTUNEORBITALPERIODDAYS: f64 = 60182.0;
pub const VERNALEQUINOXJD: f64 = 2470213.0;

#[derive(Debug, Clone, Copy, PartialEq)]
pub enum Season {
    Spring,
    Summer,
    Autumn,
    Winter,
}

pub struct SeasonalState {
    pub seasonnorth: Season,
    pub seasonsouth: Season,
    pub solardeclinationdeg: f64,
    pub daylengthhours: f64,
}

pub fn seasonat(juliandate: f64, latdeg: f64) -> SeasonalState {
    let dayssinceequinox = juliandate - VERNALEQUINOXJD;
    let yearfraction =
        ((dayssinceequinox % NEPTUNEORBITALPERIODDAYS) / NEPTUNEORBITALPERIODDAYS).rem_euclid(1.0);
    let solardeclination = AXIALTILTDEG * (2.0 * std::f64::consts::PI * yearfraction).sin();
    let northseason = match yearfraction {
        f if f < 0.25 => Season::Spring,
        f if f < 0.5 => Season::Summer,
        f if f < 0.75 => Season::Autumn,
        f if f >= 0.75 => Season::Winter,
        f => {
            debug_assert!(f.is_nan());
            Season::Winter
        }
    };
    let southseason = match northseason {
        Season::Spring => Season::Autumn,
        Season::Summer => Season::Winter,
        Season::Autumn => Season::Spring,
        Season::Winter => Season::Summer,
    };
    let latrad = latdeg.to_radians();
    let decrad = solardeclination.to_radians();
    let coshourangle = (-latrad.tan() * decrad.tan()).clamp(-1.0, 1.0);
    let hourangle = coshourangle.acos();
    let neptunedayhours = 16.11;
    let daylengthhours = neptunedayhours * hourangle / std::f64::consts::PI;
    SeasonalState {
        seasonnorth: northseason,
        seasonsouth: southseason,
        solardeclinationdeg: solardeclination,
        daylengthhours,
    }
}

pub fn subsolarpoint(juliandate: f64) -> (f64, f64) {
    let d = juliandate - J2000EPOCH;
    let meananomaly = (256.228 + 360.0 * d / NEPTUNEORBITALPERIODDAYS) % 360.0;
    let mrad = meananomaly.to_radians();
    let eoc = 0.6541 * mrad.sin() + 0.0033 * (2.0 * mrad).sin();
    let ecllon = (meananomaly + eoc + 180.0 + 44.971) % 360.0;
    let eclrad = ecllon.to_radians();
    let obliquityrad = AXIALTILTDEG.to_radians();
    let declination = (obliquityrad.sin() * eclrad.sin()).asin().to_degrees();
    let neptunesiderealdays = 0.6713;
    let gmst = (280.0 + 360.0 / neptunesiderealdays * d) % 360.0;
    let subsolarlon = -gmst % 360.0;
    (declination, subsolarlon)
}
