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
136
137
138
139

use sciforge::hub::prelude::constants::{EARTH_MASS, EARTH_RADIUS, G};
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum OrbitType {
    LEO,
    MEO,
    GEO,
    HEO,
    Custom,
}
pub struct ArtificialSatellite {
    pub name: String,
    pub masskg: f64,
    pub altitudem: f64,
    pub semimajoraxism: f64,
    pub eccentricity: f64,
    pub inclinationrad: f64,
    pub raanrad: f64,
    pub argperigeerad: f64,
    pub orbitalanglerad: f64,
    pub orbittype: OrbitType,
}
impl ArtificialSatellite {
    pub fn new(
        name: &str,
        masskg: f64,
        altitudem: f64,
        eccentricity: f64,
        inclinationrad: f64,
    ) -> Self {
        let semimajoraxism = EARTH_RADIUS + altitudem;
        let orbittype = match altitudem {
            a if a < 2000000.0 => OrbitType::LEO,
            a if a < 20200000.0 => OrbitType::MEO,
            a if (35786000.0 - a).abs() < 100000.0 => OrbitType::GEO,
            a if a > 35786000.0 => OrbitType::HEO,
            a if a.is_finite() => OrbitType::Custom,
            a => {
                debug_assert!(!a.is_finite());
                OrbitType::Custom
            }
        };
        Self {
            name: name.to_string(),
            masskg,
            altitudem,
            semimajoraxism,
            eccentricity,
            inclinationrad,
            raanrad: 0.0,
            argperigeerad: 0.0,
            orbitalanglerad: 0.0,
            orbittype,
        }
    }
    pub fn leo(name: &str, masskg: f64, altitudem: f64) -> Self {
        Self::new(name, masskg, altitudem, 0.0, 0.0)
    }
    pub fn geo(name: &str, masskg: f64) -> Self {
        Self::new(name, masskg, 35786000.0, 0.0, 0.0)
    }
    pub fn orbitalperiods(&self) -> f64 {
        2.0 * std::f64::consts::PI * (self.semimajoraxism.powi(3) / (G * EARTH_MASS)).sqrt()
    }
    pub fn orbitalvelocityms(&self) -> f64 {
        (G * EARTH_MASS / self.semimajoraxism).sqrt()
    }
    pub fn step(&mut self, dts: f64) {
        let period = self.orbitalperiods();
        let n = 2.0 * std::f64::consts::PI / period;
        let pi2 = 2.0 * std::f64::consts::PI;
        self.orbitalanglerad = (self.orbitalanglerad + n * dts) % pi2;
        let j2 = crate::J2EARTH;
        let p = self.semimajoraxism * (1.0 - self.eccentricity * self.eccentricity);
        let rratiosq = (EARTH_RADIUS / p).powi(2);
        let cosi = self.inclinationrad.cos();
        let sini = self.inclinationrad.sin();
        let raandot = -1.5 * n * j2 * rratiosq * cosi;
        self.raanrad = (self.raanrad + raandot * dts) % pi2;
        let omegadot = 1.5 * n * j2 * rratiosq * (2.0 - 2.5 * sini * sini);
        self.argperigeerad = (self.argperigeerad + omegadot * dts) % pi2;
        if self.altitudem < 1000000.0 {
            let scaleheight = 8500.0;
            let rho0 = *crate::SEALEVELAIRDENSITY;
            let rho = rho0 * (-self.altitudem / scaleheight).exp();
            let cdaoverm = 2.2 * 0.01;
            let v = self.orbitalvelocityms();
            let adrag = 0.5 * rho * v * v * cdaoverm;
            let da =
                -2.0 * self.semimajoraxism * self.semimajoraxism * adrag / (G * EARTH_MASS) * dts;
            self.semimajoraxism = (self.semimajoraxism + da).max(EARTH_RADIUS + 100000.0);
            self.altitudem = self.semimajoraxism - EARTH_RADIUS;
        }
    }
    pub fn position(&self) -> (f64, f64, f64) {
        let e = self.eccentricity;
        let nu = self.orbitalanglerad;
        let r = self.semimajoraxism * (1.0 - e * e) / (1.0 + e * nu.cos());
        let xorb = r * nu.cos();
        let yorb = r * nu.sin();
        let w = self.argperigeerad;
        let xw = xorb * w.cos() - yorb * w.sin();
        let yw = xorb * w.sin() + yorb * w.cos();
        let i = self.inclinationrad;
        let xi = xw;
        let yi = yw * i.cos();
        let zi = yw * i.sin();
        let omega = self.raanrad;
        let x = xi * omega.cos() - yi * omega.sin();
        let y = xi * omega.sin() + yi * omega.cos();
        let z = zi;
        (x, y, z)
    }
    pub fn gravityatsurface(&self) -> f64 {
        G * EARTH_MASS / self.semimajoraxism.powi(2)
    }
}
pub struct Constellation {
    pub name: String,
    pub satellites: Vec<ArtificialSatellite>,
}
impl Constellation {
    pub fn new(name: &str) -> Self {
        Self {
            name: name.to_string(),
            satellites: Vec::new(),
        }
    }
    pub fn add(&mut self, sat: ArtificialSatellite) {
        self.satellites.push(sat);
    }
    pub fn stepall(&mut self, dts: f64) {
        for sat in &mut self.satellites {
            sat.step(dts);
        }
    }
    pub fn positions(&self) -> Vec<(f64, f64, f64)> {
        self.satellites.iter().map(|s| s.position()).collect()
    }
}