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
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
140
141
142
143
144
145
146
147
148
149

use sciforge::hub::domain::astronomy::celestial::gravitational_force;
use sciforge::hub::domain::astronomy::orbits::{
    angular_momentum, apoapsis, escape_velocity, kepler_period, orbital_energy, periapsis, vis_viva,
};
use sciforge::hub::domain::common::constants::{AU, G, SOLAR_MASS};

pub const MARS_MASS: f64 = crate::MARS_MASS;
pub const MARS_RADIUS: f64 = crate::MARS_RADIUS;
pub const SEMI_MAJOR_AXIS: f64 = crate::SEMI_MAJOR_AXIS_AU * AU;
pub const ECCENTRICITY: f64 = crate::ECCENTRICITY;
pub const MEAN_ANOMALY_J2000_DEG: f64 = crate::MEAN_ANOMALY_J2000_DEG;

const MU_SUN: f64 = G * SOLAR_MASS;

pub struct MarsOrbit {
    pub semi_major_axis_m: f64,
    pub eccentricity: f64,
    pub true_anomaly_rad: f64,
    pub mean_anomaly_rad: f64,
}

impl Default for MarsOrbit {
    fn default() -> Self {
        Self::new()
    }
}

impl MarsOrbit {
    pub fn new() -> Self {
        let m0 = MEAN_ANOMALY_J2000_DEG.to_radians();
        let mut s = Self {
            semi_major_axis_m: SEMI_MAJOR_AXIS,
            eccentricity: ECCENTRICITY,
            true_anomaly_rad: 0.0,
            mean_anomaly_rad: m0,
        };
        s.solve_kepler();
        s
    }

    pub fn at_mean_anomaly(m_rad: f64) -> Self {
        let mut s = Self {
            semi_major_axis_m: SEMI_MAJOR_AXIS,
            eccentricity: ECCENTRICITY,
            true_anomaly_rad: 0.0,
            mean_anomaly_rad: m_rad,
        };
        s.solve_kepler();
        s
    }

    fn solve_kepler(&mut self) {
        let m = self.mean_anomaly_rad;
        let e = self.eccentricity;
        let mut ea = m + e * m.sin();
        for _ in 0..15 {
            let f = ea - e * ea.sin() - m;
            let fp = 1.0 - e * ea.cos();
            ea -= f / fp;
        }
        self.true_anomaly_rad = 2.0
            * f64::atan2(
                (1.0 + e).sqrt() * (ea / 2.0).sin(),
                (1.0 - e).sqrt() * (ea / 2.0).cos(),
            );
    }

    pub fn orbital_period_s(&self) -> f64 {
        kepler_period(self.semi_major_axis_m, MU_SUN)
    }

    pub fn orbital_period_days(&self) -> f64 {
        self.orbital_period_s() / 86400.0
    }

    pub fn velocity_at_distance(&self, r: f64) -> f64 {
        vis_viva(MU_SUN, r, self.semi_major_axis_m)
    }

    pub fn perihelion_m(&self) -> f64 {
        periapsis(self.semi_major_axis_m, self.eccentricity)
    }

    pub fn aphelion_m(&self) -> f64 {
        apoapsis(self.semi_major_axis_m, self.eccentricity)
    }

    pub fn specific_orbital_energy(&self) -> f64 {
        orbital_energy(MU_SUN, self.semi_major_axis_m)
    }

    pub fn specific_angular_momentum(&self) -> f64 {
        angular_momentum(MU_SUN, self.semi_major_axis_m, self.eccentricity)
    }

    pub fn escape_velocity_at_surface() -> f64 {
        escape_velocity(G * MARS_MASS, MARS_RADIUS)
    }

    pub fn gravitational_force_sun(&self) -> f64 {
        let r = self.current_radius();
        gravitational_force(SOLAR_MASS, MARS_MASS, r)
    }

    pub fn current_radius(&self) -> f64 {
        let a = self.semi_major_axis_m;
        let e = self.eccentricity;
        a * (1.0 - e * e) / (1.0 + e * self.true_anomaly_rad.cos())
    }

    pub fn mean_orbital_velocity(&self) -> f64 {
        let period = self.orbital_period_s();
        2.0 * std::f64::consts::PI * self.semi_major_axis_m / period
    }

    pub fn true_anomaly_deg(&self) -> f64 {
        self.true_anomaly_rad.to_degrees()
    }

    pub fn position(&self) -> (f64, f64, f64) {
        let r = self.current_radius();
        let nu = self.true_anomaly_rad;
        let x_orb = r * nu.cos();
        let y_orb = r * nu.sin();
        let w = crate::ARGUMENT_PERIHELION_DEG.to_radians();
        let x_w = x_orb * w.cos() - y_orb * w.sin();
        let y_w = x_orb * w.sin() + y_orb * w.cos();
        let i = crate::INCLINATION_DEG.to_radians();
        let x_i = x_w;
        let y_i = y_w * i.cos();
        let z_i = y_w * i.sin();
        let omega = crate::LONGITUDE_ASCENDING_NODE_DEG.to_radians();
        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 solar_irradiance(&self) -> f64 {
        let r_au = self.current_radius() / AU;
        sciforge::hub::domain::meteorology::radiation::solar_constant() / (r_au * r_au)
    }

    pub fn step(&mut self, dt_s: f64) {
        let n = 2.0 * std::f64::consts::PI / self.orbital_period_s();
        self.mean_anomaly_rad = (self.mean_anomaly_rad + n * dt_s) % (2.0 * std::f64::consts::PI);
        self.solve_kepler();
    }
}