diff --git a/src/bin/audio.rs b/src/bin/audio.rs
new file mode 100644
index 0000000..5eaed9e
--- /dev/null
+++ b/src/bin/audio.rs
@@ -0,0 +1,117 @@
+use symplectic_double_pendulum::*;
+
+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(())
+}