///|
/// Build a D2Q9 equilibrium vector.
pub fn equilibrium_distribution(
  rho~ : Double,
  ux~ : Double,
  uy~ : Double,
) -> Array[Double] {
  let result = Array::make(q, 0.0)
  for d in 0.. Double {
  let mut total = 0.0
  for d in 0.. Double {
  let mut total = 0.0
  for d in 0.. Double {
  let mut total = 0.0
  for d in 0.. Cell {
  let rho = distribution_mass(distribution)
  if rho <= 0.0 {
    { rho, ux: 0.0, uy: 0.0 }
  } else {
    {
      rho,
      ux: distribution_momentum_x(distribution) / rho,
      uy: distribution_momentum_y(distribution) / rho,
    }
  }
}

///|
/// Collision coefficients used by the conservative MRT-style operator.
fn collision_factor(model : CollisionModel, direction : Int) -> Double {
  match model {
    Bgk(omega) => 1.0 - omega
    Regularized(omega) =>
      if direction >= 5 {
        1.0 - omega * 0.75
      } else {
        1.0 - omega
      }
    Mrt(omega) =>
      if direction == 0 {
        1.0
      } else if direction < 5 {
        1.0 - omega
      } else {
        1.0 - (omega * 1.25).min(1.95)
      }
  }
}

///|
/// Correct a post-collision vector so its conserved moments remain exact.
fn conserve_distribution(
  candidate : Array[Double],
  target_mass : Double,
  target_momentum_x : Double,
  target_momentum_y : Double,
) -> Array[Double] {
  let mass_delta = target_mass - distribution_mass(candidate)
  let x_delta = target_momentum_x - distribution_momentum_x(candidate)
  let y_delta = target_momentum_y - distribution_momentum_y(candidate)
  for d in 0.. Array[Double] {
  let equilibrium_values = equilibrium_distribution(rho~, ux~, uy~)
  let candidate = Array::make(q, 0.0)
  for d in 0.. Array[Double] {
  let equilibrium_values = equilibrium_distribution(rho~, ux~, uy~)
  let stress = Array::make(3, 0.0)
  for d in 0.. Double {
  (1.0 / model.omega() - 0.5) / 3.0
}

///|
/// Estimate the lattice sound speed.
pub fn lattice_sound_speed() -> Double {
  (1.0 / 3.0).sqrt()
}

///|
/// Compute the lattice Mach number of a velocity vector.
pub fn mach_number(ux~ : Double, uy~ : Double) -> Double {
  (ux * ux + uy * uy).sqrt() / lattice_sound_speed()
}