use solarsystems::BodyId;
use solarsystems::config::parameters::{
DEFAULT_TIMESTEP_S, default_bodies, default_bodies_with_belts, find_body, system_barycenter,
total_angular_momentum, total_energy,
};
use solarsystems::dynamics::events::{
angular_separation, detect_close_approach, detect_conjunction, detect_opposition, is_in_shadow,
launch_window_phase, orbital_speed_circular, phase_angle_rate,
};
use solarsystems::dynamics::integrator::Integrator;
use solarsystems::dynamics::stability::{
detect_resonance, energy_drift_analysis, full_stability_analysis, propagate_centuries,
};
use solarsystems::dynamics::time::{
EphemerisTime, SECONDS_PER_DAY, SECONDS_PER_JULIAN_YEAR, 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::nbody::{
compute_accelerations, compute_accelerations_relativistic, compute_pure_gravitational,
total_system_angular_momentum, total_system_energy,
};
use solarsystems::orchestrator::pipeline::Pipeline;
use solarsystems::orchestrator::scheduler::Scheduler;
fn main() {
let bodies = default_bodies();
let planets = bodies.iter().filter(|b| b.id.is_planet()).count();
let satellites = bodies.iter().filter(|b| b.id.is_satellite()).count();
let dwarfs = bodies.iter().filter(|b| b.id.is_dwarf_planet()).count();
let asteroids = bodies.iter().filter(|b| b.id.is_asteroid()).count();
let comets = bodies.iter().filter(|b| b.id.is_comet()).count();
let belt_bodies = default_bodies_with_belts(50, 50);
let belt_asteroids = belt_bodies
.iter()
.filter(|b| b.id.is_belt_asteroid())
.count();
let kuiper_objects = belt_bodies
.iter()
.filter(|b| b.id.is_kuiper_object())
.count();
let belt_particles = belt_bodies
.iter()
.filter(|b| b.id.is_belt_particle())
.count();
let mut checksum: f64 = bodies.len() as f64
+ planets as f64
+ satellites as f64
+ dwarfs as f64
+ asteroids as f64
+ comets as f64
+ belt_bodies.len() as f64
+ belt_asteroids as f64
+ kuiper_objects as f64
+ belt_particles as f64;
let all_ids = BodyId::all();
let all_planet_ids = BodyId::all_planets();
let all_sat_ids = BodyId::all_satellites();
let all_dwarf_ids = BodyId::all_dwarf_planets();
let all_ast_ids = BodyId::all_asteroids();
let all_comet_ids = BodyId::all_comets();
checksum += all_ids.len() as f64
+ all_planet_ids.len() as f64
+ all_sat_ids.len() as f64
+ all_dwarf_ids.len() as f64
+ all_ast_ids.len() as f64
+ all_comet_ids.len() as f64;
for sid in all_sat_ids {
if let Some(p) = sid.parent() {
checksum += if p.is_planet() { 0.01 } else { 0.0 };
}
}
let sun = find_body(&bodies, BodyId::Sun).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();
checksum += earth.gravitational_parameter() * 1e-20
+ earth.surface_gravity()
+ earth.escape_velocity() * 1e-5
+ earth.mean_density() * 1e-4
+ earth.schwarzschild_radius() * 1e3
+ earth.hill_sphere(sun.mass, earth.position.distance(&sun.position)) * 1e-10
+ earth.roche_limit(3344.0) * 1e-8
+ earth.sphere_of_influence(sun.mass, earth.position.distance(&sun.position)) * 1e-10
+ earth.orbital_period(sun.mass) * 1e-8
+ earth.kinetic_energy() * 1e-35
+ earth.potential_energy(sun) * 1e-35
+ earth.angular_momentum().magnitude() * 1e-42
+ earth.specific_angular_momentum().magnitude() * 1e-16
+ earth.specific_orbital_energy(sun.mass) * 1e-10
+ earth.vis_viva_speed(sun.mass, earth.position.distance(&sun.position)) * 1e-5;
let (bary_pos, bary_vel) = system_barycenter(&bodies);
let te = total_energy(&bodies);
let tam = total_angular_momentum(&bodies);
checksum += bary_pos.magnitude() * 1e-15
+ bary_vel.magnitude() * 1e-5
+ te * 1e-40
+ tam.magnitude() * 1e-45;
let nbe = total_system_energy(&bodies);
let nbam = total_system_angular_momentum(&bodies);
checksum += nbe * 1e-40 + nbam.magnitude() * 1e-45;
let mut copy1 = bodies.clone();
compute_accelerations(&mut copy1);
let a1: f64 = copy1.iter().map(|b| b.acceleration.magnitude()).sum();
let mut copy2 = bodies.clone();
compute_accelerations_relativistic(&mut copy2);
let a2: f64 = copy2.iter().map(|b| b.acceleration.magnitude()).sum();
let mut copy3 = bodies.clone();
compute_pure_gravitational(&mut copy3);
let a3: f64 = copy3.iter().map(|b| b.acceleration.magnitude()).sum();
checksum += (a1 + a2 + a3) * 1e-10;
let earth_a = earth.position.distance(&sun.position);
let mars_a = mars.position.distance(&sun.position);
checksum += angular_separation(earth.position, mars.position, sun.position);
checksum += if detect_conjunction(sun, earth, mars, 5.0_f64.to_radians()) {
1.0
} else {
0.0
};
checksum += if detect_opposition(sun, earth, mars) {
1.0
} else {
0.0
};
checksum += if is_in_shadow(earth, sun, moon) {
1.0
} else {
0.0
};
checksum += if detect_close_approach(earth, mars, 1e12) {
1.0
} else {
0.0
};
checksum += launch_window_phase(earth, mars, sun);
checksum += orbital_speed_circular(sun.mass, earth_a) * 1e-5;
checksum += phase_angle_rate(earth_a, mars_a, sun.gravitational_parameter());
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 = extract_full(&bodies);
let earth_full = extract_by_id(&bodies, BodyId::Earth).unwrap();
let positions_flat = extract_positions_flat(&bodies);
let sv_flat = extract_state_vectors_flat(&bodies);
let full_flat = extract_full_flat(&bodies);
let pairwise = extract_pairwise_distances(&bodies);
let rel = relative_state(&bodies, BodyId::Moon, BodyId::Earth).unwrap();
checksum += states_3d.len() as f64
+ states_5d
.iter()
.take(5)
.map(|s| s.speed * 1e-5)
.sum::<f64>()
+ states_6d
.iter()
.take(5)
.map(|s| s.velocity.magnitude() * 1e-5)
.sum::<f64>()
+ states_7d
.iter()
.take(5)
.map(|s| s.mass * 1e-30)
.sum::<f64>()
+ states_8d
.iter()
.take(5)
.map(|s| s.radius * 1e-8)
.sum::<f64>()
+ full.len() as f64
+ earth_full.kinetic_energy * 1e-35
+ earth_full.distance_from_origin * 1e-12
+ positions_flat.len() as f64
+ sv_flat.len() as f64
+ full_flat.len() as f64
+ pairwise
.iter()
.map(|(_, _, d)| *d)
.fold(f64::INFINITY, f64::min)
* 1e-12
+ rel.position.magnitude() * 1e-10;
let jupiter_period = jupiter.orbital_period(sun.mass);
let saturn_period = saturn.orbital_period(sun.mass);
if let Some((p, q)) = detect_resonance(jupiter_period, saturn_period, 10) {
checksum += p as f64 + q as f64 * 0.1;
}
let dt = 3600.0;
let steps_per_year = 8766;
let years = 10;
let total = steps_per_year * years;
let mut scheduler = Scheduler::with_integrator(dt, Integrator::Yoshida4);
scheduler.output_interval_steps = steps_per_year;
let e0 = {
let mut s = Scheduler::default_system(dt);
s.pipeline.initialize();
s.pipeline.total_energy()
};
scheduler.run_with_callback(total, |_step, pipeline| {
let yr = pipeline.time.years_elapsed();
let days = yr * 365.25;
let e = pipeline.total_energy();
let drift = ((e - e0) / e0).abs();
checksum += drift * 1e10 + days * 1e-6;
let (y, m, d, h, min, sec) = julian_date_to_calendar(pipeline.time.current_jd());
checksum += (y as f64 + m as f64 + d as f64 + h as f64 + min as f64 + sec) * 1e-10;
});
checksum += scheduler.elapsed_years() * 1e-2
+ scheduler.elapsed_days() * 1e-5
+ scheduler.event_count() as f64;
let mut pipeline = Pipeline::new(
default_bodies(),
DEFAULT_TIMESTEP_S,
Integrator::VelocityVerlet,
);
pipeline.initialize();
pipeline.run(100);
let target_jd = pipeline.time.current_jd() + 1.0;
pipeline.run_until(target_jd);
let pam = pipeline.total_angular_momentum();
checksum += pipeline.total_energy() * 1e-40 + pam.magnitude() * 1e-45;
let snapshots = propagate_centuries(
default_bodies(),
DEFAULT_TIMESTEP_S,
0.01,
SECONDS_PER_DAY * 30.0,
);
let mut perturbed = default_bodies();
perturbed[0].position.x += 1.0;
let pert_snapshots =
propagate_centuries(perturbed, DEFAULT_TIMESTEP_S, 0.01, SECONDS_PER_DAY * 30.0);
let (ed_max, ed_mean) = energy_drift_analysis(&snapshots);
let stability = full_stability_analysis(&snapshots, &pert_snapshots);
checksum += ed_max
+ ed_mean
+ stability.max_energy_drift * 1e10
+ stability.lyapunov_exponent
+ stability.megno
+ if stability.is_stable { 1.0 } else { 0.0 };
let trajectory = extract_trajectory(&snapshots, BodyId::Earth);
checksum += trajectory.len() as f64;
if let Some((t, p, v)) = trajectory.last() {
checksum += t * 1e-15 + p.magnitude() * 1e-12 + v.magnitude() * 1e-5;
}
let eph = EphemerisTime::j2000();
checksum += eph.centuries_since_j2000()
+ eph.current_jd() * 1e-7
+ eph.days_elapsed() * 1e-10
+ eph.years_elapsed() * 1e-10
+ SECONDS_PER_DAY * 1e-10
+ SECONDS_PER_JULIAN_YEAR * 1e-10;
assert!(checksum.is_finite(), "checksum diverged");
assert!(checksum > 0.0, "checksum should be positive");
}