use crate::{BodyId, CelestialBody, Vec3};
use sciforge::hub::domain::common::constants::G;
use std::f64::consts::PI;
#[derive(Clone, Debug)]
pub enum CelestialEvent {
Conjunction {
body_a: BodyId,
body_b: BodyId,
time_jd: f64,
angular_separation_rad: f64,
},
Opposition {
body: BodyId,
time_jd: f64,
},
Periapsis {
body: BodyId,
time_jd: f64,
distance: f64,
},
Apoapsis {
body: BodyId,
time_jd: f64,
distance: f64,
},
Eclipse {
eclipsed: BodyId,
eclipser: BodyId,
time_jd: f64,
},
CloseApproach {
body_a: BodyId,
body_b: BodyId,
time_jd: f64,
distance: f64,
},
HillSphereEntry {
body: BodyId,
intruder: BodyId,
time_jd: f64,
},
Transit {
transiting: BodyId,
observer: BodyId,
time_jd: f64,
angular_radius_rad: f64,
},
MutualEclipse {
eclipsed: BodyId,
eclipser: BodyId,
observer: BodyId,
time_jd: f64,
},
LaunchWindow {
from: BodyId,
to: BodyId,
time_jd: f64,
phase_angle_rad: f64,
delta_v: f64,
},
NodeCrossing {
body: BodyId,
time_jd: f64,
ascending: bool,
},
SyzygyAlignment {
bodies: [BodyId; 3],
time_jd: f64,
angular_deviation_rad: f64,
},
}
pub fn angular_separation(pos_a: Vec3, pos_b: Vec3, observer: Vec3) -> f64 {
let va = pos_a - observer;
let vb = pos_b - observer;
va.angle_between(&vb)
}
pub fn detect_conjunction(
sun: &CelestialBody,
a: &CelestialBody,
b: &CelestialBody,
threshold_rad: f64,
) -> bool {
let sep = angular_separation(a.position, b.position, sun.position);
sep < threshold_rad
}
pub fn detect_opposition(sun: &CelestialBody, earth: &CelestialBody, body: &CelestialBody) -> bool {
let to_sun = sun.position - earth.position;
let to_body = body.position - earth.position;
let angle = to_sun.angle_between(&to_body);
(angle - PI).abs() < 0.01
}
pub fn detect_periapsis_crossing(
body: &CelestialBody,
central: &CelestialBody,
prev_rdot: f64,
) -> Option<f64> {
let r_vec = body.position - central.position;
let v_vec = body.velocity - central.velocity;
let rdot = r_vec.dot(&v_vec) / r_vec.magnitude();
if prev_rdot < 0.0 && rdot >= 0.0 {
Some(r_vec.magnitude())
} else {
None
}
}
pub fn detect_apoapsis_crossing(
body: &CelestialBody,
central: &CelestialBody,
prev_rdot: f64,
) -> Option<f64> {
let r_vec = body.position - central.position;
let v_vec = body.velocity - central.velocity;
let rdot = r_vec.dot(&v_vec) / r_vec.magnitude();
if prev_rdot > 0.0 && rdot <= 0.0 {
Some(r_vec.magnitude())
} else {
None
}
}
pub fn detect_close_approach(a: &CelestialBody, b: &CelestialBody, threshold_m: f64) -> bool {
a.position.distance(&b.position) < threshold_m
}
pub fn detect_hill_sphere_entry(
body: &CelestialBody,
intruder: &CelestialBody,
parent_mass: f64,
semi_major: f64,
) -> bool {
let hill = body.hill_sphere(parent_mass, semi_major);
body.position.distance(&intruder.position) < hill
}
pub fn orbital_speed_circular(central_mass: f64, distance: f64) -> f64 {
(G * central_mass / distance).sqrt()
}
pub fn is_in_shadow(body: &CelestialBody, sun: &CelestialBody, occluder: &CelestialBody) -> bool {
let to_sun = sun.position - body.position;
let to_occluder = occluder.position - body.position;
let dist_occ = to_occluder.magnitude();
let dist_sun = to_sun.magnitude();
if dist_occ > dist_sun {
return false;
}
let angle = to_sun.angle_between(&to_occluder);
let angular_radius_occ = (occluder.radius / dist_occ).asin();
let angular_radius_sun = (sun.radius / dist_sun).asin();
angle < angular_radius_occ && angular_radius_occ >= angular_radius_sun
}
pub fn detect_transit(
transiting: &CelestialBody,
star: &CelestialBody,
observer: &CelestialBody,
) -> Option<f64> {
let to_star = star.position - observer.position;
let to_body = transiting.position - observer.position;
let dist_body = to_body.magnitude();
let dist_star = to_star.magnitude();
if dist_body >= dist_star {
return None;
}
let sep = to_star.angle_between(&to_body);
let angular_radius_star = (star.radius / dist_star).asin();
if sep < angular_radius_star {
Some((transiting.radius / dist_body).asin())
} else {
None
}
}
pub fn detect_mutual_eclipse(
eclipsed: &CelestialBody,
eclipser: &CelestialBody,
observer: &CelestialBody,
) -> bool {
let to_eclipsed = eclipsed.position - observer.position;
let to_eclipser = eclipser.position - observer.position;
let d_eclipsed = to_eclipsed.magnitude();
let d_eclipser = to_eclipser.magnitude();
if d_eclipser >= d_eclipsed {
return false;
}
let sep = to_eclipsed.angle_between(&to_eclipser);
let ar_eclipser = (eclipser.radius / d_eclipser).asin();
sep < ar_eclipser
}
pub fn detect_node_crossing(
body: &CelestialBody,
central: &CelestialBody,
prev_z: f64,
) -> Option<bool> {
let rel_z = body.position.z - central.position.z;
if prev_z < 0.0 && rel_z >= 0.0 {
Some(true)
} else if prev_z > 0.0 && rel_z <= 0.0 {
Some(false)
} else {
None
}
}
pub fn detect_syzygy(
a: &CelestialBody,
b: &CelestialBody,
c: &CelestialBody,
threshold_rad: f64,
) -> Option<f64> {
let ab = b.position - a.position;
let ac = c.position - a.position;
let angle = ab.angle_between(&ac);
let deviation = angle.min((PI - angle).abs());
if deviation < threshold_rad {
Some(deviation)
} else {
None
}
}
pub fn launch_window_phase(
departure: &CelestialBody,
arrival: &CelestialBody,
sun: &CelestialBody,
) -> f64 {
let r1 = departure.position - sun.position;
let r2 = arrival.position - sun.position;
let cos_phase = r1.dot(&r2) / (r1.magnitude() * r2.magnitude());
let sin_phase = r1.cross(&r2).z;
let phase = cos_phase.clamp(-1.0, 1.0).acos();
if sin_phase < 0.0 {
2.0 * PI - phase
} else {
phase
}
}
pub fn is_launch_window(
departure: &CelestialBody,
arrival: &CelestialBody,
sun: &CelestialBody,
mu: f64,
phase_tolerance: f64,
) -> Option<f64> {
let r1 = (departure.position - sun.position).magnitude();
let r2 = (arrival.position - sun.position).magnitude();
let required = crate::orbits::transfers::launch_window_phase_angle(r1, r2, mu);
let current = launch_window_phase(departure, arrival, sun);
if (current - required).abs() < phase_tolerance {
let transfer = crate::orbits::transfers::hohmann_transfer(r1, r2, mu);
Some(transfer.total_delta_v)
} else {
None
}
}
pub fn phase_angle_rate(body1_r: f64, body2_r: f64, mu: f64) -> f64 {
let omega1 = (mu / body1_r.powi(3)).sqrt();
let omega2 = (mu / body2_r.powi(3)).sqrt();
omega1 - omega2
}