use crate::config::parameters;
pub fn resonance_semi_major_axis(p: u32, q: u32) -> f64 {
if p == 0 || q == 0 {
return 0.0;
}
let ratio = (q as f64) / (p as f64);
parameters::NEPTUNE_A * ratio.powf(2.0 / 3.0)
}
pub fn is_near_resonance(semi_major_au: f64, p: u32, q: u32, tolerance_au: f64) -> bool {
let a_res = resonance_semi_major_axis(p, q);
(semi_major_au - a_res).abs() < tolerance_au
}
pub fn resonance_width(p: u32, q: u32, eccentricity: f64) -> f64 {
let a_res = resonance_semi_major_axis(p, q);
let mass_ratio = parameters::NEPTUNE_MASS / parameters::SOLAR_MASS;
a_res * (mass_ratio * eccentricity.abs().max(0.001)).sqrt()
}
pub fn libration_period_years(p: u32, q: u32, eccentricity: f64) -> f64 {
let a_res = resonance_semi_major_axis(p, q);
let t_orb = parameters::orbital_period(a_res, parameters::SOLAR_MASS) / parameters::YEAR;
let mass_ratio = parameters::NEPTUNE_MASS / parameters::SOLAR_MASS;
if mass_ratio < 1.0e-20 || eccentricity < 1.0e-10 {
return f64::INFINITY;
}
t_orb / (mass_ratio * eccentricity).sqrt()
}
pub fn major_resonances() -> Vec<(u32, u32, f64)> {
let resonances = [(2, 3), (1, 2), (3, 5), (4, 7), (2, 5), (1, 3)];
resonances
.iter()
.map(|&(p, q)| (p, q, resonance_semi_major_axis(p, q)))
.collect()
}
pub fn closest_resonance(semi_major_au: f64, tolerance_au: f64) -> Option<(u32, u32, f64)> {
major_resonances()
.into_iter()
.filter(|&(_, _, a_res)| (semi_major_au - a_res).abs() < tolerance_au)
.min_by(|a, b| {
(semi_major_au - a.2)
.abs()
.partial_cmp(&(semi_major_au - b.2).abs())
.unwrap_or(std::cmp::Ordering::Equal)
})
}