Improved orbit calculations and added n-body problem solution

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],

solar_engine/Cargo.toml

@@ -4,3 +4,4 @@ version = "0.1.0"
 edition = "2021"
 
 [dependencies]
+rayon = "1.8"

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];
+        #[derive(Clone)]
-
+        struct State {
-        for i in 0..n {
+            position: [f64; 2],
-            for j in 0..n {
+            velocity: [f64; 2],
-                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;
-            }
         }
-        
+
+        let original_states: Vec<State> = self
+            .bodies
+            .iter()
+            .map(|b| State {
+                position: b.position,
+                velocity: b.velocity,
+            })
+            .collect();
+
+        let masses: Vec<f64> = 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::<Vec<_>>();
+
+            (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::<Vec<_>>();
+        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::<Vec<_>>();
+
+        let k2_pos = temp_states.iter().map(|s| s.velocity).collect::<Vec<_>>();
+        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::<Vec<_>>();
+        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::<Vec<_>>();
+        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.position[0] += (dt / 6.0)
-            body.velocity[1] += a[1] * dt;
+                * (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.velocity[0] += (dt / 6.0)
-            body.position[1] += body.velocity[1] * dt;
+                * (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;