venuss 0.0.2

Venus 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

use sciforge::hub::domain::astronomy::orbits::{kepler_period, kepler_velocity};
use sciforge::hub::domain::common::constants::{AU, G, SOLAR_MASS};

pub const VENUS_MASS: f64 = crate::VENUS_MASS;
pub const VENUS_RADIUS: f64 = crate::VENUS_RADIUS;
pub const SEMI_MAJOR_AXIS: f64 = crate::SEMI_MAJOR_AXIS_AU * AU;

pub struct VenusOrbit {
    pub semi_major_axis_m: f64,
    pub eccentricity: f64,
    pub inclination_deg: f64,
    pub longitude_ascending_node_deg: f64,
    pub argument_perihelion_deg: f64,
    pub mean_anomaly_rad: f64,
}

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

impl VenusOrbit {
    pub fn new() -> Self {
        Self {
            semi_major_axis_m: SEMI_MAJOR_AXIS,
            eccentricity: crate::ECCENTRICITY,
            inclination_deg: crate::INCLINATION_DEG,
            longitude_ascending_node_deg: crate::LONGITUDE_ASCENDING_NODE_DEG,
            argument_perihelion_deg: crate::ARGUMENT_PERIHELION_DEG,
            mean_anomaly_rad: crate::MEAN_ANOMALY_J2000_DEG.to_radians(),
        }
    }

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

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

    pub fn velocity_at_distance(&self, radius_m: f64) -> f64 {
        kepler_velocity(SOLAR_MASS, radius_m, self.semi_major_axis_m)
    }

    pub fn perihelion_m(&self) -> f64 {
        self.semi_major_axis_m * (1.0 - self.eccentricity)
    }

    pub fn aphelion_m(&self) -> f64 {
        self.semi_major_axis_m * (1.0 + self.eccentricity)
    }

    pub fn escape_velocity_at_surface() -> f64 {
        (2.0 * G * VENUS_MASS / VENUS_RADIUS).sqrt()
    }

    pub fn gravitational_force_sun(&self) -> f64 {
        sciforge::hub::domain::astronomy::celestial::gravitational_force(
            SOLAR_MASS,
            VENUS_MASS,
            self.current_radius(),
        )
    }

    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
    }

    pub fn true_anomaly_rad(&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 current_radius(&self) -> f64 {
        let nu = self.true_anomaly_rad();
        let e = self.eccentricity;
        self.semi_major_axis_m * (1.0 - e * e) / (1.0 + e * nu.cos())
    }

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

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

    pub fn position(&self) -> (f64, f64, f64) {
        let nu = self.true_anomaly_rad();
        let r = self.current_radius();
        let x_orb = r * nu.cos();
        let y_orb = r * nu.sin();

        let (sin_w, cos_w) = self.argument_perihelion_deg.to_radians().sin_cos();
        let (sin_o, cos_o) = self.longitude_ascending_node_deg.to_radians().sin_cos();
        let (sin_i, cos_i) = self.inclination_deg.to_radians().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 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);
    }
}