From a223c883170c5a7ab823d847a0de747f15cc3ec1 Mon Sep 17 00:00:00 2001
From: Camden Dixie O'Brien <cdo@wip.sh>
Date: Sat, 17 Oct 2026 18:45:43 +0100
Subject: [PATCH] Create binary to output energy data for symplectic vs RK4

---
 src/bin/comp_rk4.rs | 41 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)
 create mode 100644 src/bin/comp_rk4.rs

diff --git a/src/bin/comp_rk4.rs b/src/bin/comp_rk4.rs
new file mode 100644
index 0000000..74e8541
--- /dev/null
+++ b/src/bin/comp_rk4.rs
@@ -0,0 +1,41 @@
+use symplectic_double_pendulum::*;
+
+use std::ops::{Add, Mul};
+
+fn rk4<T>(y: T, dt: f64, f: impl Fn(T) -> T) -> T
+where T: Add<Output = T> + Mul<f64, Output = T> + Copy
+{
+	let k1 = f(y);
+	let k2 = f(y + k1 * dt * 0.5);
+	let k3 = f(y + k2 * dt * 0.5);
+	let k4 = f(y + k3 * dt);
+	y + (k1 + k2 * 2.0 + k3 * 2.0 + k4) * (dt / 6.0)
+}
+
+fn rk4_evolve<T>(state: &mut T, h: f64)
+where T: Hamiltonian + Add<Output = T> + Mul<f64, Output = T> + Copy
+{
+	*state = rk4(*state, h, |s| s.dots());
+}
+
+fn main() {
+	let omega = 500.0;
+	let h = 1e-3;
+
+	let mut rk4 = State {
+		q: DoublePendulum { theta: std::f64::consts::PI - 0.01, phi: 0.0 },
+		p: DoublePendulum { theta: 0.0, phi: 0.0 },
+	};
+	let mut aug = Augmented::from(rk4, omega);
+
+	println!("time,rk4,tao");
+	for i in 0..36000000u32 {
+		rk4_evolve(&mut rk4, h);
+		aug.evolve(h);
+		if i % 1024 == 0 {
+			println!(
+				"{},{},{}",
+				h * (i as f64), rk4.energy(), aug.average().energy());
+		}
+	}
+}