From 0618b81763214906818c0c84e2bfcc6ee0e7c244 Mon Sep 17 00:00:00 2001 From: Verox001 Date: Sun, 4 May 2025 19:24:03 +0200 Subject: [PATCH] Improved orbit calculations and added n-body problem solution --- simulator/src/main.rs | 2 +- solar_engine/Cargo.toml | 1 + solar_engine/src/simulator.rs | 153 ++++++++++++++++++++++++++-------- 3 files changed, 120 insertions(+), 36 deletions(-) diff --git a/simulator/src/main.rs b/simulator/src/main.rs index 6f5dcbb..d4a10d0 100644 --- a/simulator/src/main.rs +++ b/simulator/src/main.rs @@ -287,7 +287,7 @@ impl<'a> State<'a> { index_count: circle_indices.len() as u32, }); - let mut sim = Simulator::new(60.0 * 60.0 * 24.0); + let mut sim = Simulator::new(43200.0); sim.add_body(Body { name: "Sun".to_string(), position: [0.0, 0.0], diff --git a/solar_engine/Cargo.toml b/solar_engine/Cargo.toml index 4877b9c..cca4314 100644 --- a/solar_engine/Cargo.toml +++ b/solar_engine/Cargo.toml @@ -4,3 +4,4 @@ version = "0.1.0" edition = "2021" [dependencies] +rayon = "1.8" \ No newline at end of file diff --git a/solar_engine/src/simulator.rs b/solar_engine/src/simulator.rs index 46cf029..9c85b6b 100644 --- a/solar_engine/src/simulator.rs +++ b/solar_engine/src/simulator.rs @@ -1,4 +1,6 @@ +use std::sync::Mutex; use crate::body::Body; +use rayon::prelude::*; const G: f64 = 6.67430e-11; @@ -28,48 +30,129 @@ impl Simulator { } pub fn step(&mut self) { - let n = self.bodies.len(); let dt = self.timestep; + let n = self.bodies.len(); - let mut accelerations = vec![[0.0, 0.0]; n]; - - for i in 0..n { - for j in 0..n { - if i == j { - continue; - } - - let (bi, bj) = (&self.bodies[i], &self.bodies[j]); - - let dx = bj.position[0] - bi.position[0]; - let dy = bj.position[1] - bi.position[1]; - let dist_sq = distance_squared(bi.position, bj.position); - let dist = dist_sq.sqrt(); - - if dist < 1e-3 { - continue; - } - - let force = G * bi.mass * bj.mass / dist_sq; - let accel = force / bi.mass; - - let ax = accel * dx / dist; - let ay = accel * dy / dist; - - accelerations[i][0] += ax; - accelerations[i][1] += ay; - } + #[derive(Clone)] + struct State { + position: [f64; 2], + velocity: [f64; 2], } - + + let original_states: Vec = self + .bodies + .iter() + .map(|b| State { + position: b.position, + velocity: b.velocity, + }) + .collect(); + + let masses: Vec = self.bodies.iter().map(|b| b.mass).collect(); + + fn compute_accelerations(states: &[State], masses: &[f64]) -> Vec<[f64; 2]> { + let n = states.len(); + let accels = (0..n).map(|_| Mutex::new([0.0, 0.0])).collect::>(); + + (0..n).into_par_iter().for_each(|i| { + for j in (i + 1)..n { + let dx = states[j].position[0] - states[i].position[0]; + let dy = states[j].position[1] - states[i].position[1]; + let dist_sq = dx * dx + dy * dy; + let dist = dist_sq.sqrt(); + + if dist < 1e-3 { + continue; + } + + let force = G * masses[i] * masses[j] / dist_sq; + let ax = force * dx / dist; + let ay = force * dy / dist; + + { + let mut a_i_lock = accels[i].lock().unwrap(); + a_i_lock[0] += ax / masses[i]; + a_i_lock[1] += ay / masses[i]; + } + { + let mut a_j_lock = accels[j].lock().unwrap(); + a_j_lock[0] -= ax / masses[j]; + a_j_lock[1] -= ay / masses[j]; + } + } + }); + + accels + .into_iter() + .map(|mutex| mutex.into_inner().unwrap()) + .collect() + } + + // RK4 Stufen + let k1_pos = original_states.iter().map(|s| s.velocity).collect::>(); + let k1_vel = compute_accelerations(&original_states, &masses); + + let mut temp_states = original_states + .iter() + .enumerate() + .map(|(i, s)| State { + position: [ + s.position[0] + k1_pos[i][0] * dt / 2.0, + s.position[1] + k1_pos[i][1] * dt / 2.0, + ], + velocity: [ + s.velocity[0] + k1_vel[i][0] * dt / 2.0, + s.velocity[1] + k1_vel[i][1] * dt / 2.0, + ], + }) + .collect::>(); + + let k2_pos = temp_states.iter().map(|s| s.velocity).collect::>(); + let k2_vel = compute_accelerations(&temp_states, &masses); + + for i in 0..n { + temp_states[i].position[0] = original_states[i].position[0] + k2_pos[i][0] * dt / 2.0; + temp_states[i].position[1] = original_states[i].position[1] + k2_pos[i][1] * dt / 2.0; + temp_states[i].velocity[0] = original_states[i].velocity[0] + k2_vel[i][0] * dt / 2.0; + temp_states[i].velocity[1] = original_states[i].velocity[1] + k2_vel[i][1] * dt / 2.0; + } + + let k3_pos = temp_states.iter().map(|s| s.velocity).collect::>(); + let k3_vel = compute_accelerations(&temp_states, &masses); + + for i in 0..n { + temp_states[i].position[0] = original_states[i].position[0] + k3_pos[i][0] * dt; + temp_states[i].position[1] = original_states[i].position[1] + k3_pos[i][1] * dt; + temp_states[i].velocity[0] = original_states[i].velocity[0] + k3_vel[i][0] * dt; + temp_states[i].velocity[1] = original_states[i].velocity[1] + k3_vel[i][1] * dt; + } + + let k4_pos = temp_states.iter().map(|s| s.velocity).collect::>(); + let k4_vel = compute_accelerations(&temp_states, &masses); + + // Finale Updates for i in 0..n { - let a = accelerations[i]; let body = &mut self.bodies[i]; + let old_position = body.position; - body.velocity[0] += a[0] * dt; - body.velocity[1] += a[1] * dt; + body.position[0] += (dt / 6.0) + * (k1_pos[i][0] + 2.0 * k2_pos[i][0] + 2.0 * k3_pos[i][0] + k4_pos[i][0]); + body.position[1] += (dt / 6.0) + * (k1_pos[i][1] + 2.0 * k2_pos[i][1] + 2.0 * k3_pos[i][1] + k4_pos[i][1]); - body.position[0] += body.velocity[0] * dt; - body.position[1] += body.velocity[1] * dt; + body.velocity[0] += (dt / 6.0) + * (k1_vel[i][0] + 2.0 * k2_vel[i][0] + 2.0 * k3_vel[i][0] + k4_vel[i][0]); + body.velocity[1] += (dt / 6.0) + * (k1_vel[i][1] + 2.0 * k2_vel[i][1] + 2.0 * k3_vel[i][1] + k4_vel[i][1]); + + // Bewegung loggen + let dx = body.position[0] - old_position[0]; + let dy = body.position[1] - old_position[1]; + let movement = (dx * dx + dy * dy).sqrt(); + + if i == 1 { + println!("Earth moved {:.6} meters", movement); + } } self.time += dt;