use sciforge::hub::domain::common::constants::AU;
use solarsystems::config::parameters::{
DEFAULT_TIMESTEP_S, GRAVITATIONAL_SOFTENING, J2_EARTH, J2_JUPITER, J2_MARS, J2_SATURN, J2_SUN,
MU_SUN, default_bodies, default_bodies_with_belts, find_body, find_body_mut, system_barycenter,
total_angular_momentum, total_energy,
};
use solarsystems::dynamics::events::{
angular_separation, detect_apoapsis_crossing, detect_close_approach, detect_conjunction,
detect_hill_sphere_entry, detect_opposition, detect_periapsis_crossing, is_in_shadow,
launch_window_phase, orbital_speed_circular,
};
use solarsystems::dynamics::integrator::{Integrator, adaptive_step};
use solarsystems::dynamics::stability::{
angular_momentum_drift_analysis, detect_resonance, energy_drift_analysis,
full_stability_analysis, lyapunov_exponent, megno_indicator, propagate_centuries,
};
use solarsystems::dynamics::time::{
DAYS_PER_JULIAN_CENTURY, EphemerisTime, J2000_EPOCH_JD, SECONDS_PER_DAY, SECONDS_PER_HOUR,
SECONDS_PER_JULIAN_YEAR, SECONDS_PER_MINUTE, SECONDS_PER_SIDEREAL_YEAR,
calendar_to_julian_date, julian_date_to_calendar,
};
use solarsystems::exports::state::{
extract_3d, extract_5d, extract_6d, extract_7d, extract_8d, extract_by_id, extract_full,
extract_full_flat, extract_pairwise_distances, extract_positions_flat,
extract_state_vectors_flat, extract_trajectory, relative_state,
};
use solarsystems::gravity::interactions::{
gravitational_acceleration, j2_acceleration, mutual_gravitational_force,
relativistic_correction, tidal_acceleration, tidal_force_magnitude, total_acceleration,
};
use solarsystems::gravity::nbody::{
barycenter, compute_accelerations, compute_accelerations_relativistic,
compute_pure_gravitational, total_system_angular_momentum, total_system_energy,
};
use solarsystems::orbits::keplerian::{
OrbitalElements, hohmann_transfer_dv, solve_kepler, synodic_period,
};
use solarsystems::orbits::perturbations::{
atmospheric_drag_acceleration, j2_secular_precession_argp, j2_secular_precession_raan,
kozai_lidov_max_eccentricity, lense_thirring_precession, mercury_gr_precession_rate,
poynting_robertson_drag, solar_radiation_pressure_accel, third_body_perturbation,
tidal_orbital_decay_rate, yarkovsky_acceleration,
};
use solarsystems::orbits::transfers::{
bielliptic_transfer, burn_time, characteristic_energy, circular_velocity,
combined_plane_change_dv, delta_v_budget, escape_velocity, flight_path_angle,
gravity_assist_delta_v, gravity_assist_max_deflection, hohmann_transfer,
hohmann_vs_bielliptic_threshold, hyperbolic_excess_velocity, lambert_tof_estimate,
launch_window_next, launch_window_phase_angle, oberth_effect_dv, plane_change_dv,
tsiolkovsky_mass_ratio, v_infinity_from_c3, vis_viva,
};
use solarsystems::orchestrator::pipeline::Pipeline;
use solarsystems::orchestrator::scheduler::Scheduler;
use solarsystems::{BodyId, CelestialBody, Vec3};
struct SolarSystemDiagnostics {
step: usize,
elapsed_days: f64,
elapsed_years: f64,
total_energy: f64,
initial_energy: f64,
angular_momentum: Vec3,
angular_momentum_mag: f64,
initial_am_mag: f64,
energy_drift: f64,
am_drift: f64,
event_count: usize,
body_count: usize,
barycenter_pos: Vec3,
barycenter_vel: Vec3,
conjunction_count: usize,
opposition_count: usize,
periapsis_count: usize,
close_approach_count: usize,
mercury_precession_rate: f64,
earth_j2_raan_rate: f64,
earth_j2_argp_rate: f64,
earth_sun_force: f64,
moon_tidal_force: f64,
jupiter_tidal_on_earth: Vec3,
mercury_j2_accel: Vec3,
mercury_gr_correction: Vec3,
sun_grav_on_earth: Vec3,
hohmann_earth_mars: f64,
bielliptic_earth_jupiter: f64,
earth_mars_synodic: f64,
earth_mars_launch_phase: f64,
earth_mars_tof: f64,
jupiter_gravity_assist_dv: f64,
jupiter_max_deflection: f64,
plane_change_5deg: f64,
oberth_at_sun: f64,
leo_circular_v: f64,
earth_escape_v: f64,
vis_viva_mars_orbit: f64,
mass_ratio_delta_v_10: f64,
launch_window_next_days: f64,
mercury_kepler_solution: f64,
earth_orbital_elements: OrbitalElements,
mercury_orbital_period: f64,
kozai_max_ecc_60deg: f64,
srp_accel_at_1au: f64,
earth_surface_gravity: f64,
jupiter_escape_velocity: f64,
sun_schwarzschild: f64,
earth_hill_sphere: f64,
earth_roche_limit: f64,
earth_soi: f64,
mars_orbital_period: f64,
full_state_count: usize,
earth_full_state_ke: f64,
earth_full_state_dist: f64,
calendar_year: i32,
calendar_month: u32,
calendar_day: u32,
earth_in_shadow: bool,
earth_mars_angular_sep: f64,
earth_mars_conjunction: bool,
mars_opposition: bool,
earth_mars_launch_window_phase: f64,
jupiter_saturn_resonance: Option<(u32, u32)>,
stability_max_energy_drift: f64,
stability_lyapunov: f64,
stability_megno: f64,
stability_is_stable: bool,
trajectory_point_count: usize,
trajectory_last_pos_mag: f64,
trajectory_last_vel_mag: f64,
centuries_since_j2000: f64,
ephemeris_jd: f64,
default_bodies_count: usize,
softening_val: f64,
j2_constants_sum: f64,
time_constants_sum: f64,
periapsis_crossing: Option<f64>,
apoapsis_crossing: Option<f64>,
close_approach_earth_mars: bool,
hill_entry_moon_mars: bool,
orbital_speed_earth: f64,
adaptive_step_dt: f64,
energy_drift_max: f64,
energy_drift_mean: f64,
am_drift_max: f64,
am_drift_mean: f64,
lyapunov: f64,
megno: f64,
calendar_roundtrip_jd: f64,
eph_from_jd_centuries: f64,
extract_3d_count: usize,
extract_5d_sum: f64,
extract_6d_sum: f64,
extract_7d_sum: f64,
extract_8d_sum: f64,
positions_flat_len: usize,
state_vectors_flat_len: usize,
full_flat_len: usize,
pairwise_min_dist: f64,
relative_moon_earth_pos_mag: f64,
total_accel_earth_mag: f64,
nbody_system_energy: f64,
nbody_system_am_mag: f64,
earth_oe_periapsis: f64,
earth_oe_apoapsis: f64,
earth_oe_semi_latus_rectum: f64,
earth_oe_mean_motion: f64,
earth_oe_eccentric_anomaly: f64,
earth_oe_mean_anomaly: f64,
earth_oe_radius_at_nu: f64,
earth_oe_velocity_at_r: f64,
roundtrip_state_pos_mag: f64,
advanced_oe_true_anomaly: f64,
hohmann_dv_keplerian: f64,
true_anomaly_from_mean_val: f64,
from_body_state_sma: f64,
atmospheric_drag: f64,
third_body_pert_mag: f64,
yarkovsky: f64,
poynting_robertson: f64,
lense_thirring: f64,
tidal_decay: f64,
bielliptic_threshold: f64,
c3_energy: f64,
v_inf_from_c3: f64,
hyperbolic_excess: f64,
flight_path: f64,
burn_time_val: f64,
scheduler_elapsed_years: f64,
scheduler_event_count: usize,
pipeline_run_energy: f64,
vec3_mag_sq: f64,
vec3_dot: f64,
vec3_cross_mag: f64,
vec3_angle: f64,
vec3_rotated_mag: f64,
vec3_lerp_mag: f64,
earth_mean_density: f64,
earth_potential_energy: f64,
earth_angular_momentum_mag: f64,
earth_spec_am_mag: f64,
earth_spec_orbital_energy: f64,
earth_vis_viva_speed: f64,
body_id_counts: f64,
body_id_parent_sum: f64,
body_id_type_count: f64,
prev_rdot_earth: f64,
checksum: f64,
}
impl SolarSystemDiagnostics {
fn recompute_checksum(&mut self) {
self.checksum = self.step as f64
+ self.elapsed_days * 1e-3
+ self.elapsed_years * 1e-1
+ self.total_energy * 1e-40
+ self.angular_momentum_mag * 1e-45
+ self.energy_drift
+ self.am_drift
+ self.event_count as f64
+ self.body_count as f64
+ self.barycenter_pos.magnitude() * 1e-15
+ self.barycenter_vel.magnitude() * 1e-5
+ self.conjunction_count as f64
+ self.opposition_count as f64
+ self.periapsis_count as f64
+ self.close_approach_count as f64
+ self.mercury_precession_rate
+ self.earth_j2_raan_rate
+ self.earth_j2_argp_rate
+ self.earth_sun_force * 1e-25
+ self.moon_tidal_force * 1e-10
+ self.jupiter_tidal_on_earth.magnitude() * 1e10
+ self.mercury_j2_accel.magnitude() * 1e8
+ self.mercury_gr_correction.magnitude() * 1e8
+ self.sun_grav_on_earth.magnitude() * 1e6
+ self.hohmann_earth_mars * 1e-4
+ self.bielliptic_earth_jupiter * 1e-4
+ self.earth_mars_synodic * 1e-8
+ self.earth_mars_launch_phase
+ self.earth_mars_tof * 1e-8
+ self.jupiter_gravity_assist_dv * 1e-4
+ self.jupiter_max_deflection
+ self.plane_change_5deg * 1e-4
+ self.oberth_at_sun * 1e-4
+ self.leo_circular_v * 1e-4
+ self.earth_escape_v * 1e-4
+ self.vis_viva_mars_orbit * 1e-4
+ self.mass_ratio_delta_v_10
+ self.launch_window_next_days * 1e-7
+ self.mercury_kepler_solution
+ self.earth_orbital_elements.semi_major_axis * 1e-12
+ self.earth_orbital_elements.eccentricity
+ self.mercury_orbital_period * 1e-8
+ self.kozai_max_ecc_60deg
+ self.srp_accel_at_1au * 1e6
+ self.earth_surface_gravity
+ self.jupiter_escape_velocity * 1e-5
+ self.sun_schwarzschild * 1e-4
+ self.earth_hill_sphere * 1e-10
+ self.earth_roche_limit * 1e-7
+ self.earth_soi * 1e-10
+ self.mars_orbital_period * 1e-8
+ self.full_state_count as f64
+ self.earth_full_state_ke * 1e-36
+ self.earth_full_state_dist * 1e-12
+ self.calendar_year as f64
+ self.calendar_month as f64
+ self.calendar_day as f64
+ if self.earth_in_shadow { 1.0 } else { 0.0 }
+ self.earth_mars_angular_sep
+ if self.earth_mars_conjunction {
1.0
} else {
0.0
}
+ if self.mars_opposition { 1.0 } else { 0.0 }
+ self.earth_mars_launch_window_phase
+ self
.jupiter_saturn_resonance
.map_or(0.0, |(p, q)| p as f64 + q as f64 * 0.1)
+ self.stability_max_energy_drift * 1e10
+ self.stability_lyapunov
+ self.stability_megno
+ if self.stability_is_stable { 1.0 } else { 0.0 }
+ self.trajectory_point_count as f64
+ self.trajectory_last_pos_mag * 1e-12
+ self.trajectory_last_vel_mag * 1e-5
+ self.centuries_since_j2000
+ self.ephemeris_jd * 1e-7
+ self.default_bodies_count as f64
+ self.softening_val * 1e-3
+ self.j2_constants_sum
+ self.time_constants_sum * 1e-8
+ self.periapsis_crossing.unwrap_or(0.0) * 1e-12
+ self.apoapsis_crossing.unwrap_or(0.0) * 1e-12
+ if self.close_approach_earth_mars {
1.0
} else {
0.0
}
+ if self.hill_entry_moon_mars { 1.0 } else { 0.0 }
+ self.orbital_speed_earth * 1e-5
+ self.adaptive_step_dt * 1e-5
+ self.energy_drift_max * 1e10
+ self.energy_drift_mean * 1e10
+ self.am_drift_max * 1e10
+ self.am_drift_mean * 1e10
+ self.lyapunov
+ self.megno
+ self.calendar_roundtrip_jd * 1e-7
+ self.eph_from_jd_centuries
+ self.extract_3d_count as f64
+ self.extract_5d_sum * 1e-15
+ self.extract_6d_sum * 1e-15
+ self.extract_7d_sum * 1e-15
+ self.extract_8d_sum * 1e-15
+ self.positions_flat_len as f64
+ self.state_vectors_flat_len as f64
+ self.full_flat_len as f64
+ self.pairwise_min_dist * 1e-12
+ self.relative_moon_earth_pos_mag * 1e-10
+ self.total_accel_earth_mag * 1e5
+ self.nbody_system_energy * 1e-40
+ self.nbody_system_am_mag * 1e-45
+ self.earth_oe_periapsis * 1e-12
+ self.earth_oe_apoapsis * 1e-12
+ self.earth_oe_semi_latus_rectum * 1e-12
+ self.earth_oe_mean_motion * 1e8
+ self.earth_oe_eccentric_anomaly
+ self.earth_oe_mean_anomaly
+ self.earth_oe_radius_at_nu * 1e-12
+ self.earth_oe_velocity_at_r * 1e-5
+ self.roundtrip_state_pos_mag * 1e-12
+ self.advanced_oe_true_anomaly
+ self.hohmann_dv_keplerian * 1e-4
+ self.true_anomaly_from_mean_val
+ self.from_body_state_sma * 1e-12
+ self.atmospheric_drag * 1e10
+ self.third_body_pert_mag * 1e5
+ self.yarkovsky * 1e15
+ self.poynting_robertson * 1e15
+ self.lense_thirring * 1e10
+ self.tidal_decay * 1e20
+ self.bielliptic_threshold * 1e-12
+ self.c3_energy * 1e-8
+ self.v_inf_from_c3 * 1e-4
+ self.hyperbolic_excess * 1e-4
+ self.flight_path
+ self.burn_time_val * 1e-3
+ self.scheduler_elapsed_years * 1e-3
+ self.scheduler_event_count as f64
+ self.pipeline_run_energy * 1e-40
+ self.vec3_mag_sq * 1e-25
+ self.vec3_dot * 1e-25
+ self.vec3_cross_mag * 1e-25
+ self.vec3_angle
+ self.vec3_rotated_mag * 1e-12
+ self.vec3_lerp_mag * 1e-12
+ self.earth_mean_density * 1e-4
+ self.earth_potential_energy * 1e-35
+ self.earth_angular_momentum_mag * 1e-40
+ self.earth_spec_am_mag * 1e-15
+ self.earth_spec_orbital_energy * 1e-10
+ self.earth_vis_viva_speed * 1e-5
+ self.body_id_counts
+ self.body_id_parent_sum
+ self.body_id_type_count
+ self.prev_rdot_earth * 1e-15;
}
}
fn count_events(
events: &[solarsystems::dynamics::events::CelestialEvent],
) -> (usize, usize, usize, usize) {
use solarsystems::dynamics::events::CelestialEvent;
let mut conj = 0usize;
let mut opp = 0usize;
let mut peri = 0usize;
let mut close = 0usize;
for e in events {
match e {
CelestialEvent::Conjunction { .. } => conj += 1,
CelestialEvent::Opposition { .. } => opp += 1,
CelestialEvent::Periapsis { .. } => peri += 1,
CelestialEvent::Apoapsis { .. } => peri += 1,
CelestialEvent::CloseApproach { .. } => close += 1,
_ => {}
}
}
(conj, opp, peri, close)
}
fn compute_initial_metrics(bodies: &[CelestialBody]) -> SolarSystemDiagnostics {
let sun = find_body(bodies, BodyId::Sun).unwrap();
let mercury = find_body(bodies, BodyId::Mercury).unwrap();
let earth = find_body(bodies, BodyId::Earth).unwrap();
let mars = find_body(bodies, BodyId::Mars).unwrap();
let jupiter = find_body(bodies, BodyId::Jupiter).unwrap();
let saturn = find_body(bodies, BodyId::Saturn).unwrap();
let moon = find_body(bodies, BodyId::Moon).unwrap();
let mercury_a = mercury.position.distance(&sun.position);
let earth_a = earth.position.distance(&sun.position);
let mars_a = mars.position.distance(&sun.position);
let jupiter_a = jupiter.position.distance(&sun.position);
let moon_dist = moon.position.distance(&earth.position);
let mercury_e = 0.2056;
let earth_e = 0.0167;
let earth_inc = 0.0_f64.to_radians();
let softening_val = GRAVITATIONAL_SOFTENING;
let j2_constants_sum = J2_EARTH + J2_JUPITER + J2_SATURN + J2_MARS + J2_SUN;
let default_bodies_count = default_bodies().len();
let time_constants_sum = SECONDS_PER_MINUTE
+ SECONDS_PER_HOUR
+ SECONDS_PER_DAY
+ SECONDS_PER_JULIAN_YEAR
+ SECONDS_PER_SIDEREAL_YEAR
+ DAYS_PER_JULIAN_CENTURY
+ J2000_EPOCH_JD;
let roundtrip_jd = calendar_to_julian_date(2000, 1, 1, 12, 0, 0.0);
let eph_from_jd = EphemerisTime::from_jd(roundtrip_jd);
let eph_from_jd_centuries = eph_from_jd.centuries_since_j2000();
let mercury_precession_rate =
mercury_gr_precession_rate(mercury_a, mercury_e, MU_SUN) + J2_SUN * 1e-20;
let earth_j2_raan_rate =
j2_secular_precession_raan(J2_EARTH, earth.radius, earth_a, earth_e, earth_inc, MU_SUN);
let earth_j2_argp_rate =
j2_secular_precession_argp(J2_EARTH, earth.radius, earth_a, earth_e, earth_inc, MU_SUN);
let atmospheric_drag = atmospheric_drag_acceleration(1e-12, 29780.0, 2.2, 0.01);
let third_body_pert = third_body_perturbation(earth.position, jupiter.position, jupiter.mass);
let solar_luminosity = solarsystems::stars::suns::luminosity();
let yarkovsky_val = yarkovsky_acceleration(0.01, 0.1, 500.0, 2000.0, AU, solar_luminosity);
let pr_drag = poynting_robertson_drag(solar_luminosity, 0.001, 3000.0, AU, 100.0);
let lense_thirring_val = lense_thirring_precession(1.9e41, earth_a, MU_SUN);
let tidal_decay =
tidal_orbital_decay_rate(moon.mass, earth.mass, earth.radius, moon_dist, 12.0, 0.3);
let kozai_max_ecc_60deg = kozai_lidov_max_eccentricity(60.0_f64.to_radians());
let srp_accel_at_1au = solar_radiation_pressure_accel(solar_luminosity, 0.01, 1.0, AU);
let earth_sun_force = mutual_gravitational_force(earth, sun);
let moon_tidal_force = tidal_force_magnitude(moon.mass, moon_dist, earth.radius);
let jupiter_tidal_on_earth =
tidal_acceleration(earth.position, jupiter.mass, jupiter.position, sun.position);
let mercury_j2_accel = j2_acceleration(mercury, sun);
let mercury_gr_correction = relativistic_correction(mercury, sun);
let sun_grav_on_earth = gravitational_acceleration(earth, sun);
let total_accel_earth = total_acceleration(earth, bodies);
let nbody_system_energy = total_system_energy(bodies);
let nbody_system_am = total_system_angular_momentum(bodies);
let hohmann = hohmann_transfer(earth_a, mars_a, MU_SUN);
let bielliptic = bielliptic_transfer(earth_a, jupiter_a, jupiter_a * 2.0, MU_SUN);
let earth_mars_synodic = synodic_period(
earth.orbital_period(sun.mass),
mars.orbital_period(sun.mass),
);
let earth_mars_launch_phase = launch_window_phase_angle(earth_a, mars_a, MU_SUN);
let earth_mars_tof = lambert_tof_estimate(
earth.position - sun.position,
mars.position - sun.position,
MU_SUN,
true,
);
let v_inf_jupiter = 5500.0;
let jupiter_periapsis = jupiter.radius + 200_000.0;
let jupiter_gravity_assist_dv =
gravity_assist_delta_v(v_inf_jupiter, jupiter.mass, jupiter_periapsis);
let jupiter_max_deflection =
gravity_assist_max_deflection(v_inf_jupiter, jupiter.mass, jupiter_periapsis);
let plane_change_5deg = plane_change_dv(29_780.0, 5.0_f64.to_radians());
let v_circ_leo = circular_velocity(earth.gravitational_parameter(), earth.radius + 400_000.0);
let oberth_at_sun = oberth_effect_dv(617_000.0, 1000.0);
let leo_circular_v =
circular_velocity(earth.gravitational_parameter(), earth.radius + 400_000.0);
let earth_escape_v = escape_velocity(earth.gravitational_parameter(), earth.radius + 400_000.0);
let vis_viva_mars_orbit = vis_viva(MU_SUN, mars_a, mars_a);
let combined_dv = combined_plane_change_dv(v_circ_leo, earth_escape_v, 28.5_f64.to_radians());
let mass_ratio_delta_v_10 = tsiolkovsky_mass_ratio(10_000.0, 350.0);
let budget = delta_v_budget(&[hohmann.total_delta_v, combined_dv]);
let launch_window_next_days =
launch_window_next(0.0, earth_mars_launch_phase, earth_mars_synodic);
let bielliptic_threshold = hohmann_vs_bielliptic_threshold(earth_a, MU_SUN);
let c3_energy = characteristic_energy(v_inf_jupiter);
let v_inf_from_c3_val = v_infinity_from_c3(c3_energy);
let v_departure = (v_circ_leo * v_circ_leo + v_inf_jupiter * v_inf_jupiter).sqrt();
let hyperbolic_excess_val = hyperbolic_excess_velocity(v_departure, v_circ_leo);
let earth_h = earth.specific_angular_momentum().magnitude();
let flight_path_val = flight_path_angle(earth_a, earth.velocity.magnitude(), earth_h);
let burn_time_val = burn_time(hohmann.delta_v1, 200_000.0, 500_000.0);
let mercury_kepler_solution = solve_kepler(1.5, 0.2056);
let earth_oe = OrbitalElements::from_state_vectors(
earth.position - sun.position,
earth.velocity - sun.velocity,
MU_SUN,
);
let mercury_orbital_period = mercury.orbital_period(sun.mass);
let earth_oe_periapsis = earth_oe.periapsis();
let earth_oe_apoapsis = earth_oe.apoapsis();
let earth_oe_slr = earth_oe.semi_latus_rectum();
let earth_oe_mean_motion = earth_oe.mean_motion(MU_SUN);
let earth_oe_ecc_anomaly = earth_oe.eccentric_anomaly();
let earth_oe_mean_anomaly = earth_oe.mean_anomaly();
let earth_oe_radius_at_nu = earth_oe.radius_at_true_anomaly(earth_oe.true_anomaly);
let earth_oe_vel_at_r = earth_oe.velocity_at_radius(earth_a, MU_SUN);
let (rt_pos, rt_vel) = earth_oe.to_state_vectors(MU_SUN);
let roundtrip_state_pos_mag = rt_pos.magnitude() + rt_vel.magnitude() * 1e-10;
let mut advanced_oe = earth_oe.clone();
advanced_oe.advance(86400.0, MU_SUN);
let advanced_oe_true_anomaly = advanced_oe.true_anomaly;
let (hdv1, hdv2) = hohmann_transfer_dv(earth_a, mars_a, MU_SUN);
let hohmann_dv_keplerian = hdv1 + hdv2;
let ta_from_mean = OrbitalElements::true_anomaly_from_mean(1.5, 0.2);
let from_body_oe = OrbitalElements::from_body_state(
earth.position - sun.position,
earth.velocity - sun.velocity,
sun.mass,
);
let earth_surface_gravity = earth.surface_gravity();
let jupiter_escape_velocity = jupiter.escape_velocity();
let sun_schwarzschild = sun.schwarzschild_radius();
let earth_hill_sphere = earth.hill_sphere(sun.mass, earth_a);
let earth_roche_limit = earth.roche_limit(3344.0);
let earth_soi = earth.sphere_of_influence(sun.mass, earth_a);
let mars_orbital_period = mars.orbital_period(sun.mass);
let earth_mean_density = earth.mean_density();
let earth_pe = earth.potential_energy(sun);
let earth_am = earth.angular_momentum();
let earth_spec_am = earth.specific_angular_momentum();
let earth_spec_oe = earth.specific_orbital_energy(sun.mass);
let earth_vv_speed = earth.vis_viva_speed(sun.mass, earth_a);
let earth_pos = earth.position;
let mars_pos = mars.position;
let vec3_mag_sq = earth_pos.magnitude_squared();
let vec3_normalized = earth_pos.normalize();
let vec3_dot = earth_pos.dot(&mars_pos);
let vec3_cross = earth_pos.cross(&mars_pos);
let vec3_angle = earth_pos.angle_between(&mars_pos);
let vec3_rotated = earth_pos.rotate_x(0.1).rotate_y(0.2).rotate_z(0.3);
let vec3_lerp = earth_pos.lerp(&mars_pos, 0.5);
let all_planets = BodyId::all_planets();
let all_sats = BodyId::all_satellites();
let all_dwarfs = BodyId::all_dwarf_planets();
let all_asteroids = BodyId::all_asteroids();
let all_comets = BodyId::all_comets();
let all_bodies = BodyId::all();
let body_id_counts = all_planets.len() as f64
+ all_sats.len() as f64
+ all_dwarfs.len() as f64
+ all_asteroids.len() as f64
+ all_comets.len() as f64
+ all_bodies.len() as f64;
let mut parent_sum = 0.0_f64;
for bid in all_sats {
if let Some(p) = bid.parent() {
parent_sum += if p.is_planet() { 1.0 } else { 0.0 };
}
}
let mut type_count = 0.0_f64;
for bid in all_bodies {
type_count += if bid.is_planet() { 1.0 } else { 0.0 };
type_count += if bid.is_satellite() { 0.1 } else { 0.0 };
type_count += if bid.is_dwarf_planet() { 0.01 } else { 0.0 };
type_count += if bid.is_asteroid() { 0.001 } else { 0.0 };
type_count += if bid.is_comet() { 0.0001 } else { 0.0 };
}
let states_3d = extract_3d(bodies);
let states_5d = extract_5d(bodies);
let states_6d = extract_6d(bodies);
let states_7d = extract_7d(bodies);
let states_8d = extract_8d(bodies);
let full_states = extract_full(bodies);
let earth_full = extract_by_id(bodies, BodyId::Earth).unwrap();
let positions_flat = extract_positions_flat(bodies);
let state_vectors_flat = extract_state_vectors_flat(bodies);
let full_flat = extract_full_flat(bodies);
let pairwise = extract_pairwise_distances(bodies);
let pairwise_min = pairwise
.iter()
.map(|(_, _, d)| *d)
.fold(f64::INFINITY, f64::min);
let rel_moon_earth = relative_state(bodies, BodyId::Moon, BodyId::Earth).unwrap();
let extract_5d_sum: f64 = states_5d
.iter()
.take(10)
.map(|s| s.speed * 1e-5 + s.mass * 1e-30)
.sum();
let extract_6d_sum: f64 = states_6d
.iter()
.take(10)
.map(|s| s.velocity.magnitude() * 1e-5)
.sum();
let extract_7d_sum: f64 = states_7d.iter().take(10).map(|s| s.mass * 1e-30).sum();
let extract_8d_sum: f64 = states_8d.iter().take(10).map(|s| s.radius * 1e-8).sum();
let earth_vel_rel = earth.velocity - sun.velocity;
let earth_pos_rel = earth.position - sun.position;
let prev_rdot_earth = earth_pos_rel.dot(&earth_vel_rel);
let periapsis_crossing = detect_periapsis_crossing(earth, sun, prev_rdot_earth);
let apoapsis_crossing = detect_apoapsis_crossing(earth, sun, prev_rdot_earth);
let close_approach_earth_mars = detect_close_approach(earth, mars, 1e12);
let hill_entry = detect_hill_sphere_entry(earth, moon, sun.mass, earth_a);
let orbital_speed_earth = orbital_speed_circular(sun.mass, earth_a);
let earth_in_shadow = is_in_shadow(earth, sun, moon);
let earth_mars_angular_sep = angular_separation(earth.position, mars.position, sun.position);
let earth_mars_conjunction = detect_conjunction(sun, earth, mars, 5.0_f64.to_radians());
let mars_opposition = detect_opposition(sun, earth, mars);
let earth_mars_launch_phase_dynamic = launch_window_phase(earth, mars, sun);
let jupiter_period = jupiter.orbital_period(sun.mass);
let saturn_period = saturn.orbital_period(sun.mass);
let jupiter_saturn_resonance = detect_resonance(jupiter_period, saturn_period, 10);
let ref_snapshots =
propagate_centuries(bodies.to_vec(), DEFAULT_TIMESTEP_S, 0.01, 86400.0 * 30.0);
let mut perturbed = bodies.to_vec();
perturbed[0].position.x += 1.0;
let pert_snapshots = propagate_centuries(perturbed, DEFAULT_TIMESTEP_S, 0.01, 86400.0 * 30.0);
let stability = full_stability_analysis(&ref_snapshots, &pert_snapshots);
let (ed_max, ed_mean) = energy_drift_analysis(&ref_snapshots);
let (amd_max, amd_mean) = angular_momentum_drift_analysis(&ref_snapshots);
let lyapunov_val = lyapunov_exponent(&ref_snapshots, &pert_snapshots);
let megno_val = megno_indicator(&ref_snapshots, &pert_snapshots);
let earth_trajectory = extract_trajectory(&ref_snapshots, BodyId::Earth);
let traj_count = earth_trajectory.len();
let (traj_last_pos_mag, traj_last_vel_mag) = earth_trajectory
.last()
.map(|(t, p, v)| (p.magnitude() + t * 1e-15, v.magnitude()))
.unwrap_or((0.0, 0.0));
let eph = EphemerisTime::j2000();
let centuries_since = eph.centuries_since_j2000();
let eph_jd = eph.current_jd();
let mut adaptive_bodies = bodies[..5].to_vec();
let adaptive_dt = adaptive_step(
&mut adaptive_bodies,
DEFAULT_TIMESTEP_S,
1e-10,
&Integrator::VelocityVerlet,
);
let mut scheduler = Scheduler::with_integrator(DEFAULT_TIMESTEP_S, Integrator::Yoshida4);
scheduler.run(10);
let scheduler_elapsed_years = scheduler.elapsed_years();
let scheduler_event_count = scheduler.event_count();
let mut run_pipeline = Pipeline::new(
default_bodies(),
DEFAULT_TIMESTEP_S,
Integrator::VelocityVerlet,
);
run_pipeline.initialize();
run_pipeline.run(5);
let pipeline_run_energy = run_pipeline.total_energy();
let mut nbody_copy = bodies[..10].to_vec();
compute_accelerations(&mut nbody_copy);
let accel_sum_1: f64 = nbody_copy.iter().map(|b| b.acceleration.magnitude()).sum();
let mut nbody_copy_2 = bodies[..10].to_vec();
compute_accelerations_relativistic(&mut nbody_copy_2);
let accel_sum_2: f64 = nbody_copy_2
.iter()
.map(|b| b.acceleration.magnitude())
.sum();
let mut nbody_copy_3 = bodies[..10].to_vec();
compute_pure_gravitational(&mut nbody_copy_3);
let accel_sum_3: f64 = nbody_copy_3
.iter()
.map(|b| b.acceleration.magnitude())
.sum();
let (bary_pos, bary_vel) = system_barycenter(bodies);
let total_e = total_energy(bodies);
let total_am = total_angular_momentum(bodies);
SolarSystemDiagnostics {
step: 0,
elapsed_days: 0.0,
elapsed_years: 0.0,
total_energy: total_e,
initial_energy: total_e,
angular_momentum: total_am,
angular_momentum_mag: total_am.magnitude(),
initial_am_mag: total_am.magnitude(),
energy_drift: 0.0,
am_drift: 0.0,
event_count: 0,
body_count: bodies.len(),
barycenter_pos: bary_pos,
barycenter_vel: bary_vel,
conjunction_count: 0,
opposition_count: 0,
periapsis_count: 0,
close_approach_count: 0,
mercury_precession_rate,
earth_j2_raan_rate,
earth_j2_argp_rate,
earth_sun_force,
moon_tidal_force,
jupiter_tidal_on_earth,
mercury_j2_accel,
mercury_gr_correction,
sun_grav_on_earth,
hohmann_earth_mars: hohmann.total_delta_v,
bielliptic_earth_jupiter: bielliptic.total_delta_v,
earth_mars_synodic,
earth_mars_launch_phase,
earth_mars_tof,
jupiter_gravity_assist_dv,
jupiter_max_deflection,
plane_change_5deg,
oberth_at_sun,
leo_circular_v,
earth_escape_v,
vis_viva_mars_orbit,
mass_ratio_delta_v_10: mass_ratio_delta_v_10 + budget * 1e-10,
launch_window_next_days,
mercury_kepler_solution,
earth_orbital_elements: earth_oe,
mercury_orbital_period,
kozai_max_ecc_60deg,
srp_accel_at_1au,
earth_surface_gravity,
jupiter_escape_velocity,
sun_schwarzschild,
earth_hill_sphere,
earth_roche_limit,
earth_soi,
mars_orbital_period,
full_state_count: full_states.len(),
earth_full_state_ke: earth_full.kinetic_energy,
earth_full_state_dist: earth_full.distance_from_origin,
calendar_year: 2000,
calendar_month: 1,
calendar_day: 1,
earth_in_shadow,
earth_mars_angular_sep,
earth_mars_conjunction,
mars_opposition,
earth_mars_launch_window_phase: earth_mars_launch_phase_dynamic,
jupiter_saturn_resonance,
stability_max_energy_drift: stability.max_energy_drift,
stability_lyapunov: stability.lyapunov_exponent,
stability_megno: stability.megno,
stability_is_stable: stability.is_stable,
trajectory_point_count: traj_count,
trajectory_last_pos_mag: traj_last_pos_mag,
trajectory_last_vel_mag: traj_last_vel_mag,
centuries_since_j2000: centuries_since,
ephemeris_jd: eph_jd,
default_bodies_count,
softening_val,
j2_constants_sum,
time_constants_sum,
periapsis_crossing,
apoapsis_crossing,
close_approach_earth_mars,
hill_entry_moon_mars: hill_entry,
orbital_speed_earth,
adaptive_step_dt: adaptive_dt,
energy_drift_max: ed_max,
energy_drift_mean: ed_mean,
am_drift_max: amd_max,
am_drift_mean: amd_mean,
lyapunov: lyapunov_val,
megno: megno_val,
calendar_roundtrip_jd: roundtrip_jd,
eph_from_jd_centuries,
extract_3d_count: states_3d.len(),
extract_5d_sum,
extract_6d_sum,
extract_7d_sum,
extract_8d_sum,
positions_flat_len: positions_flat.len(),
state_vectors_flat_len: state_vectors_flat.len(),
full_flat_len: full_flat.len(),
pairwise_min_dist: pairwise_min,
relative_moon_earth_pos_mag: rel_moon_earth.position.magnitude(),
total_accel_earth_mag: total_accel_earth.magnitude(),
nbody_system_energy: nbody_system_energy
+ (accel_sum_1 + accel_sum_2 + accel_sum_3) * 1e-20,
nbody_system_am_mag: nbody_system_am.magnitude(),
earth_oe_periapsis,
earth_oe_apoapsis,
earth_oe_semi_latus_rectum: earth_oe_slr,
earth_oe_mean_motion,
earth_oe_eccentric_anomaly: earth_oe_ecc_anomaly,
earth_oe_mean_anomaly,
earth_oe_radius_at_nu,
earth_oe_velocity_at_r: earth_oe_vel_at_r,
roundtrip_state_pos_mag,
advanced_oe_true_anomaly,
hohmann_dv_keplerian,
true_anomaly_from_mean_val: ta_from_mean,
from_body_state_sma: from_body_oe.semi_major_axis + vec3_normalized.magnitude() * 1e-20,
atmospheric_drag,
third_body_pert_mag: third_body_pert.magnitude(),
yarkovsky: yarkovsky_val,
poynting_robertson: pr_drag,
lense_thirring: lense_thirring_val,
tidal_decay,
bielliptic_threshold,
c3_energy,
v_inf_from_c3: v_inf_from_c3_val,
hyperbolic_excess: hyperbolic_excess_val,
flight_path: flight_path_val,
burn_time_val,
scheduler_elapsed_years,
scheduler_event_count,
pipeline_run_energy,
vec3_mag_sq,
vec3_dot,
vec3_cross_mag: vec3_cross.magnitude(),
vec3_angle,
vec3_rotated_mag: vec3_rotated.magnitude(),
vec3_lerp_mag: vec3_lerp.magnitude(),
earth_mean_density,
earth_potential_energy: earth_pe,
earth_angular_momentum_mag: earth_am.magnitude(),
earth_spec_am_mag: earth_spec_am.magnitude(),
earth_spec_orbital_energy: earth_spec_oe,
earth_vis_viva_speed: earth_vv_speed,
body_id_counts,
body_id_parent_sum: parent_sum,
body_id_type_count: type_count,
prev_rdot_earth,
checksum: 0.0,
}
}
fn main() {
let bodies = default_bodies_with_belts(200, 200);
let mut pipeline = Pipeline::new(bodies, DEFAULT_TIMESTEP_S, Integrator::LeapfrogKDK);
pipeline.initialize();
let mut diagnostics = compute_initial_metrics(&pipeline.bodies);
diagnostics.initial_energy = diagnostics.total_energy;
diagnostics.initial_am_mag = diagnostics.angular_momentum_mag;
let mut step: usize = 0;
let target_dt = std::time::Duration::from_secs_f64(1.0 / 60.0);
let mut last = std::time::Instant::now();
loop {
let now = std::time::Instant::now();
let real_dt = now.duration_since(last).as_secs_f64();
last = now;
diagnostics.checksum += real_dt * 1e-15;
pipeline.step();
step += 1;
let (bary_pos, bary_vel) = barycenter(&pipeline.bodies);
let te = pipeline.total_energy();
let am = pipeline.total_angular_momentum();
let (conj, opp, peri, close) = count_events(&pipeline.events);
diagnostics.step = step;
diagnostics.elapsed_days = pipeline.time.days_elapsed();
diagnostics.elapsed_years = pipeline.time.years_elapsed();
diagnostics.total_energy = te;
diagnostics.angular_momentum = am;
diagnostics.angular_momentum_mag = am.magnitude();
diagnostics.energy_drift =
(te - diagnostics.initial_energy).abs() / diagnostics.initial_energy.abs();
diagnostics.am_drift =
(am.magnitude() - diagnostics.initial_am_mag).abs() / diagnostics.initial_am_mag.abs();
diagnostics.event_count = pipeline.events.len();
diagnostics.body_count = pipeline.bodies.len();
diagnostics.barycenter_pos = bary_pos;
diagnostics.barycenter_vel = bary_vel;
diagnostics.conjunction_count = conj;
diagnostics.opposition_count = opp;
diagnostics.periapsis_count = peri;
diagnostics.close_approach_count = close;
let jd = pipeline.time.current_jd();
let (year, month, day, hour, minute, second) = julian_date_to_calendar(jd);
diagnostics.checksum += (hour as f64 + minute as f64 + second) * 1e-15;
diagnostics.calendar_year = year;
diagnostics.calendar_month = month;
diagnostics.calendar_day = day;
{
let sun = find_body(&pipeline.bodies, BodyId::Sun).unwrap();
let earth = find_body(&pipeline.bodies, BodyId::Earth).unwrap();
diagnostics.periapsis_crossing =
detect_periapsis_crossing(earth, sun, diagnostics.prev_rdot_earth);
diagnostics.apoapsis_crossing =
detect_apoapsis_crossing(earth, sun, diagnostics.prev_rdot_earth);
let rel_pos = earth.position - sun.position;
let rel_vel = earth.velocity - sun.velocity;
diagnostics.prev_rdot_earth = rel_pos.dot(&rel_vel);
}
if step.is_multiple_of(1000) {
let sun = find_body(&pipeline.bodies, BodyId::Sun).unwrap();
let earth = find_body(&pipeline.bodies, BodyId::Earth).unwrap();
let mars = find_body(&pipeline.bodies, BodyId::Mars).unwrap();
let moon = find_body(&pipeline.bodies, BodyId::Moon).unwrap();
let jupiter = find_body(&pipeline.bodies, BodyId::Jupiter).unwrap();
diagnostics.sun_grav_on_earth = gravitational_acceleration(earth, sun);
diagnostics.mercury_j2_accel =
j2_acceleration(find_body(&pipeline.bodies, BodyId::Mercury).unwrap(), sun);
diagnostics.mercury_gr_correction =
relativistic_correction(find_body(&pipeline.bodies, BodyId::Mercury).unwrap(), sun);
diagnostics.jupiter_tidal_on_earth =
tidal_acceleration(earth.position, jupiter.mass, jupiter.position, sun.position);
diagnostics.earth_sun_force = mutual_gravitational_force(earth, sun);
diagnostics.moon_tidal_force = tidal_force_magnitude(
moon.mass,
moon.position.distance(&earth.position),
earth.radius,
);
diagnostics.total_accel_earth_mag =
total_acceleration(earth, &pipeline.bodies).magnitude();
diagnostics.earth_in_shadow = is_in_shadow(earth, sun, moon);
diagnostics.earth_mars_angular_sep =
angular_separation(earth.position, mars.position, sun.position);
diagnostics.earth_mars_conjunction =
detect_conjunction(sun, earth, mars, 5.0_f64.to_radians());
diagnostics.mars_opposition = detect_opposition(sun, earth, mars);
diagnostics.earth_mars_launch_window_phase = launch_window_phase(earth, mars, sun);
diagnostics.close_approach_earth_mars = detect_close_approach(earth, mars, 1e12);
let earth_a = earth.position.distance(&sun.position);
diagnostics.hill_entry_moon_mars =
detect_hill_sphere_entry(earth, moon, sun.mass, earth_a);
diagnostics.orbital_speed_earth = orbital_speed_circular(sun.mass, earth_a);
let oe = OrbitalElements::from_state_vectors(
earth.position - sun.position,
earth.velocity - sun.velocity,
MU_SUN,
);
diagnostics.earth_orbital_elements = oe.clone();
diagnostics.earth_oe_periapsis = oe.periapsis();
diagnostics.earth_oe_apoapsis = oe.apoapsis();
diagnostics.earth_oe_semi_latus_rectum = oe.semi_latus_rectum();
diagnostics.earth_oe_mean_motion = oe.mean_motion(MU_SUN);
diagnostics.earth_oe_eccentric_anomaly = oe.eccentric_anomaly();
diagnostics.earth_oe_mean_anomaly = oe.mean_anomaly();
diagnostics.earth_oe_radius_at_nu = oe.radius_at_true_anomaly(oe.true_anomaly);
diagnostics.earth_oe_velocity_at_r = oe.velocity_at_radius(earth_a, MU_SUN);
let (rt_p, rt_v) = oe.to_state_vectors(MU_SUN);
diagnostics.roundtrip_state_pos_mag = rt_p.magnitude() + rt_v.magnitude() * 1e-10;
let mut adv = oe.clone();
adv.advance(86400.0, MU_SUN);
diagnostics.advanced_oe_true_anomaly = adv.true_anomaly;
let earth_full = extract_by_id(&pipeline.bodies, BodyId::Earth).unwrap();
diagnostics.earth_full_state_ke = earth_full.kinetic_energy;
diagnostics.earth_full_state_dist = earth_full.distance_from_origin;
diagnostics.full_state_count = extract_full(&pipeline.bodies).len();
diagnostics.extract_3d_count = extract_3d(&pipeline.bodies).len();
diagnostics.positions_flat_len = extract_positions_flat(&pipeline.bodies).len();
diagnostics.state_vectors_flat_len = extract_state_vectors_flat(&pipeline.bodies).len();
diagnostics.full_flat_len = extract_full_flat(&pipeline.bodies).len();
let pw = extract_pairwise_distances(&pipeline.bodies);
diagnostics.pairwise_min_dist =
pw.iter().map(|(_, _, d)| *d).fold(f64::INFINITY, f64::min);
let rel = relative_state(&pipeline.bodies, BodyId::Moon, BodyId::Earth).unwrap();
diagnostics.relative_moon_earth_pos_mag = rel.position.magnitude();
let s5 = extract_5d(&pipeline.bodies);
diagnostics.extract_5d_sum = s5
.iter()
.take(10)
.map(|s| s.speed * 1e-5 + s.mass * 1e-30)
.sum();
let s6 = extract_6d(&pipeline.bodies);
diagnostics.extract_6d_sum = s6
.iter()
.take(10)
.map(|s| s.velocity.magnitude() * 1e-5)
.sum();
let s7 = extract_7d(&pipeline.bodies);
diagnostics.extract_7d_sum = s7.iter().take(10).map(|s| s.mass * 1e-30).sum();
let s8 = extract_8d(&pipeline.bodies);
diagnostics.extract_8d_sum = s8.iter().take(10).map(|s| s.radius * 1e-8).sum();
diagnostics.earth_mean_density = earth.mean_density();
diagnostics.earth_potential_energy = earth.potential_energy(sun);
diagnostics.earth_angular_momentum_mag = earth.angular_momentum().magnitude();
diagnostics.earth_spec_am_mag = earth.specific_angular_momentum().magnitude();
diagnostics.earth_spec_orbital_energy = earth.specific_orbital_energy(sun.mass);
diagnostics.earth_vis_viva_speed = earth.vis_viva_speed(sun.mass, earth_a);
let ep = earth.position;
let mp = mars.position;
diagnostics.vec3_mag_sq = ep.magnitude_squared();
diagnostics.vec3_dot = ep.dot(&mp);
diagnostics.vec3_cross_mag = ep.cross(&mp).magnitude();
diagnostics.vec3_angle = ep.angle_between(&mp);
diagnostics.vec3_rotated_mag = ep.rotate_x(0.1).rotate_y(0.2).rotate_z(0.3).magnitude();
diagnostics.vec3_lerp_mag = ep.lerp(&mp, 0.5).magnitude();
diagnostics.nbody_system_energy = total_system_energy(&pipeline.bodies);
diagnostics.nbody_system_am_mag =
total_system_angular_momentum(&pipeline.bodies).magnitude();
let mut eph_now = EphemerisTime::j2000();
eph_now.advance(pipeline.time.elapsed_s);
diagnostics.centuries_since_j2000 = eph_now.centuries_since_j2000();
diagnostics.ephemeris_jd = eph_now.current_jd();
let third_bp = third_body_perturbation(earth.position, jupiter.position, jupiter.mass);
diagnostics.third_body_pert_mag = third_bp.magnitude();
diagnostics.tidal_decay = tidal_orbital_decay_rate(
moon.mass,
earth.mass,
earth.radius,
moon.position.distance(&earth.position),
12.0,
0.3,
);
}
if step.is_multiple_of(10000) {
let mut bodies_copy = pipeline.bodies.clone();
if let Some(belt) = find_body_mut(&mut bodies_copy, BodyId::BeltAsteroid(0)) {
belt.position.x += 1e-20;
diagnostics.checksum += belt.position.x * 1e-30;
}
let mut small = pipeline.bodies[..8].to_vec();
compute_accelerations(&mut small);
let a1: f64 = small.iter().map(|b| b.acceleration.magnitude()).sum();
let mut small2 = pipeline.bodies[..8].to_vec();
compute_accelerations_relativistic(&mut small2);
let a2: f64 = small2.iter().map(|b| b.acceleration.magnitude()).sum();
let mut small3 = pipeline.bodies[..8].to_vec();
compute_pure_gravitational(&mut small3);
let a3: f64 = small3.iter().map(|b| b.acceleration.magnitude()).sum();
diagnostics.checksum += (a1 + a2 + a3) * 1e-20;
let mut adaptive_bodies = pipeline.bodies[..5].to_vec();
diagnostics.adaptive_step_dt = adaptive_step(
&mut adaptive_bodies,
DEFAULT_TIMESTEP_S,
1e-10,
&Integrator::VelocityVerlet,
);
}
diagnostics.recompute_checksum();
let elapsed = now.elapsed();
if elapsed < target_dt {
std::thread::sleep(target_dt - elapsed);
}
}
}