///|
/// A reusable D2Q9 simulation with selectable collision and outer boundaries.
pub(all) struct Simulation {
  size : Size
  options : SimulationOptions
  f : Array[Double]
  next : Array[Double]
  solid : Array[Bool]
  mut step_count : Int
} derive(Debug)

///|
/// Allocate an equilibrium advanced simulation.
pub fn Simulation::new(
  size~ : Size,
  options~ : SimulationOptions,
) -> Simulation {
  let cells = size.width.max(0) * size.height.max(0)
  let f = Array::make(cells * q, 0.0)
  for id in 0.. Size {
  self.size
}

///|
/// Return whether integer coordinates are inside the simulation.
pub fn Simulation::contains(self : Simulation, x : Int, y : Int) -> Bool {
  x >= 0 && y >= 0 && x < self.size.width && y < self.size.height
}

///|
/// Convert coordinates and direction to a distribution index.
pub fn Simulation::distribution_index(
  self : Simulation,
  x : Int,
  y : Int,
  d : Int,
) -> Int {
  (y * self.size.width + x) * q + d
}

///|
/// Set a solid flag inside the simulation.
pub fn Simulation::set_solid(
  self : Simulation,
  x~ : Int,
  y~ : Int,
  solid~ : Bool,
) -> Unit {
  if self.contains(x, y) {
    self.solid[y * self.size.width + x] = solid
  }
}

///|
/// Read a solid flag; outside coordinates are solid.
pub fn Simulation::is_solid(self : Simulation, x~ : Int, y~ : Int) -> Bool {
  if self.contains(x, y) {
    self.solid[y * self.size.width + x]
  } else {
    true
  }
}

///|
/// Reset a cell to an equilibrium state.
pub fn Simulation::set_cell(
  self : Simulation,
  x~ : Int,
  y~ : Int,
  rho~ : Double,
  ux~ : Double,
  uy~ : Double,
) -> Unit {
  if self.contains(x, y) {
    let values = equilibrium_distribution(rho~, ux~, uy~)
    let base = self.distribution_index(x, y, 0)
    for d in 0.. Unit {
  for y in 0.. Cell {
  if self.is_solid(x~, y~) {
    { rho: 1.0, ux: 0.0, uy: 0.0 }
  } else if !self.contains(x, y) {
    { rho: 0.0, ux: 0.0, uy: 0.0 }
  } else {
    let base = self.distribution_index(x, y, 0)
    distribution_cell(self.f[base:base + q])
  }
}

///|
/// Return the total fluid mass.
pub fn Simulation::mass(self : Simulation) -> Double {
  let mut total = 0.0
  for y in 0.. Array[Cell] {
  let result = Array::new()
  for y in 0.. Stability {
  let mut minimum = 1000000000.0
  let mut maximum = -1000000000.0
  let mut max_speed = 0.0
  let mut mass = 0.0
  let mut has_nan = false
  let mut fluid_cells = 0
  for y in 0.. self.options.density_floor &&
    mach_number(ux=max_speed, uy=0.0) < self.options.max_mach &&
    self.options.model.omega() > 0.0 &&
    self.options.model.omega() < 2.0,
  }
}

///|
/// Advance several steps.
pub fn Simulation::run(self : Simulation, steps : Int) -> Unit {
  for _ in 0.. Int {
  if size <= 0 {
    0
  } else {
    let remainder = value % size
    if remainder < 0 {
      remainder + size
    } else {
      remainder
    }
  }
}

///|
/// Resolve a streamed destination according to the outer boundary mode.
fn Simulation::destination(
  self : Simulation,
  x : Int,
  y : Int,
) -> (Int, Int, Bool) {
  if self.contains(x, y) {
    (x, y, true)
  } else {
    match self.options.boundary {
      Periodic =>
        (
          wrap_coordinate(x, self.size.width),
          wrap_coordinate(y, self.size.height),
          true,
        )
      BounceBack => (x, y, false)
      Open => (x, y, false)
    }
  }
}

///|
/// Stream one direction into a destination or reflect it locally.
fn Simulation::stream(
  self : Simulation,
  x : Int,
  y : Int,
  d : Int,
  value : Double,
) -> Unit {
  let (nx, ny, valid) = self.destination(x + cx[d], y + cy[d])
  if valid && !self.is_solid(x=nx, y=ny) {
    self.next[self.distribution_index(nx, ny, d)] = value
  } else {
    self.next[self.distribution_index(x, y, opposite[d])] = value
  }
}

///|
/// Advance one collision/stream/boundary step.
pub fn Simulation::step(self : Simulation) -> Unit {
  self.next.fill(0.0)
  for y in 0.. Array[Double] {
  let yy = y.clamp(min=0, max=self.size.height - 1)
  let result = Array::new()
  for x in 0.. Array[Double] {
  let xx = x.clamp(min=0, max=self.size.width - 1)
  let result = Array::new()
  for y in 0..