mercurys 0.0.3

Mercury celestial simulation crate for the MilkyWay SolarSystem workspace
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

use crate::AXIAL_TILT_RAD;
use crate::ECCENTRICITY;
use crate::SEMI_MAJOR_AXIS_M;
use sciforge::hub::domain::meteorology::radiation::solar_zenith_angle;

pub struct SolarPosition {
    pub elevation_rad: f64,
    pub azimuth_rad: f64,
    pub distance_au: f64,
}

const AU: f64 = 1.496e11;

impl SolarPosition {
    pub fn compute(
        lat_deg: f64,
        lon_deg: f64,
        true_anomaly_rad: f64,
        subsolar_lon_deg: f64,
    ) -> Self {
        let r = SEMI_MAJOR_AXIS_M * (1.0 - ECCENTRICITY * ECCENTRICITY)
            / (1.0 + ECCENTRICITY * true_anomaly_rad.cos());
        let distance_au = r / AU;

        let solar_declination = AXIAL_TILT_RAD.sin() * true_anomaly_rad.sin();
        let dec_rad = solar_declination.asin();

        let hour_angle_deg = lon_deg - subsolar_lon_deg;
        let hour_angle_rad = hour_angle_deg.to_radians();

        let lat_rad = lat_deg.to_radians();
        let sin_elev =
            lat_rad.sin() * dec_rad.sin() + lat_rad.cos() * dec_rad.cos() * hour_angle_rad.cos();
        let elevation_rad = sin_elev.asin();

        let cos_az = (dec_rad.sin() - lat_rad.sin() * sin_elev)
            / (lat_rad.cos() * elevation_rad.cos() + 1e-30);
        let mut azimuth_rad = cos_az.clamp(-1.0, 1.0).acos();
        if hour_angle_rad.sin() > 0.0 {
            azimuth_rad = 2.0 * std::f64::consts::PI - azimuth_rad;
        }

        Self {
            elevation_rad,
            azimuth_rad,
            distance_au,
        }
    }

    pub fn elevation_deg(&self) -> f64 {
        self.elevation_rad.to_degrees()
    }

    pub fn is_above_horizon(&self) -> bool {
        self.elevation_rad > 0.0
    }

    pub fn irradiance(&self) -> f64 {
        let base = sciforge::hub::domain::meteorology::radiation::solar_constant();
        base / (self.distance_au * self.distance_au)
    }
}

pub fn zenith_angle(lat_deg: f64, dec_deg: f64, hour_angle_deg: f64) -> f64 {
    solar_zenith_angle(lat_deg, dec_deg, hour_angle_deg)
}