///|
pub(all) struct VlasovConfig {
  grid : Grid1D
  velocity_min : Double
  velocity_max : Double
  velocity_bins : Int
  dt : Double
} derive(Debug, ToJson)

///|
pub(all) struct VlasovState {
  config : VlasovConfig
  distribution : Array[Double]
  charge_density : Array[Double]
  time : Double
} derive(Debug, ToJson)

///|
pub fn VlasovConfig::new(
  grid : Grid1D,
  velocity_min~ : Double,
  velocity_max~ : Double,
  velocity_bins~ : Int,
  dt~ : Double,
) -> VlasovConfig {
  { grid, velocity_min, velocity_max, velocity_bins, dt }
}

///|
pub fn VlasovConfig::dv(config : VlasovConfig) -> Double {
  (config.velocity_max - config.velocity_min) / config.velocity_bins.to_double()
}

///|
pub fn VlasovConfig::velocity(config : VlasovConfig, bin : Int) -> Double {
  config.velocity_min + (bin.to_double() + 0.5) * config.dv()
}

///|
pub fn maxwellian(v : Double, thermal_speed : Double) -> Double {
  let norm = 1.0 / ((2.0 * pi).sqrt() * thermal_speed)
  norm * @math.exp(-0.5 * v * v / (thermal_speed * thermal_speed))
}

///|
pub fn make_landau_initial(
  config : VlasovConfig,
  density0~ : Double,
  thermal_speed~ : Double,
  perturbation? : Double = 0.01,
  mode? : Int = 1,
) -> VlasovState {
  let nx = config.grid.cells
  let nv = config.velocity_bins
  let dist = Array::make(nx * nv, 0.0)
  for ix in 0.. Array[Double] {
  let rho = zeros(config.grid.cells)
  let dv = config.dv()
  for ix in 0.. VlasovState {
  let c = state.config
  let nx = c.grid.cells
  let nv = c.velocity_bins
  let next = Array::make(nx * nv, 0.0)
  for ix in 0..