///|
pub(all) struct Moments1D {
  density : Array[Double]
  momentum : Array[Double]
  temperature_like : Array[Double]
} derive(Debug, ToJson)

///|
pub fn particle_number_density(
  grid : Grid1D,
  particles : ArrayView[Particle],
) -> Array[Double] {
  let density = zeros(grid.cells)
  for p in particles {
    let y = grid.wrap(p.x)
    let scaled = y / grid.dx - 0.5
    let i0 = @math.floor(scaled).to_int()
    let frac = scaled - i0.to_double()
    let left = if i0 < 0 { grid.cells - 1 } else { i0 % grid.cells }
    let right = (left + 1) % grid.cells
    let n = p.weight / grid.dx
    density[left] = density[left] + n * (1.0 - frac)
    density[right] = density[right] + n * frac
  }
  density
}

///|
pub fn particle_current_density(
  grid : Grid1D,
  particles : ArrayView[Particle],
) -> Array[Double] {
  let current = zeros(grid.cells)
  for p in particles {
    let y = grid.wrap(p.x)
    let scaled = y / grid.dx - 0.5
    let i0 = @math.floor(scaled).to_int()
    let frac = scaled - i0.to_double()
    let left = if i0 < 0 { grid.cells - 1 } else { i0 % grid.cells }
    let right = (left + 1) % grid.cells
    let j = p.charge * p.weight * p.v / grid.dx
    current[left] = current[left] + j * (1.0 - frac)
    current[right] = current[right] + j * frac
  }
  current
}

///|
pub fn total_charge(
  grid : Grid1D,
  charge_density : ArrayView[Double],
) -> Double {
  charge_density.fold(init=0.0, fn(acc, rho) { acc + rho * grid.dx })
}

///|
pub fn total_particle_charge(particles : ArrayView[Particle]) -> Double {
  particles.fold(init=0.0, fn(acc, p) { acc + p.charge * p.weight })
}

///|
pub fn center_of_mass(grid : Grid1D, particles : ArrayView[Particle]) -> Double {
  let weighted = particles.fold(init=(0.0, 0.0), fn(acc, p) {
    let (sum_x, sum_w) = acc
    (sum_x + grid.wrap(p.x) * p.weight, sum_w + p.weight)
  })
  let (sum_x, sum_w) = weighted
  if sum_w == 0.0 {
    0.0
  } else {
    sum_x / sum_w
  }
}

///|
pub fn velocity_mean(particles : ArrayView[Particle]) -> Double {
  let weighted = particles.fold(init=(0.0, 0.0), fn(acc, p) {
    let (sum_v, sum_w) = acc
    (sum_v + p.v * p.weight, sum_w + p.weight)
  })
  let (sum_v, sum_w) = weighted
  if sum_w == 0.0 {
    0.0
  } else {
    sum_v / sum_w
  }
}

///|
pub fn velocity_variance(particles : ArrayView[Particle]) -> Double {
  let u = velocity_mean(particles)
  let weighted = particles.fold(init=(0.0, 0.0), fn(acc, p) {
    let (sum, sum_w) = acc
    let dv = p.v - u
    (sum + dv * dv * p.weight, sum_w + p.weight)
  })
  let (sum, sum_w) = weighted
  if sum_w == 0.0 {
    0.0
  } else {
    sum / sum_w
  }
}

///|
pub fn vlasov_moments(state : VlasovState) -> Moments1D {
  let c = state.config
  let density = zeros(c.grid.cells)
  let momentum = zeros(c.grid.cells)
  let temperature_like = zeros(c.grid.cells)
  let dv = c.dv()
  for ix in 0..