earths 0.0.4

High-fidelity Earth simulation engine — orbit, atmosphere, geology, hydrology, biosphere, terrain, lighting, rendering, satellites, and temporal systems with full scientific coupling
Documentation
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)
}