use crate::{CelestialBody, Vec3};
use sciforge::hub::domain::common::constants::G;
#[derive(Clone, Debug)]
pub struct StabilityMetrics {
pub max_energy_drift: f64,
pub mean_energy_drift: f64,
pub max_angular_momentum_drift: f64,
pub mean_angular_momentum_drift: f64,
pub lyapunov_exponent: f64,
pub megno: f64,
pub is_stable: bool,
}
fn separation(a: &CelestialBody, b: &CelestialBody) -> f64 {
a.position.distance(&b.position)
}
fn total_energy(bodies: &[CelestialBody]) -> f64 {
let kinetic: f64 = bodies.iter().map(|b| b.kinetic_energy()).sum();
let mut potential = 0.0;
for i in 0..bodies.len() {
for j in (i + 1)..bodies.len() {
let r = separation(&bodies[i], &bodies[j]);
if r > 0.0 {
potential -= G * bodies[i].mass * bodies[j].mass / r;
}
}
}
kinetic + potential
}
fn total_angular_momentum(bodies: &[CelestialBody]) -> Vec3 {
let mut l = Vec3::ZERO;
for b in bodies {
let li = b.position.cross(&(b.velocity * b.mass));
l += li;
}
l
}
pub fn energy_drift_analysis(snapshots: &[(f64, Vec<CelestialBody>)]) -> (f64, f64) {
if snapshots.is_empty() {
return (0.0, 0.0);
}
let e0 = total_energy(&snapshots[0].1);
if e0.abs() < 1e-30 {
return (0.0, 0.0);
}
let mut max_drift = 0.0_f64;
let mut sum_drift = 0.0;
for (_, bodies) in snapshots.iter().skip(1) {
let e = total_energy(bodies);
let drift = ((e - e0) / e0).abs();
max_drift = max_drift.max(drift);
sum_drift += drift;
}
let mean_drift = if snapshots.len() > 1 {
sum_drift / (snapshots.len() - 1) as f64
} else {
0.0
};
(max_drift, mean_drift)
}
pub fn angular_momentum_drift_analysis(snapshots: &[(f64, Vec<CelestialBody>)]) -> (f64, f64) {
if snapshots.is_empty() {
return (0.0, 0.0);
}
let l0 = total_angular_momentum(&snapshots[0].1);
let l0_mag = l0.magnitude();
if l0_mag < 1e-30 {
return (0.0, 0.0);
}
let mut max_drift = 0.0_f64;
let mut sum_drift = 0.0;
for (_, bodies) in snapshots.iter().skip(1) {
let l = total_angular_momentum(bodies);
let drift = (l - l0).magnitude() / l0_mag;
max_drift = max_drift.max(drift);
sum_drift += drift;
}
let mean_drift = if snapshots.len() > 1 {
sum_drift / (snapshots.len() - 1) as f64
} else {
0.0
};
(max_drift, mean_drift)
}
pub fn lyapunov_exponent(
snapshots_ref: &[(f64, Vec<CelestialBody>)],
snapshots_perturbed: &[(f64, Vec<CelestialBody>)],
) -> f64 {
if snapshots_ref.len() < 2 || snapshots_perturbed.len() < 2 {
return 0.0;
}
let n = snapshots_ref.len().min(snapshots_perturbed.len());
let d0 = initial_phase_distance(&snapshots_ref[0].1, &snapshots_perturbed[0].1);
if d0 < 1e-30 {
return 0.0;
}
let t0 = snapshots_ref[0].0;
let mut sum = 0.0;
let mut count = 0;
for i in 1..n {
let d = initial_phase_distance(&snapshots_ref[i].1, &snapshots_perturbed[i].1);
let dt = snapshots_ref[i].0 - t0;
if dt > 0.0 && d > 0.0 {
sum += (d / d0).ln() / dt;
count += 1;
}
}
if count > 0 { sum / count as f64 } else { 0.0 }
}
fn initial_phase_distance(a: &[CelestialBody], b: &[CelestialBody]) -> f64 {
let n = a.len().min(b.len());
let mut sum_sq = 0.0;
for i in 0..n {
let dp = a[i].position - b[i].position;
let dv = a[i].velocity - b[i].velocity;
sum_sq += dp.magnitude_squared() + dv.magnitude_squared();
}
sum_sq.sqrt()
}
pub fn megno_indicator(
snapshots_ref: &[(f64, Vec<CelestialBody>)],
snapshots_perturbed: &[(f64, Vec<CelestialBody>)],
) -> f64 {
if snapshots_ref.len() < 2 || snapshots_perturbed.len() < 2 {
return 0.0;
}
let n = snapshots_ref.len().min(snapshots_perturbed.len());
let t0 = snapshots_ref[0].0;
let d0 = initial_phase_distance(&snapshots_ref[0].1, &snapshots_perturbed[0].1);
if d0 < 1e-30 {
return 0.0;
}
let mut y_sum = 0.0;
let mut count = 0;
for i in 1..n {
let d = initial_phase_distance(&snapshots_ref[i].1, &snapshots_perturbed[i].1);
let t = snapshots_ref[i].0 - t0;
if t > 0.0 && d > 0.0 {
let d_dot_approx = if i > 1 {
let d_prev =
initial_phase_distance(&snapshots_ref[i - 1].1, &snapshots_perturbed[i - 1].1);
let dt_step = snapshots_ref[i].0 - snapshots_ref[i - 1].0;
if dt_step > 0.0 {
(d - d_prev) / dt_step
} else {
0.0
}
} else {
(d - d0) / t
};
y_sum += t * d_dot_approx / d;
count += 1;
}
}
if count == 0 {
return 0.0;
}
let total_time = snapshots_ref[n - 1].0 - t0;
if total_time <= 0.0 {
return 0.0;
}
2.0 * y_sum / (count as f64 * total_time)
}
pub fn hill_stability_criterion(
m_central: f64,
m1: f64,
m2: f64,
a1: f64,
a2: f64,
e1: f64,
e2: f64,
) -> bool {
let mu1 = m1 / m_central;
let mu2 = m2 / m_central;
let alpha = a1 / a2;
let gamma = mu1 * mu2 / ((mu1 + mu2) * (mu1 + mu2));
let delta = (a2 * (1.0 - e2) - a1 * (1.0 + e1)) / a1;
let critical = 2.0 * 3.0_f64.sqrt() * (gamma * alpha).powf(1.0 / 3.0);
delta > critical
}
pub fn mardling_aarseth_stability(
a_outer: f64,
e_outer: f64,
a_inner: f64,
q_out: f64,
q_in: f64,
) -> bool {
let ratio = a_outer * (1.0 - e_outer) / a_inner;
let mass_factor = (1.0 + q_out).powf(2.0 / 5.0) * (1.0 + q_in).powf(1.0 / 5.0);
let critical =
2.8 * mass_factor * ((1.0 + e_outer) / (1.0 - e_outer).powf(1.2)).powf(2.0 / 5.0);
ratio > critical
}
pub fn full_stability_analysis(
snapshots_ref: &[(f64, Vec<CelestialBody>)],
snapshots_perturbed: &[(f64, Vec<CelestialBody>)],
) -> StabilityMetrics {
let (max_e, mean_e) = energy_drift_analysis(snapshots_ref);
let (max_l, mean_l) = angular_momentum_drift_analysis(snapshots_ref);
let lce = lyapunov_exponent(snapshots_ref, snapshots_perturbed);
let meg = megno_indicator(snapshots_ref, snapshots_perturbed);
let stable = max_e < 1e-6 && max_l < 1e-6 && lce < 1e-4 && (meg - 2.0).abs() < 0.5;
StabilityMetrics {
max_energy_drift: max_e,
mean_energy_drift: mean_e,
max_angular_momentum_drift: max_l,
mean_angular_momentum_drift: mean_l,
lyapunov_exponent: lce,
megno: meg,
is_stable: stable,
}
}
pub fn propagate_centuries(
mut bodies: Vec<CelestialBody>,
dt: f64,
centuries: f64,
snapshot_interval: f64,
) -> Vec<(f64, Vec<CelestialBody>)> {
let total_time = centuries * 365.25 * 86400.0 * 100.0;
let steps = (total_time / dt) as u64;
let snap_steps = (snapshot_interval / dt).max(1.0) as u64;
let mut snapshots = Vec::new();
snapshots.push((0.0, bodies.clone()));
for step in 1..=steps {
leapfrog_step(&mut bodies, dt);
if step % snap_steps == 0 {
let t = step as f64 * dt;
snapshots.push((t, bodies.clone()));
}
}
snapshots
}
fn leapfrog_step(bodies: &mut [CelestialBody], dt: f64) {
let n = bodies.len();
let half_dt = 0.5 * dt;
for body in bodies.iter_mut() {
body.velocity += body.acceleration * half_dt;
}
for body in bodies.iter_mut() {
body.position += body.velocity * dt;
}
let mut accels = vec![Vec3::ZERO; n];
for i in 0..n {
for j in 0..n {
if i == j {
continue;
}
let r = bodies[j].position - bodies[i].position;
let dist = r.magnitude();
if dist > 0.0 {
let a = G * bodies[j].mass / (dist * dist * dist);
accels[i] += r * a;
}
}
}
for i in 0..n {
bodies[i].acceleration = accels[i];
bodies[i].velocity += bodies[i].acceleration * half_dt;
}
}
pub fn detect_resonance(period1: f64, period2: f64, max_order: u32) -> Option<(u32, u32)> {
let ratio = period1 / period2;
let mut best: Option<(u32, u32, f64)> = None;
for p in 1..=max_order {
for q in 1..=max_order {
let r = p as f64 / q as f64;
let err = (ratio - r).abs() / ratio;
if err < 0.01 {
match best {
Some((_, _, prev_err)) if err < prev_err => {
best = Some((p, q, err));
}
None => {
best = Some((p, q, err));
}
_ => {}
}
}
}
}
best.map(|(p, q, _)| (p, q))
}
pub fn secular_frequency(snapshots: &[(f64, Vec<CelestialBody>)], body_idx: usize) -> f64 {
if snapshots.len() < 3 || body_idx >= snapshots[0].1.len() {
return 0.0;
}
let mut crossings = Vec::new();
let mut prev_z = snapshots[0].1[body_idx].position.z;
for (t, bodies) in snapshots.iter().skip(1) {
let z = bodies[body_idx].position.z;
if prev_z < 0.0 && z >= 0.0 {
crossings.push(*t);
}
prev_z = z;
}
if crossings.len() < 2 {
return 0.0;
}
let total_span = crossings[crossings.len() - 1] - crossings[0];
let num_periods = (crossings.len() - 1) as f64;
if total_span <= 0.0 {
0.0
} else {
2.0 * std::f64::consts::PI * num_periods / total_span
}
}