marss 0.0.3

Mars 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
67
68
69
70
71
72
73
74
75

pub fn solar_declination_deg(ls_deg: f64) -> f64 {
    let eps = crate::AXIAL_TILT_DEG.to_radians();
    let ls = ls_deg.to_radians();
    (eps.sin() * ls.sin()).asin().to_degrees()
}

pub fn equation_of_time_minutes(ls_deg: f64) -> f64 {
    let e = crate::ECCENTRICITY;
    let ls = ls_deg.to_radians();
    let peri = crate::ARGUMENT_PERIHELION_DEG.to_radians();
    let eot_rad = -2.0 * e * (ls - peri).sin() + 1.25 * e * e * (2.0 * (ls - peri)).sin();
    eot_rad.to_degrees() * 4.0
}

pub struct SolarPosition {
    pub azimuth_deg: f64,
    pub elevation_deg: f64,
    pub distance_au: f64,
    pub direction: [f64; 3],
}

impl SolarPosition {
    pub fn compute(ls_deg: f64, local_time_h: f64, lat_deg: f64, lon_deg: f64) -> Self {
        if !lon_deg.is_finite() {
            return Self {
                azimuth_deg: f64::NAN,
                elevation_deg: f64::NAN,
                distance_au: f64::NAN,
                direction: [f64::NAN; 3],
            };
        }

        let decl = solar_declination_deg(ls_deg);
        let eot = equation_of_time_minutes(ls_deg);

        let hour_angle = (local_time_h - 12.0) * 15.0 + eot / 4.0;
        let ha = hour_angle.to_radians();
        let lat = lat_deg.to_radians();
        let dec = decl.to_radians();

        let sin_elev = lat.sin() * dec.sin() + lat.cos() * dec.cos() * ha.cos();
        let elevation = sin_elev.asin();

        let cos_az = (dec.sin() - lat.sin() * sin_elev) / (lat.cos() * elevation.cos().max(1e-10));
        let mut azimuth = cos_az.clamp(-1.0, 1.0).acos();
        if ha > 0.0 {
            azimuth = 2.0 * std::f64::consts::PI - azimuth;
        }

        let az = azimuth;
        let el = elevation;
        let direction = [-az.sin() * el.cos(), el.sin(), -az.cos() * el.cos()];

        let ls = ls_deg.to_radians();
        let e = crate::ECCENTRICITY;
        let peri = crate::ARGUMENT_PERIHELION_DEG.to_radians();
        let nu = ls - peri;
        let r_au = crate::SEMI_MAJOR_AXIS_AU * (1.0 - e * e) / (1.0 + e * nu.cos());

        Self {
            azimuth_deg: azimuth.to_degrees(),
            elevation_deg: elevation.to_degrees(),
            distance_au: r_au,
            direction,
        }
    }

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

    pub fn distance_m(&self) -> f64 {
        self.distance_au * sciforge::hub::domain::common::constants::AU
    }
}