From c3c82ae06405fa2c7bacac1d03d148897321c32b Mon Sep 17 00:00:00 2001
From: eug-vs <eugene@eug-vs.xyz>
Date: Sun, 15 Dec 2024 23:00:51 +0100
Subject: feat: use Runge Kutta ODE solver

---
 physics/src/solver/rk4.rs | 37 +++++++++++++++++++++++++++++++++++++
 1 file changed, 37 insertions(+)
 create mode 100644 physics/src/solver/rk4.rs

(limited to 'physics/src/solver/rk4.rs')

diff --git a/physics/src/solver/rk4.rs b/physics/src/solver/rk4.rs
new file mode 100644
index 0000000..224bec5
--- /dev/null
+++ b/physics/src/solver/rk4.rs
@@ -0,0 +1,37 @@
+use crate::{algebra::Scalar, particle_system::ParticleSystem};
+use super::{PhaseSpace, Solver};
+
+// Runge-Kutta 4th order
+impl Solver for ParticleSystem {
+    fn step(&mut self, dt: Scalar) {
+        let state = self.collect_phase_space();
+
+        // k1 = f(t, y)
+        let k1 = self.compute_derivative();
+
+        // k2 = f(t + dt/2, y + dt * k1 / 2)
+        self.scatter_phase_space(&PhaseSpace {
+            0: state.0.clone() + k1.0.clone() * dt / 2.0,
+        });
+        let k2 = self.compute_derivative();
+
+
+        // k3 = f(t + dt/2, y + dt * k2 / 2)
+        self.scatter_phase_space(&PhaseSpace {
+            0: state.0.clone() + k2.0.clone() * dt / 2.0,
+        });
+        let k3 = self.compute_derivative();
+
+        // k4 = f(t + dt, y + dt * k3)
+        self.scatter_phase_space(&PhaseSpace {
+            0: state.0.clone() + k3.0.clone() * dt
+        });
+        let k4 = self.compute_derivative();
+
+        self.scatter_phase_space(&PhaseSpace {
+            0: state.0.clone() + (k1.0 + 2.0 * k2.0 + 2.0 * k3.0 + k4.0) * dt / 6.0
+        });
+
+        self.t += dt;
+    }
+}
-- 
cgit v1.2.3