add leadfinder, add double pendulum

double_pendulum/.gitignore

@@ -0,0 +1,2 @@
+target
+out

double_pendulum/Cargo.lock

@@ -0,0 +1,61 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 4
+
+[[package]]
+name = "crossbeam-deque"
+version = "0.8.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9dd111b7b7f7d55b72c0a6ae361660ee5853c9af73f70c3c2ef6858b950e2e51"
+dependencies = [
+ "crossbeam-epoch",
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-epoch"
+version = "0.9.18"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e"
+dependencies = [
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-utils"
+version = "0.8.21"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d0a5c400df2834b80a4c3327b3aad3a4c4cd4de0629063962b03235697506a28"
+
+[[package]]
+name = "double_pendulum"
+version = "0.1.0"
+dependencies = [
+ "rayon",
+]
+
+[[package]]
+name = "either"
+version = "1.15.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719"
+
+[[package]]
+name = "rayon"
+version = "1.11.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "368f01d005bf8fd9b1206fb6fa653e6c4a81ceb1466406b81792d87c5677a58f"
+dependencies = [
+ "either",
+ "rayon-core",
+]
+
+[[package]]
+name = "rayon-core"
+version = "1.13.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91"
+dependencies = [
+ "crossbeam-deque",
+ "crossbeam-utils",
+]

double_pendulum/Cargo.toml

@@ -0,0 +1,7 @@
+[package]
+name = "double_pendulum"
+version = "0.1.0"
+edition = "2021"
+
+[dependencies]
+rayon = "1.10"

double_pendulum/src/main.rs

@@ -0,0 +1,244 @@
+use rayon::prelude::*;
+use std::f32::consts::PI;
+use std::io::Write;
+use std::process::{Command, Stdio};
+
+const WIDTH: usize = 2000;
+const HEIGHT: usize = 2000;
+
+const FPS: usize = 1;
+const DURATION: usize = 25;
+const TOTAL_FRAMES: usize = FPS * DURATION;
+
+const SUBSTEPS_PER_FRAME: usize = 100;
+
+const DT: f32 = 0.002;
+
+const EPS: f32 = 1e-6;
+
+const OUTPUT: &str = "double_pendulum_lyapunov.mp4";
+
+const G: f32 = 9.81;
+
+#[derive(Clone, Copy)]
+struct State {
+    th1: f32,
+    th2: f32,
+    w1: f32,
+    w2: f32,
+}
+
+#[derive(Clone, Copy)]
+struct Cell {
+    a: State,
+    b: State,
+    lyap: f32,
+}
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let mut cells = initialize();
+
+    let mut frame = vec![0u8; WIDTH * HEIGHT * 3];
+
+    let mut ffmpeg = Command::new("ffmpeg")
+        .args([
+            "-y",
+            "-f",
+            "rawvideo",
+            "-pix_fmt",
+            "rgb24",
+            "-s",
+            &format!("{}x{}", WIDTH, HEIGHT),
+            "-r",
+            &FPS.to_string(),
+            "-i",
+            "-",
+            "-c:v",
+            "libx264",
+            "-preset",
+            "slow",
+            "-crf",
+            "18",
+            "-pix_fmt",
+            "yuv420p",
+            OUTPUT,
+        ])
+        .stdin(Stdio::piped())
+        .spawn()?;
+
+    let stdin = ffmpeg.stdin.as_mut().unwrap();
+
+    for frame_i in 0..TOTAL_FRAMES {
+        for _ in 0..SUBSTEPS_PER_FRAME {
+            cells.par_iter_mut().for_each(|c| {
+                c.a = rk4(c.a);
+                c.b = rk4(c.b);
+
+                let dth1 = wrap(c.a.th1 - c.b.th1);
+                let dth2 = wrap(c.a.th2 - c.b.th2);
+
+                let dist = (dth1 * dth1 + dth2 * dth2).sqrt();
+
+                if dist > 0.0 {
+                    c.lyap += (dist / EPS).ln();
+
+                    let scale = EPS / dist;
+
+                    c.b.th1 = c.a.th1 + dth1 * scale;
+                    c.b.th2 = c.a.th2 + dth2 * scale;
+                    c.b.w1 = c.a.w1;
+                    c.b.w2 = c.a.w2;
+                }
+            });
+        }
+
+        render(&cells, &mut frame);
+
+        stdin.write_all(&frame)?;
+
+        println!("frame {}/{}", frame_i + 1, TOTAL_FRAMES);
+    }
+
+    drop(ffmpeg.stdin.take());
+
+    ffmpeg.wait()?;
+
+    Ok(())
+}
+
+fn initialize() -> Vec<Cell> {
+    let mut v = vec![
+        Cell {
+            a: State {
+                th1: 0.0,
+                th2: 0.0,
+                w1: 0.0,
+                w2: 0.0
+            },
+            b: State {
+                th1: 0.0,
+                th2: 0.0,
+                w1: 0.0,
+                w2: 0.0
+            },
+            lyap: 0.0
+        };
+        WIDTH * HEIGHT
+    ];
+
+    v.par_chunks_mut(WIDTH).enumerate().for_each(|(y, row)| {
+        let fy = y as f32 / (HEIGHT - 1) as f32;
+
+        for (x, c) in row.iter_mut().enumerate() {
+            let fx = x as f32 / (WIDTH - 1) as f32;
+
+            let phi1 = -PI + fx * 2.0 * PI;
+            let phi2 = -PI + fy * 2.0 * PI;
+
+            let a = State {
+                th1: phi1,
+                th2: phi2,
+                w1: 0.0,
+                w2: 0.0,
+            };
+
+            let b = State {
+                th1: phi1 + EPS,
+                th2: phi2 + EPS,
+                w1: 0.0,
+                w2: 0.0,
+            };
+
+            *c = Cell { a, b, lyap: 0.0 };
+        }
+    });
+
+    v
+}
+
+fn render(cells: &[Cell], frame: &mut [u8]) {
+    frame
+        .par_chunks_mut(WIDTH * 3)
+        .enumerate()
+        .for_each(|(y, row)| {
+            let src = &cells[y * WIDTH..(y + 1) * WIDTH];
+
+            for (x, c) in src.iter().enumerate() {
+                let i = x * 3;
+
+                let val = (c.lyap * 0.02).tanh();
+
+                let v = (val * 255.0) as u8;
+
+                row[i] = v;
+                row[i + 1] = v;
+                row[i + 2] = v;
+            }
+        });
+}
+
+fn rk4(s: State) -> State {
+    let k1 = deriv(s);
+
+    let s2 = add(s, k1, 0.5 * DT);
+    let k2 = deriv(s2);
+
+    let s3 = add(s, k2, 0.5 * DT);
+    let k3 = deriv(s3);
+
+    let s4 = add(s, k3, DT);
+    let k4 = deriv(s4);
+
+    State {
+        th1: s.th1 + DT * (k1.th1 + 2.0 * k2.th1 + 2.0 * k3.th1 + k4.th1) / 6.0,
+        th2: s.th2 + DT * (k1.th2 + 2.0 * k2.th2 + 2.0 * k3.th2 + k4.th2) / 6.0,
+        w1: s.w1 + DT * (k1.w1 + 2.0 * k2.w1 + 2.0 * k3.w1 + k4.w1) / 6.0,
+        w2: s.w2 + DT * (k1.w2 + 2.0 * k2.w2 + 2.0 * k3.w2 + k4.w2) / 6.0,
+    }
+}
+
+#[derive(Copy, Clone)]
+struct Deriv {
+    th1: f32,
+    th2: f32,
+    w1: f32,
+    w2: f32,
+}
+
+fn deriv(s: State) -> Deriv {
+    let d = s.th1 - s.th2;
+
+    let sin = d.sin();
+    let cos = d.cos();
+
+    let denom = 3.0 - (2.0 * d).cos();
+
+    let w1dot = (-9.81 * (3.0 * s.th1.sin() + (s.th1 - 2.0 * s.th2).sin())
+        - 2.0 * sin * (s.w2 * s.w2 + s.w1 * s.w1 * cos))
+        / denom;
+
+    let w2dot =
+        (2.0 * sin * (2.0 * s.w1 * s.w1 + 2.0 * 9.81 * s.th1.cos() + s.w2 * s.w2 * cos)) / denom;
+
+    Deriv {
+        th1: s.w1,
+        th2: s.w2,
+        w1: w1dot,
+        w2: w2dot,
+    }
+}
+
+fn add(s: State, k: Deriv, h: f32) -> State {
+    State {
+        th1: s.th1 + k.th1 * h,
+        th2: s.th2 + k.th2 * h,
+        w1: s.w1 + k.w1 * h,
+        w2: s.w2 + k.w2 * h,
+    }
+}
+
+fn wrap(a: f32) -> f32 {
+    let tau = 2.0 * PI;
+
+    (a + PI).rem_euclid(tau) - PI
+}