use crate::CelestialBody;
use crate::dynamics::events::CelestialEvent;
use crate::dynamics::integrator::Integrator;
use crate::dynamics::time::EphemerisTime;
use crate::gravity::nbody::compute_accelerations;
use sciforge::hub::domain::common::constants::{AU, SOLAR_MASS};
use std::time::{Duration, Instant};
#[derive(Debug, Clone, Copy)]
pub struct RunBudget {
pub max_steps: Option<usize>,
pub max_runtime_ms: Option<u64>,
pub yield_every: usize,
}
impl Default for RunBudget {
fn default() -> Self {
Self {
max_steps: None,
max_runtime_ms: None,
yield_every: 1024,
}
}
}
pub struct Pipeline {
pub bodies: Vec<CelestialBody>,
pub time: EphemerisTime,
pub integrator: Integrator,
pub dt: f64,
pub events: Vec<CelestialEvent>,
prev_rdots: Vec<f64>,
}
impl Pipeline {
pub fn new(bodies: Vec<CelestialBody>, dt: f64, integrator: Integrator) -> Self {
let prev_rdots = vec![0.0; bodies.len()];
Self {
bodies,
time: EphemerisTime::j2000(),
integrator,
dt,
events: Vec::new(),
prev_rdots,
}
}
pub fn initialize(&mut self) {
compute_accelerations(&mut self.bodies);
self.update_rdots();
}
pub fn step(&mut self) {
self.integrator.step(&mut self.bodies, self.dt);
self.time.advance(self.dt);
self.detect_events();
self.update_rdots();
}
pub fn run(&mut self, steps: usize) {
let _ = self.run_with_budget(steps, RunBudget::default());
}
pub fn run_with_budget(&mut self, steps: usize, budget: RunBudget) -> usize {
self.initialize();
let capped_steps = budget.max_steps.map_or(steps, |m| m.min(steps));
let deadline = budget
.max_runtime_ms
.map(|ms| Instant::now() + Duration::from_millis(ms));
let mut executed = 0usize;
for idx in 0..capped_steps {
if let Some(t) = deadline
&& Instant::now() >= t
{
break;
}
self.step();
executed += 1;
if budget.yield_every > 0 && (idx + 1).is_multiple_of(budget.yield_every) {
std::thread::yield_now();
}
}
executed
}
pub fn run_until(&mut self, target_jd: f64) {
let _ = self.run_until_with_budget(target_jd, RunBudget::default());
}
pub fn run_until_with_budget(&mut self, target_jd: f64, budget: RunBudget) -> usize {
self.initialize();
let deadline = budget
.max_runtime_ms
.map(|ms| Instant::now() + Duration::from_millis(ms));
let mut executed = 0usize;
while self.time.current_jd() < target_jd {
if let Some(max_steps) = budget.max_steps
&& executed >= max_steps
{
break;
}
if let Some(t) = deadline
&& Instant::now() >= t
{
break;
}
self.step();
executed += 1;
if budget.yield_every > 0 && executed.is_multiple_of(budget.yield_every) {
std::thread::yield_now();
}
}
executed
}
fn update_rdots(&mut self) {
for (i, body) in self.bodies.iter().enumerate() {
if let Some(parent_id) = body.id.parent()
&& let Some(parent) = self.bodies.iter().find(|b| b.id == parent_id)
{
let r = body.position - parent.position;
let v = body.velocity - parent.velocity;
self.prev_rdots[i] = r.dot(&v) / r.magnitude();
}
}
}
fn detect_events(&mut self) {
use crate::BodyId;
use crate::dynamics::events::*;
let snapshot = self.bodies.clone();
let jd = self.time.current_jd();
for (i, body) in snapshot.iter().enumerate() {
if let Some(parent_id) = body.id.parent()
&& let Some(parent) = snapshot.iter().find(|b| b.id == parent_id)
{
if let Some(dist) = detect_periapsis_crossing(body, parent, self.prev_rdots[i]) {
self.events.push(CelestialEvent::Periapsis {
body: body.id,
time_jd: jd,
distance: dist,
});
}
if let Some(dist) = detect_apoapsis_crossing(body, parent, self.prev_rdots[i]) {
self.events.push(CelestialEvent::Apoapsis {
body: body.id,
time_jd: jd,
distance: dist,
});
}
}
}
let mut trackable: Vec<BodyId> = Vec::new();
trackable.extend_from_slice(BodyId::all_planets());
trackable.extend_from_slice(BodyId::all_dwarf_planets());
trackable.extend_from_slice(BodyId::all_asteroids());
trackable.extend_from_slice(BodyId::all_comets());
for (pi, &pa) in trackable.iter().enumerate() {
for &pb in &trackable[pi + 1..] {
if let (Some(a), Some(b)) = (
snapshot.iter().find(|b| b.id == pa),
snapshot.iter().find(|b| b.id == pb),
) {
let dist = a.position.distance(&b.position);
let sun_mass = snapshot
.iter()
.find(|b| b.id == BodyId::Sun)
.map_or(SOLAR_MASS, |s| s.mass);
let a_r = a.position.magnitude().max(1.0);
let b_r = b.position.magnitude().max(1.0);
let a_hill = a.hill_sphere(sun_mass, a_r);
let b_hill = b.hill_sphere(sun_mass, b_r);
let a_soi = a.sphere_of_influence(sun_mass, a_r);
let b_soi = b.sphere_of_influence(sun_mass, b_r);
let threshold = (a_hill + b_hill)
.max(0.25 * (a_soi + b_soi))
.clamp(0.005 * AU, 0.5 * AU);
if dist < threshold {
self.events.push(CelestialEvent::CloseApproach {
body_a: pa,
body_b: pb,
time_jd: jd,
distance: dist,
});
}
}
}
}
}
pub fn total_energy(&self) -> f64 {
crate::config::parameters::total_energy(&self.bodies)
}
pub fn total_angular_momentum(&self) -> crate::Vec3 {
crate::config::parameters::total_angular_momentum(&self.bodies)
}
}