///|
/// Extract density values in row-major order.
pub fn Lattice::rho_field(self : Lattice) -> Array[Double] {
  let out = Array::new(capacity=self.size.width * self.size.height)
  for y in 0.. Array[Double] {
  let out = Array::new(capacity=self.size.width * self.size.height)
  for y in 0.. Array[Double] {
  let out = Array::new(capacity=self.size.width * self.size.height)
  for y in 0.. Array[Double] {
  let out = Array::new(capacity=self.size.width * self.size.height)
  for y in 0.. Array[Double] {
  let out = Array::new(capacity=self.size.width * self.size.height)
  for y in 0.. Double {
  if self.is_solid(x, y) {
    0.0
  } else {
    let xl = (x - 1).clamp(min=0, max=self.size.width - 1)
    let xr = (x + 1).clamp(min=0, max=self.size.width - 1)
    let yb = (y - 1).clamp(min=0, max=self.size.height - 1)
    let yt = (y + 1).clamp(min=0, max=self.size.height - 1)
    let dx = (xr - xl).max(1).to_double()
    let dy = (yt - yb).max(1).to_double()
    let duy_dx = (self.cell(xr, y).uy - self.cell(xl, y).uy) / dx
    let dux_dy = (self.cell(x, yt).ux - self.cell(x, yb).ux) / dy
    duy_dx - dux_dy
  }
}

///|
/// Sample `ux` along one horizontal row.
pub fn Lattice::horizontal_profile(self : Lattice, y~ : Int) -> Array[Double] {
  let yy = y.clamp(min=0, max=self.size.height - 1)
  let out = Array::new(capacity=self.size.width)
  for x in 0.. Array[Double] {
  let xx = x.clamp(min=0, max=self.size.width - 1)
  let out = Array::new(capacity=self.size.height)
  for y in 0..