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

use sciforge::hub::domain::common::constants::G;

pub const PHOBOS_MASS: f64 = 1.0659e16;
pub const PHOBOS_RADIUS: f64 = 11_267.0;
pub const PHOBOS_DIMENSIONS: (f64, f64, f64) = (27_000.0, 22_000.0, 18_000.0);
pub const PHOBOS_SEMI_MAJOR_AXIS: f64 = 9_376_000.0;
pub const PHOBOS_ECCENTRICITY: f64 = 0.0151;
pub const PHOBOS_INCLINATION_DEG: f64 = 1.093;
pub const PHOBOS_ORBITAL_PERIOD_S: f64 = 27_554.0;
pub const PHOBOS_DENSITY: f64 = 1_876.0;
pub const PHOBOS_ALBEDO: f64 = 0.071;
pub const PHOBOS_SURFACE_GRAVITY: f64 = 0.0057;

pub struct PhobosState {
    pub mass_kg: f64,
    pub semi_major_axis_m: f64,
    pub eccentricity: f64,
    pub inclination_rad: f64,
    pub mean_anomaly_rad: f64,
    pub raan_rad: f64,
    pub arg_perigee_rad: f64,
}

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

impl PhobosState {
    pub fn new() -> Self {
        Self {
            mass_kg: PHOBOS_MASS,
            semi_major_axis_m: PHOBOS_SEMI_MAJOR_AXIS,
            eccentricity: PHOBOS_ECCENTRICITY,
            inclination_rad: PHOBOS_INCLINATION_DEG.to_radians(),
            mean_anomaly_rad: 0.0,
            raan_rad: 0.0,
            arg_perigee_rad: 0.0,
        }
    }

    pub fn orbital_period_s(&self) -> f64 {
        let mu = G * crate::MARS_MASS;
        2.0 * std::f64::consts::PI * (self.semi_major_axis_m.powi(3) / mu).sqrt()
    }

    pub fn orbital_velocity_ms(&self) -> f64 {
        (G * crate::MARS_MASS / self.semi_major_axis_m).sqrt()
    }

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

        let j2 = crate::J2_MARS;
        let p = self.semi_major_axis_m * (1.0 - self.eccentricity * self.eccentricity);
        let r_ratio_sq = (crate::MARS_EQUATORIAL_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;
    }

    fn eccentric_anomaly(&self) -> f64 {
        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;
        }
        ea
    }

    fn true_anomaly(&self) -> f64 {
        let ea = self.eccentric_anomaly();
        let e = self.eccentricity;
        2.0 * f64::atan2(
            (1.0 + e).sqrt() * (ea / 2.0).sin(),
            (1.0 - e).sqrt() * (ea / 2.0).cos(),
        )
    }

    pub fn distance(&self) -> f64 {
        let nu = self.true_anomaly();
        let e = self.eccentricity;
        self.semi_major_axis_m * (1.0 - e * e) / (1.0 + e * nu.cos())
    }

    pub fn position(&self) -> (f64, f64, f64) {
        let nu = self.true_anomaly();
        let r = self.distance();
        let (sin_nu, cos_nu) = nu.sin_cos();

        let x_orb = r * cos_nu;
        let y_orb = r * sin_nu;

        let (sin_w, cos_w) = self.arg_perigee_rad.sin_cos();
        let (sin_o, cos_o) = self.raan_rad.sin_cos();
        let (sin_i, cos_i) = self.inclination_rad.sin_cos();

        let x = (cos_w * cos_o - sin_w * sin_o * cos_i) * x_orb
            + (-sin_w * cos_o - cos_w * sin_o * cos_i) * y_orb;
        let y = (cos_w * sin_o + sin_w * cos_o * cos_i) * x_orb
            + (-sin_w * sin_o + cos_w * cos_o * cos_i) * y_orb;
        let z = (sin_w * sin_i) * x_orb + (cos_w * sin_i) * y_orb;
        (x, y, z)
    }

    pub fn escape_velocity(&self) -> f64 {
        (2.0 * G * PHOBOS_MASS / PHOBOS_RADIUS).sqrt()
    }
}

pub fn roche_limit_rigid() -> f64 {
    let rho_mars =
        crate::MARS_MASS / (4.0 / 3.0 * std::f64::consts::PI * crate::MARS_RADIUS.powi(3));
    crate::MARS_RADIUS * (2.0 * rho_mars / PHOBOS_DENSITY).cbrt()
}

pub fn apparent_angular_size_rad() -> f64 {
    use sciforge::hub::domain::astronomy::celestial::apparent_angular_size;
    apparent_angular_size(PHOBOS_RADIUS, PHOBOS_SEMI_MAJOR_AXIS - crate::MARS_RADIUS)
}

pub fn orbital_decay_rate_m_per_year() -> f64 {
    1.8e-2
}

pub fn stickney_crater_diameter_m() -> f64 {
    9_000.0
}