earths 0.0.1

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

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 mass_kg: f64,
    pub altitude_m: f64,
    pub semi_major_axis_m: f64,
    pub eccentricity: f64,
    pub inclination_rad: f64,
    pub raan_rad: f64,
    pub arg_perigee_rad: f64,
    pub orbital_angle_rad: f64,
    pub orbit_type: OrbitType,
}
impl ArtificialSatellite {
    pub fn new(
        name: &str,
        mass_kg: f64,
        altitude_m: f64,
        eccentricity: f64,
        inclination_rad: f64,
    ) -> Self {
        let semi_major_axis_m = EARTH_RADIUS + altitude_m;
        let orbit_type = match altitude_m {
            a if a < 2_000_000.0 => OrbitType::LEO,
            a if a < 20_200_000.0 => OrbitType::MEO,
            a if (35_786_000.0 - a).abs() < 100_000.0 => OrbitType::GEO,
            a if a > 35_786_000.0 => OrbitType::HEO,
            _ => OrbitType::Custom,
        };
        Self {
            name: name.to_string(),
            mass_kg,
            altitude_m,
            semi_major_axis_m,
            eccentricity,
            inclination_rad,
            raan_rad: 0.0,
            arg_perigee_rad: 0.0,
            orbital_angle_rad: 0.0,
            orbit_type,
        }
    }
    pub fn leo(name: &str, mass_kg: f64, altitude_m: f64) -> Self {
        Self::new(name, mass_kg, altitude_m, 0.0, 0.0)
    }
    pub fn geo(name: &str, mass_kg: f64) -> Self {
        Self::new(name, mass_kg, 35_786_000.0, 0.0, 0.0)
    }
    pub fn orbital_period_s(&self) -> f64 {
        2.0 * std::f64::consts::PI * (self.semi_major_axis_m.powi(3) / (G * EARTH_MASS)).sqrt()
    }
    pub fn orbital_velocity_ms(&self) -> f64 {
        (G * EARTH_MASS / self.semi_major_axis_m).sqrt()
    }
    pub fn step(&mut self, dt_s: f64) {
        let period = self.orbital_period_s();
        let n = 2.0 * std::f64::consts::PI / period;
        let pi2 = 2.0 * std::f64::consts::PI;
        self.orbital_angle_rad = (self.orbital_angle_rad + n * dt_s) % pi2;
        let j2 = crate::J2_EARTH;
        let p = self.semi_major_axis_m * (1.0 - self.eccentricity * self.eccentricity);
        let r_ratio_sq = (EARTH_RADIUS / p).powi(2);
        let cos_i = self.inclination_rad.cos();
        let sin_i = self.inclination_rad.sin();
        let raan_dot = -1.5 * n * j2 * r_ratio_sq * cos_i;
        self.raan_rad = (self.raan_rad + raan_dot * dt_s) % pi2;
        let omega_dot = 1.5 * n * j2 * r_ratio_sq * (2.0 - 2.5 * sin_i * sin_i);
        self.arg_perigee_rad = (self.arg_perigee_rad + omega_dot * dt_s) % pi2;
        if self.altitude_m < 1_000_000.0 {
            let scale_height = 8_500.0;
            let rho_0 = *crate::SEA_LEVEL_AIR_DENSITY;
            let rho = rho_0 * (-self.altitude_m / scale_height).exp();
            let cd_a_over_m = 2.2 * 0.01;
            let v = self.orbital_velocity_ms();
            let a_drag = 0.5 * rho * v * v * cd_a_over_m;
            let da = -2.0 * self.semi_major_axis_m * self.semi_major_axis_m * a_drag
                / (G * EARTH_MASS)
                * dt_s;
            self.semi_major_axis_m = (self.semi_major_axis_m + da).max(EARTH_RADIUS + 100_000.0);
            self.altitude_m = self.semi_major_axis_m - EARTH_RADIUS;
        }
    }
    pub fn position(&self) -> (f64, f64, f64) {
        let e = self.eccentricity;
        let nu = self.orbital_angle_rad;
        let r = self.semi_major_axis_m * (1.0 - e * e) / (1.0 + e * nu.cos());
        let x_orb = r * nu.cos();
        let y_orb = r * nu.sin();
        let w = self.arg_perigee_rad;
        let x_w = x_orb * w.cos() - y_orb * w.sin();
        let y_w = x_orb * w.sin() + y_orb * w.cos();
        let i = self.inclination_rad;
        let x_i = x_w;
        let y_i = y_w * i.cos();
        let z_i = y_w * i.sin();
        let omega = self.raan_rad;
        let x = x_i * omega.cos() - y_i * omega.sin();
        let y = x_i * omega.sin() + y_i * omega.cos();
        let z = z_i;
        (x, y, z)
    }
    pub fn gravity_at_surface(&self) -> f64 {
        G * EARTH_MASS / self.semi_major_axis_m.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 step_all(&mut self, dt_s: f64) {
        for sat in &mut self.satellites {
            sat.step(dt_s);
        }
    }
    pub fn positions(&self) -> Vec<(f64, f64, f64)> {
        self.satellites.iter().map(|s| s.position()).collect()
    }
}