earths 0.0.3

High-fidelity Earth simulation engine — orbit, atmosphere, geology, hydrology, biosphere, terrain, lighting, rendering, satellites, and temporal systems with full scientific coupling
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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

use sciforge::hub::domain::common::constants::EARTH_RADIUS;
use sciforge::hub::prelude::constants::N_A;
use sciforge::hub::prelude::constants::elements::atomic_mass;
fn rayleighbetargb() -> [f64; 3] {
    let mn2 = 2.0 * atomic_mass(7) * 1e-3;
    let mo2 = 2.0 * atomic_mass(8) * 1e-3;
    let mair = 0.78084 * mn2 + 0.20946 * mo2 + 0.00934 * atomic_mass(18) * 1e-3;
    let ndensity = N_A * *crate::SEALEVELAIRDENSITY / mair;
    let nminus1 = 2.9e-4;
    let nsqminus1sq = (2.0 * nminus1) * (2.0 * nminus1);
    let coeff = 8.0 * std::f64::consts::PI.powi(3) * nsqminus1sq / (3.0 * ndensity);
    let lambdas = [680e-9f64, 550e-9, 440e-9];
    [
        coeff / lambdas[0].powi(4),
        coeff / lambdas[1].powi(4),
        coeff / lambdas[2].powi(4),
    ]
}
pub struct AtmosphereParams {
    pub planetradiusm: f64,
    pub atmosphereheightm: f64,
    pub rayleighcoefficients: [f64; 3],
    pub rayleighscaleheightm: f64,
    pub miecoefficient: f64,
    pub miescaleheightm: f64,
    pub mieasymmetry: f64,
    pub sunintensity: f64,
    pub numsamples: u32,
    pub numlightsamples: u32,
}
impl Default for AtmosphereParams {
    fn default() -> Self {
        let beta = rayleighbetargb();
        Self {
            planetradiusm: EARTH_RADIUS,
            atmosphereheightm: 100000.0,
            rayleighcoefficients: beta,
            rayleighscaleheightm: 8500.0,
            miecoefficient: 21e-6,
            miescaleheightm: 1200.0,
            mieasymmetry: 0.76,
            sunintensity: 22.0,
            numsamples: 16,
            numlightsamples: 8,
        }
    }
}
impl AtmosphereParams {
    pub fn rayleighdensity(&self, altitudem: f64) -> f64 {
        (-altitudem / self.rayleighscaleheightm).exp()
    }
    pub fn miedensity(&self, altitudem: f64) -> f64 {
        (-altitudem / self.miescaleheightm).exp()
    }
    pub fn scatterrayleigh(&self, costheta: f64) -> f64 {
        3.0 / (16.0 * std::f64::consts::PI) * (1.0 + costheta * costheta)
    }
    pub fn scattermie(&self, costheta: f64) -> f64 {
        let g = self.mieasymmetry;
        let g2 = g * g;
        let num = 3.0 * (1.0 - g2) * (1.0 + costheta * costheta);
        let den =
            8.0 * std::f64::consts::PI * (2.0 + g2) * (1.0 + g2 - 2.0 * g * costheta).powf(1.5);
        num / den
    }
    pub fn skycolor(&self, sundir: [f64; 3], viewdir: [f64; 3], cameraaltitudem: f64) -> [f64; 3] {
        let costheta = sundir[0] * viewdir[0] + sundir[1] * viewdir[1] + sundir[2] * viewdir[2];
        let rayleighphase = self.scatterrayleigh(costheta);
        let miephase = self.scattermie(costheta);
        const LAMBDAR: f64 = 680.0;
        const LAMBDAG: f64 = 550.0;
        const LAMBDAB: f64 = 440.0;
        let inv4r = (LAMBDAG / LAMBDAR).powi(4);
        let inv4g = 1.0;
        let inv4b = (LAMBDAG / LAMBDAB).powi(4);
        let nsteps = self.numsamples as usize;
        let atmotop = self.planetradiusm + self.atmosphereheightm;
        let raylen = atmotop - (self.planetradiusm + cameraaltitudem);
        let steplen = raylen / nsteps as f64;
        let mut color = [0.0f64; 3];
        let mut opticalr = 0.0;
        let mut opticalg = 0.0;
        let mut opticalb = 0.0;
        for i in 0..nsteps {
            let h = cameraaltitudem + steplen * (i as f64 + 0.5);
            let rhoray = self.rayleighdensity(h);
            let rhomie = self.miedensity(h);
            opticalr += (self.rayleighcoefficients[0] * inv4r * rhoray
                + self.miecoefficient * rhomie)
                * steplen;
            opticalg += (self.rayleighcoefficients[1] * inv4g * rhoray
                + self.miecoefficient * rhomie)
                * steplen;
            opticalb += (self.rayleighcoefficients[2] * inv4b * rhoray
                + self.miecoefficient * rhomie)
                * steplen;
            let nlight = self.numlightsamples as usize;
            let lightstep = (atmotop - (self.planetradiusm + h)) / nlight as f64;
            let mut sunr = 0.0;
            let mut sung = 0.0;
            let mut sunb = 0.0;
            for j in 0..nlight {
                let lh = h + lightstep * (j as f64 + 0.5);
                let lr = self.rayleighdensity(lh);
                let lm = self.miedensity(lh);
                sunr += (self.rayleighcoefficients[0] * inv4r * lr + self.miecoefficient * lm)
                    * lightstep;
                sung += (self.rayleighcoefficients[1] * inv4g * lr + self.miecoefficient * lm)
                    * lightstep;
                sunb += (self.rayleighcoefficients[2] * inv4b * lr + self.miecoefficient * lm)
                    * lightstep;
            }
            let attenr = (-sunr - opticalr).exp();
            let atteng = (-sung - opticalg).exp();
            let attenb = (-sunb - opticalb).exp();
            color[0] += (self.rayleighcoefficients[0] * inv4r * rhoray * rayleighphase
                + self.miecoefficient * rhomie * miephase)
                * attenr
                * steplen;
            color[1] += (self.rayleighcoefficients[1] * inv4g * rhoray * rayleighphase
                + self.miecoefficient * rhomie * miephase)
                * atteng
                * steplen;
            color[2] += (self.rayleighcoefficients[2] * inv4b * rhoray * rayleighphase
                + self.miecoefficient * rhomie * miephase)
                * attenb
                * steplen;
        }
        [
            self.sunintensity * color[0],
            self.sunintensity * color[1],
            self.sunintensity * color[2],
        ]
    }
}