Add license and README

.gitignore

@@ -1 +1,3 @@
 /target
+*.csv
+*.wav

LICENSE.md

@@ -0,0 +1,108 @@
+# Komorebi License
+
+Version 2.0.0
+
+## Acceptance
+
+In order to get any license under these terms, you must agree
+to them as both strict obligations and conditions to all
+your licenses.
+
+## Copyright License
+
+The licensor grants you a copyright license for the software
+to do everything you might do with the software that would
+otherwise infringe the licensor's copyright in it for any
+permitted purpose. However, you may only distribute the source
+code of the software according to the [Distribution License](
+#distribution-license), you may only make changes according
+to the [Changes License](#changes-license), and you may not
+otherwise distribute the software or new works based on the
+software.
+
+## Distribution License
+
+The licensor grants you an additional copyright license to
+distribute copies of the source code of the software. Your
+license to distribute covers distributing the source code of
+the software with changes permitted by the [Changes License](
+#changes-license).
+
+## Changes License
+
+The licensor grants you an additional copyright license to
+make changes for any permitted purpose.
+
+## Patent License
+
+The licensor grants you a patent license for the software that
+covers patent claims the licensor can license, or becomes able
+to license, that you would infringe by using the software.
+
+## Personal Uses
+
+Personal use for research, experiment, and testing for
+the benefit of public knowledge, personal study, private
+entertainment, hobby projects, amateur pursuits, or religious
+observance, without any anticipated commercial application,
+is use for a permitted purpose.
+
+## Fair Use
+
+You may have "fair use" rights for the software under the
+law. These terms do not limit them.
+
+## No Other Rights
+
+These terms do not allow you to sublicense or transfer any of
+your licenses to anyone else, or prevent the licensor from
+granting licenses to anyone else.  These terms do not imply
+any other licenses.
+
+## Patent Defense
+
+If you make any written claim that the software infringes or
+contributes to infringement of any patent, your patent license
+for the software granted under these terms ends immediately. If
+your company makes such a claim, your patent license ends
+immediately for work on behalf of your company.
+
+## Violations
+
+The first time you are notified in writing that you have
+violated any of these terms, or done anything with the software
+not covered by your licenses, your licenses can nonetheless
+continue if you come into full compliance with these terms,
+and take practical steps to correct past violations, within
+32 days of receiving notice.  Otherwise, all your licenses
+end immediately.
+
+## No Liability
+
+***As far as the law allows, the software comes as is, without
+any warranty or condition, and the licensor will not be liable
+to you for any damages arising out of these terms or the use
+or nature of the software, under any kind of legal claim.***
+
+## Definitions
+
+The **licensor** is the individual or entity offering these
+terms, and the **software** is the software the licensor makes
+available under these terms.
+
+**You** refers to the individual or entity agreeing to these
+terms.
+
+**Your company** is any legal entity, sole proprietorship,
+or other kind of organization that you work for, plus all
+organizations that have control over, are under the control of,
+or are under common control with that organization.  **Control**
+means ownership of substantially all the assets of an entity,
+or the power to direct its management and policies by vote,
+contract, or otherwise.  Control can be direct or indirect.
+
+**Your licenses** are all the licenses granted to you for the
+software under these terms.
+
+**Use** means anything you do with the software requiring one
+of your licenses.

README.md

@@ -0,0 +1,23 @@
+# Symplectic Double Pendulum Simulation
+
+A Hamitonian-based double pendulum simulation, using the augmented
+phase space method decribed in [Molei Tao's paper on the
+subject](https://arxiv.org/pdf/1609.02212).
+
+## Building and Running
+
+The simulation is written in Rust and uses Cargo as its build system.
+The Rust toolchain (including Cargo) is most easily installed with
+[Rustup](https://rustup.rs/).
+
+Given that you have a toolchain installed, a real-time animation of
+the simulation can be ran with:
+
+```
+cargo run --bin anim
+```
+
+## License
+
+Copyright © Camden Dixie O'Brien, available under the Komorebi 2.0.0
+license.  See [LICENSE.md](LICENSE.md) for details.

plot.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+import sys
+
+import matplotlib.pyplot as plt
+import pandas as pd
+
+
+def main():
+    src = sys.argv[1] if len(sys.argv) > 1 else sys.stdin
+    df = pd.read_csv(src)
+
+    fig, ax = plt.subplots(figsize=(10, 5))
+    ax.plot(df["time"], df["tao"], label="Tao (symplectic)", linewidth=1)
+    ax.plot(df["time"], df["rk4"], label="RK4", linewidth=1)
+
+    ax.set_xlabel("Time")
+    ax.set_ylabel("Energy")
+    ax.set_title("Double pendulum energy: RK4 vs Tao's method")
+    ax.legend()
+    ax.grid(True, alpha=0.3)
+
+    fig.tight_layout()
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

src/bin/anim.rs

@@ -0,0 +1,44 @@
+use dpsim::*;
+
+use macroquad::prelude::*;
+
+fn draw(dp: &DoublePendulum, fg: Color) {
+	let rod_length = 120.0;
+
+	let pivot = vec2(screen_width() / 2.0, screen_height() / 2.0);
+	let bob1 = vec2(
+		pivot.x + rod_length * (dp.theta.sin() as f32),
+		pivot.y + rod_length * (dp.theta.cos() as f32));
+	let bob2 = vec2(
+		bob1.x + rod_length * ((dp.theta + dp.phi).sin() as f32),
+		bob1.y + rod_length * ((dp.theta + dp.phi).cos() as f32));
+
+	draw_line(pivot.x, pivot.y, bob1.x, bob1.y, 2.0, fg);
+	draw_line(bob1.x, bob1.y, bob2.x, bob2.y, 2.0, fg);
+
+	draw_circle(bob1.x, bob1.y, 8.0, fg);
+	draw_circle(bob2.x, bob2.y, 8.0, fg);
+}
+
+#[macroquad::main("Double Pendulum")]
+async fn main() {
+	let omega = 500.0;
+	let h = 1e-3;
+
+	let mut state = Augmented::from(State {
+		q: DoublePendulum { theta: std::f64::consts::PI - 0.01, phi: 0.0 },
+		p: DoublePendulum { theta: 0.0, phi: 0.0 },
+	}, omega);
+
+	loop {
+		clear_background(WHITE);
+
+		let t = get_frame_time() as f64;
+		for _ in 0..((t / h) as i32) {
+			state.evolve(h);
+		}
+		draw(&state.average().q, BLACK);
+
+		next_frame().await
+	}
+}

src/bin/audio.rs

@@ -0,0 +1,117 @@
+use dpsim::*;
+
+use std::fs::File;
+use std::io::{BufWriter, Write};
+
+const RATE: u32 = 44100;
+const CHANNELS: u16 = 2;
+const DEPTH: u16 = 16;
+
+const DURATION: u32 = 60;
+
+fn write_wav_header<O: Write>(out: &mut O, n: u32) -> std::io::Result<()> {
+	let block_align = CHANNELS * DEPTH / 8;
+	let byte_rate = RATE * block_align as u32;
+	let data_size = n * block_align as u32;
+	let file_size = data_size + 36;
+
+	out.write_all(b"RIFF")?;
+	out.write_all(&file_size.to_le_bytes())?;
+	out.write_all(b"WAVE")?;
+
+	out.write_all(b"fmt ")?;
+	out.write_all(&16u32.to_le_bytes())?;
+	out.write_all(&1u16.to_le_bytes())?;
+	out.write_all(&CHANNELS.to_le_bytes())?;
+	out.write_all(&RATE.to_le_bytes())?;
+	out.write_all(&byte_rate.to_le_bytes())?;
+	out.write_all(&block_align.to_le_bytes())?;
+	out.write_all(&DEPTH.to_le_bytes())?;
+
+	out.write_all(b"data")?;
+	out.write_all(&data_size.to_le_bytes())?;
+
+	Ok(())
+}
+
+fn write_samples<O: Write>(
+	out: &mut O, samples: (u16, u16)) -> std::io::Result<()> {
+	out.write_all(&samples.0.to_le_bytes())?;
+	out.write_all(&samples.1.to_le_bytes())?;
+	Ok(())
+}
+
+fn bob1_x(state: &State<DoublePendulum>) -> f64 {
+	state.q.theta.sin()
+}
+
+fn bob1_y(state: &State<DoublePendulum>) -> f64 {
+	state.q.theta.cos()
+}
+
+fn diff_x(state: &State<DoublePendulum>) -> f64 {
+	(state.q.theta + state.q.phi).sin()
+}
+
+fn diff_y(state: &State<DoublePendulum>) -> f64 {
+	(state.q.theta + state.q.phi).cos()
+}
+
+fn bob2_x(state: &State<DoublePendulum>) -> f64 {
+	bob1_x(state) + diff_x(state)
+}
+
+fn bob2_y(state: &State<DoublePendulum>) -> f64 {
+	bob1_y(state) + diff_y(state)
+}
+
+fn sample(x: f64, min: f64, max: f64) -> u16 {
+	let norm = (x - min) / (max - min);
+	(u16::MAX as f64 * norm) as u16
+}
+
+fn bob1_samples(state: &State<DoublePendulum>) -> (u16, u16) {
+	let (min, max) = (-1.0, 1.0);
+	let (x, y) = (bob1_x(state), bob1_y(state));
+	(sample(x, min, max), sample(y, min, max))
+}
+
+fn diff_samples(state: &State<DoublePendulum>) -> (u16, u16) {
+	let (min, max) = (-1.0, 1.0);
+	let (x, y) = (diff_x(state), diff_y(state));
+	(sample(x, min, max), sample(y, min, max))
+}
+
+fn bob2_samples(state: &State<DoublePendulum>) -> (u16, u16) {
+	let (min, max) = (-2.0, 2.0);
+	let (x, y) = (bob2_x(state), bob2_y(state));
+	(sample(x, min, max), sample(y, min, max))
+}
+
+fn main() -> std::io::Result<()> {
+	let mut bob1 = BufWriter::new(File::create("bob1.wav")?);
+	let mut bob2 = BufWriter::new(File::create("bob2.wav")?);
+	let mut diff = BufWriter::new(File::create("diff.wav")?);
+
+	let n = RATE * DURATION;
+	write_wav_header(&mut bob1, n)?;
+	write_wav_header(&mut bob2, n)?;
+	write_wav_header(&mut diff, n)?;
+
+	let omega = 500.0;
+	let h = 1e-3;
+
+	let mut state = Augmented::from(State {
+		q: DoublePendulum { theta: std::f64::consts::PI - 0.01, phi: 0.0 },
+		p: DoublePendulum { theta: 0.0, phi: 0.0 },
+	}, omega);
+	for _ in 0..n {
+		state.evolve(h);
+
+		write_samples(&mut bob1, bob1_samples(&state.average()))?;
+		write_samples(&mut bob2, bob2_samples(&state.average()))?;
+		write_samples(&mut diff, diff_samples(&state.average()))?;
+	}
+
+	Ok(())
+}

src/bin/comp_rk4.rs

@@ -0,0 +1,41 @@
+use dpsim::*;
+
+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());
+		}
+	}
+}