solarsystems 0.0.1

N-body solar system engine — gravitational dynamics, orbital mechanics, perturbations, event detection, and full celestial orchestration
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

use crate::Vec3;
use sciforge::hub::domain::common::constants::G;
use std::f64::consts::PI;

pub fn j2_secular_precession_raan(j2: f64, re: f64, a: f64, e: f64, inc: f64, mu: f64) -> f64 {
    let n = (mu / a.powi(3)).sqrt();
    let p = a * (1.0 - e * e);
    -1.5 * n * j2 * (re / p).powi(2) * inc.cos()
}

pub fn j2_secular_precession_argp(j2: f64, re: f64, a: f64, e: f64, inc: f64, mu: f64) -> f64 {
    let n = (mu / a.powi(3)).sqrt();
    let p = a * (1.0 - e * e);
    0.75 * n * j2 * (re / p).powi(2) * (5.0 * inc.cos().powi(2) - 1.0)
}

pub fn mercury_gr_precession_rate(a: f64, e: f64, mu: f64) -> f64 {
    let c = 2.99792458e8;
    let t = 2.0 * PI * (a.powi(3) / mu).sqrt();
    6.0 * PI * mu / (a * (1.0 - e * e) * c * c * t)
}

pub fn solar_radiation_pressure_accel(
    luminosity: f64,
    area_to_mass: f64,
    reflectivity: f64,
    distance: f64,
) -> f64 {
    let c = 2.99792458e8;
    let flux = luminosity / (4.0 * PI * distance * distance);
    (1.0 + reflectivity) * flux * area_to_mass / c
}

pub fn atmospheric_drag_acceleration(
    density: f64,
    velocity_mag: f64,
    cd: f64,
    area_to_mass: f64,
) -> f64 {
    -0.5 * cd * area_to_mass * density * velocity_mag * velocity_mag
}

pub fn third_body_perturbation(body_pos: Vec3, perturber_pos: Vec3, perturber_mass: f64) -> Vec3 {
    let r_body = body_pos - perturber_pos;
    let r_body_mag = r_body.magnitude();
    let r_pert_mag = perturber_pos.magnitude();
    -(r_body / r_body_mag.powi(3) + perturber_pos / r_pert_mag.powi(3)) * (G * perturber_mass)
}

pub fn yarkovsky_acceleration(
    thermal_conductivity: f64,
    albedo: f64,
    radius: f64,
    density: f64,
    distance_to_sun: f64,
    solar_luminosity: f64,
) -> f64 {
    let c = 2.99792458e8;
    let flux = solar_luminosity / (4.0 * PI * distance_to_sun * distance_to_sun);
    let absorbed = (1.0 - albedo) * flux;
    let area = PI * radius * radius;
    let mass = (4.0 / 3.0) * PI * radius.powi(3) * density;
    let asymmetry = 0.5 * (thermal_conductivity / (thermal_conductivity + 1.0));
    asymmetry * absorbed * area / (mass * c)
}

pub fn poynting_robertson_drag(
    solar_luminosity: f64,
    particle_radius: f64,
    particle_density: f64,
    distance: f64,
    velocity_radial: f64,
) -> f64 {
    let c = 2.99792458e8;
    let area = PI * particle_radius * particle_radius;
    let mass = (4.0 / 3.0) * PI * particle_radius.powi(3) * particle_density;
    let flux = solar_luminosity / (4.0 * PI * distance * distance);
    -flux * area * velocity_radial / (mass * c * c)
}

pub fn kozai_lidov_max_eccentricity(inclination: f64) -> f64 {
    (1.0 - (5.0 / 3.0) * inclination.cos().powi(2)).sqrt()
}

pub fn lense_thirring_precession(angular_momentum_central: f64, a: f64, mu: f64) -> f64 {
    let c = 2.99792458e8;
    2.0 * G * angular_momentum_central / (c * c * a.powi(3) * (mu / a.powi(3)).sqrt())
}

pub fn tidal_orbital_decay_rate(
    satellite_mass: f64,
    central_mass: f64,
    central_radius: f64,
    semi_major_axis: f64,
    quality_factor: f64,
    love_number: f64,
) -> f64 {
    let mu_ratio = satellite_mass / central_mass;
    let r_ratio = central_radius / semi_major_axis;
    let n = (G * central_mass / semi_major_axis.powi(3)).sqrt();
    -3.0 * love_number * mu_ratio * r_ratio.powi(5) * n * semi_major_axis / quality_factor
}

pub fn dm_perihelion_precession(orbital_radius: f64, orbital_velocity: f64) -> f64 {
    crate::darkmatter::perturbations::perihelion_precession_rate(orbital_radius, orbital_velocity)
}

pub fn dm_orbital_velocity_correction(orbital_radius: f64, m_central: f64) -> f64 {
    crate::darkmatter::perturbations::orbital_velocity_correction(orbital_radius, m_central)
}

pub fn dm_secular_eccentricity_change(a: f64, m_central: f64, e: f64, age_s: f64) -> f64 {
    crate::darkmatter::perturbations::secular_eccentricity_change(a, m_central, e, age_s)
}

pub fn dm_dynamical_friction_on_planet(m_planet: f64, v_planet: f64) -> f64 {
    crate::darkmatter::perturbations::dynamical_friction_on_planet(m_planet, v_planet)
}

pub fn dm_orbital_decay_timescale(m_planet: f64, v_planet: f64) -> f64 {
    crate::darkmatter::perturbations::orbital_decay_timescale(m_planet, v_planet)
}

pub fn dm_acceleration_at_radius(r: f64) -> f64 {
    crate::darkmatter::perturbations::dm_acceleration_at(r)
}