earths 0.0.1

High-fidelity Earth simulation engine — orbit, atmosphere, geology, hydrology, biosphere, terrain, lighting, rendering, satellites, and temporal systems with full scientific coupling
Documentation
use crate::temporal::epoch::J2000_EPOCH;
use sciforge::hub::prelude::constants::AU;
pub const EARTH_AXIAL_TILT_DEG: f64 = 23.4393;
pub struct SolarPosition {
    pub azimuth_deg: f64,
    pub elevation_deg: f64,
    pub distance_au: f64,
    pub direction: [f64; 3],
}
impl SolarPosition {
    pub fn compute(julian_date: f64, lat_deg: f64, lon_deg: f64) -> Self {
        let d = julian_date - J2000_EPOCH;
        let mean_anomaly = (357.5291 + 0.98560028 * d) % 360.0;
        let m_rad = mean_anomaly.to_radians();
        let equation_of_center =
            1.9148 * m_rad.sin() + 0.0200 * (2.0 * m_rad).sin() + 0.0003 * (3.0 * m_rad).sin();
        let l_venus = (3.176_146_7 + 1_021.328_554_4 * d / 36525.0) % (2.0 * std::f64::consts::PI);
        let l_jupiter = (0.599_546 + 145.270_285_3 * d / 36525.0) % (2.0 * std::f64::consts::PI);
        let perturbation = 0.00313 * (2.0 * l_venus - m_rad).sin()
            + 0.00178 * (3.0 * m_rad).sin()
            + 0.00077 * (l_jupiter - m_rad).sin();
        let ecliptic_lon =
            (mean_anomaly + equation_of_center + perturbation + 180.0 + 102.9372) % 360.0;
        let ecl_rad = ecliptic_lon.to_radians();
        let obliquity_rad = EARTH_AXIAL_TILT_DEG.to_radians();
        let declination = (obliquity_rad.sin() * ecl_rad.sin()).asin();
        let right_ascension = (obliquity_rad.cos() * ecl_rad.sin()).atan2(ecl_rad.cos());
        let sidereal_time = (280.16 + 360.9856235 * d + lon_deg) % 360.0;
        let hour_angle = sidereal_time.to_radians() - right_ascension;
        let lat_rad = lat_deg.to_radians();
        let altitude = (lat_rad.sin() * declination.sin()
            + lat_rad.cos() * declination.cos() * hour_angle.cos())
        .asin();
        let azimuth = (declination.sin() - altitude.sin() * lat_rad.sin())
            .atan2(altitude.cos() * lat_rad.cos());
        let dir_x = altitude.cos() * azimuth.sin();
        let dir_y = altitude.sin();
        let dir_z = altitude.cos() * azimuth.cos();
        Self {
            azimuth_deg: azimuth.to_degrees(),
            elevation_deg: altitude.to_degrees(),
            distance_au: 1.0 - 0.0167 * m_rad.cos(),
            direction: [dir_x, dir_y, dir_z],
        }
    }
    pub fn distance_m(&self) -> f64 {
        self.distance_au * AU
    }
    pub fn is_above_horizon(&self) -> bool {
        self.elevation_deg > 0.0
    }
}