use crate::lighting::solar_position::EARTHAXIALTILTDEG;
use crate::temporal::epoch::J2000EPOCH;
pub const AXIALTILTDEG: f64 = EARTHAXIALTILTDEG;
pub const TROPICALYEARDAYS: f64 = 365.24219;
pub const VERNALEQUINOXJD: f64 = 2451623.80;
#[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 % TROPICALYEARDAYS) / TROPICALYEARDAYS;
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 daylengthhours = 24.0 * 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 = (357.5291 + 0.98560028 * d) % 360.0;
let mrad = meananomaly.to_radians();
let eoc = 1.9148 * mrad.sin() + 0.02 * (2.0 * mrad).sin();
let ecllon = (meananomaly + eoc + 180.0 + 102.9372) % 360.0;
let eclrad = ecllon.to_radians();
let obliquityrad = AXIALTILTDEG.to_radians();
let declination = (obliquityrad.sin() * eclrad.sin()).asin().to_degrees();
let eotminutes = -7.655 * mrad.sin() + 9.873 * (2.0 * eclrad - mrad).sin();
let solarnoonoffset = -eotminutes / 60.0 * 15.0;
let gmst = (280.46061837 + 360.98564736629 * d) % 360.0;
let subsolarlon = -(gmst + solarnoonoffset) % 360.0;
(declination, subsolarlon)
}