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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113

use sciforge::hub::domain::common::constants::{G, R_GAS};
use sciforge::hub::domain::meteorology::radiation::solar_constant;

pub struct AtmosphericLayer {
    pub name: &'static str,
    pub base_altitude_m: f64,
    pub top_altitude_m: f64,
    pub base_temperature_k: f64,
    pub lapse_rate_k_per_m: f64,
}

impl AtmosphericLayer {
    pub fn temperature_at(&self, altitude_m: f64) -> f64 {
        let h = (altitude_m - self.base_altitude_m)
            .max(0.0)
            .min(self.top_altitude_m - self.base_altitude_m);
        self.base_temperature_k - self.lapse_rate_k_per_m * h
    }

    pub fn pressure_at(&self, altitude_m: f64) -> f64 {
        let t = self.temperature_at(altitude_m);
        let h = altitude_m - self.base_altitude_m;
        let g0 = G * crate::MARS_MASS / (crate::MARS_RADIUS * crate::MARS_RADIUS);
        let m = crate::ATMOSPHERE_MOLAR_MASS;
        if self.lapse_rate_k_per_m.abs() < 1e-10 {
            let p_base = barometric_pressure(self.base_altitude_m);
            p_base * (-g0 * m * h / (R_GAS * self.base_temperature_k)).exp()
        } else {
            let p_base = barometric_pressure(self.base_altitude_m);
            p_base * (t / self.base_temperature_k).powf(g0 * m / (R_GAS * self.lapse_rate_k_per_m))
        }
    }

    pub fn density_at(&self, altitude_m: f64) -> f64 {
        let p = self.pressure_at(altitude_m);
        let t = self.temperature_at(altitude_m);
        p * crate::ATMOSPHERE_MOLAR_MASS / (R_GAS * t)
    }
}

pub fn troposphere() -> AtmosphericLayer {
    AtmosphericLayer {
        name: "Troposphere",
        base_altitude_m: 0.0,
        top_altitude_m: 40_000.0,
        base_temperature_k: crate::MEAN_TEMPERATURE_K,
        lapse_rate_k_per_m: crate::LAPSE_RATE,
    }
}

pub fn mesosphere() -> AtmosphericLayer {
    AtmosphericLayer {
        name: "Mesosphere",
        base_altitude_m: 40_000.0,
        top_altitude_m: 100_000.0,
        base_temperature_k: 110.0,
        lapse_rate_k_per_m: 0.0,
    }
}

pub fn thermosphere() -> AtmosphericLayer {
    AtmosphericLayer {
        name: "Thermosphere",
        base_altitude_m: 100_000.0,
        top_altitude_m: 230_000.0,
        base_temperature_k: 110.0,
        lapse_rate_k_per_m: -0.0015,
    }
}

pub fn barometric_pressure(altitude_m: f64) -> f64 {
    let g0 = G * crate::MARS_MASS / (crate::MARS_RADIUS * crate::MARS_RADIUS);
    let m = crate::ATMOSPHERE_MOLAR_MASS;
    let t0 = crate::MEAN_TEMPERATURE_K;
    let p0 = crate::SURFACE_PRESSURE_PA;
    let lapse = crate::LAPSE_RATE;

    if altitude_m <= 40_000.0 {
        let t = t0 - lapse * altitude_m;
        p0 * (t / t0).powf(g0 * m / (R_GAS * lapse))
    } else {
        let p_tropo = barometric_pressure(40_000.0);
        let t_tropo = 110.0;
        p_tropo * (-g0 * m * (altitude_m - 40_000.0) / (R_GAS * t_tropo)).exp()
    }
}

pub fn scale_height() -> f64 {
    let g0 = G * crate::MARS_MASS / (crate::MARS_RADIUS * crate::MARS_RADIUS);
    R_GAS * crate::MEAN_TEMPERATURE_K / (crate::ATMOSPHERE_MOLAR_MASS * g0)
}

pub fn mean_solar_irradiance() -> f64 {
    let d_au = crate::SEMI_MAJOR_AXIS_AU;
    solar_constant() / (d_au * d_au)
}

pub fn effective_temperature_k() -> f64 {
    let s = mean_solar_irradiance();
    let albedo = crate::BOND_ALBEDO;
    ((s * (1.0 - albedo)) / (4.0 * sciforge::hub::domain::common::constants::SIGMA_SB)).powf(0.25)
}

pub fn co2_condenses_at(temperature_k: f64) -> bool {
    temperature_k <= crate::CO2_FROST_POINT_K
}

pub fn seasonal_dust_optical_depth(areocentric_longitude_deg: f64) -> f64 {
    let ls = areocentric_longitude_deg;
    let base = crate::MEAN_DUST_OPTICAL_DEPTH;
    let peak_factor = (-(ls - 250.0).powi(2) / 2000.0).exp();
    base + 2.0 * base * peak_factor
}