diff --git a/Cargo.toml b/Cargo.toml
index 2cddf90..f24ac2c 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -12,3 +12,4 @@ decorum = "0.3"
 enum_dispatch = "0.3"
 ndarray = "0.13"
 piecewise-linear = "0.1"
+y4m = "0.7"
diff --git a/examples/em_reflection.rs b/examples/em_reflection.rs
index 11be6a4..6a11917 100644
--- a/examples/em_reflection.rs
+++ b/examples/em_reflection.rs
@@ -1,5 +1,5 @@
 use coremem::{consts, mat, SimState};
-use coremem::render::ColorTermRenderer as Renderer;
+use coremem::render;
 use std::{thread, time};
 
 fn main() {
@@ -26,6 +26,7 @@ fn main() {
     }
 
     let mut step = 0u64;
+    let mut renderer = render::Y4MRenderer::new("test.y4m");
     loop {
         step += 1;
         let imp = if step < 50 { 
@@ -40,7 +41,7 @@ fn main() {
         for x in 0..width {
             state.impulse_bz(x, 20, (imp / 3.0e8) as _);
         }
-        Renderer.render(&state);
+        renderer.render(&state);
         state.step();
         thread::sleep(time::Duration::from_millis(67));
     }
diff --git a/src/render.rs b/src/render.rs
index 7c5e277..7e93695 100644
--- a/src/render.rs
+++ b/src/render.rs
@@ -1,6 +1,9 @@
 use ansi_term::Color::RGB;
 use crate::{consts, Material as _, SimState};
 use std::fmt::Write as _;
+use std::fs::File;
+use std::path::PathBuf;
+use y4m::{Colorspace, encode, Encoder, Frame, Ratio};
 
 pub struct NumericTermRenderer;
 
@@ -63,6 +66,57 @@ impl ColorTermRenderer {
             }
             write!(&mut buf, "\n");
         }
-        println!("{}\ntime: {:.3e}", buf, state.time());
+        println!("{}\ntime: {:.3e} (fr {})", buf, state.time(), state.step_no());
+    }
+}
+
+pub struct Y4MRenderer {
+    out_path: PathBuf,
+    encoder: Option<Encoder<File>>,
+}
+
+impl Y4MRenderer {
+    pub fn new<S: Into<PathBuf>>(output: S) -> Self {
+        Self {
+            out_path: output.into(),
+            encoder: None,
+        }
+    }
+    pub fn render(&mut self, state: &SimState) {
+        if self.encoder.is_none() {
+            let writer = File::create(&self.out_path).unwrap();
+            self.encoder = Some(encode(state.width(), state.height(), Ratio::new(30, 1))
+                .with_colorspace(Colorspace::C444)
+                .write_header(writer)
+                .unwrap()
+            );
+        }
+
+        let mut pix_y = Vec::new();
+        let mut pix_u = Vec::new();
+        let mut pix_v = Vec::new();
+        for y in 0..state.height() {
+            for x in 0..state.width() {
+                let cell = state.get(x, y);
+                //let r = norm_color(cell.bz() * consts::C);
+                //let r = 0;
+                let r = norm_color(cell.mat().mz()*1.0e-2);
+                let b = (55.0*cell.mat().conductivity()).min(255.0) as u8;
+                //let b = 0;
+                //let b = norm_color(cell.ey());
+                //let g = 0;
+                //let g = norm_color(cell.ex());
+                //let g = norm_color(curl(cell.ex(), cell.ey()));
+                let g = norm_color((cell.bz() * 1.0e4).into());
+                //let g = norm_color(cell.ey().into());
+                pix_y.push(r);
+                pix_u.push(g);
+                pix_v.push(b);
+            }
+        }
+
+        let frame = Frame::new([&*pix_y, &*pix_u, &*pix_v], None);
+        let enc = self.encoder.as_mut().unwrap();
+        enc.write_frame(&frame).unwrap()
     }
 }
diff --git a/src/sim.rs b/src/sim.rs
index 4f02af6..b6e4d86 100644
--- a/src/sim.rs
+++ b/src/sim.rs
@@ -24,6 +24,10 @@ impl SimState {
         (self.timestep() * self.step_no as f64).into()
     }
 
+    pub fn step_no(&self) -> u64 {
+        self.step_no
+    }
+
     pub fn step(&mut self) {
         use consts::real::*;
         let half_time_step = HALF() * self.timestep();
@@ -193,10 +197,12 @@ impl<M: Material> Cell<M> {
 
         // XXX not obvious that bz_to_hz is sensible
         let delta_hz_y = self.hz() - self.bz_to_hz(up.bz());
+        // let delta_hz_y = self.hz() - up.hz();
         let ex_rhs = self.state.ex*(ONE() - sigma/EPS0()*delta_t) + TWO()*delta_t/EPS0()*delta_hz_y/feature_size;
         let ex_next = ex_rhs / (ONE() + sigma/EPS0()*delta_t);
 
         let delta_hz_x = self.hz() - self.bz_to_hz(left.bz());
+        // let delta_hz_x = self.hz() - left.hz();
         let ey_rhs = self.state.ey*(ONE() - sigma/EPS0()*delta_t) - TWO()*delta_t/EPS0()*delta_hz_x/feature_size;
         let ey_next = ey_rhs / (ONE() + sigma/EPS0()*delta_t);