///|
pub(all) struct PoissonResult {
  potential : Array[Double]
  residual : Double
  iterations : Int
  converged : Bool
} derive(Debug, ToJson)

///|
pub fn poisson_source(charge_density : ArrayView[Double]) -> Array[Double] {
  let average = mean(charge_density)
  charge_density.map(fn(value) { -(value - average) / vacuum_permittivity })
}

///|
pub fn poisson_residual(
  grid : Grid1D,
  potential : ArrayView[Double],
  source : ArrayView[Double],
) -> Double {
  if potential.length() != grid.cells ||
    source.length() != grid.cells ||
    grid.cells == 0 {
    0.0
  } else {
    let mut total = 0.0
    let mut scale = 0.0
    for i in 0.. PoissonResult {
  let source = poisson_source(charge_density)
  let potential = zeros(grid.cells)
  let mut residual = poisson_residual(grid, potential, source)
  let mut iterations = 0
  let limit = if max_iterations < 0 { 0 } else { max_iterations }
  while iterations < limit && residual > tolerance {
    for i in 0.. PoissonResult {
  let source = poisson_source(charge_density)
  let current = zeros(grid.cells)
  let next = zeros(grid.cells)
  let mut residual = poisson_residual(grid, current, source)
  let mut iterations = 0
  while iterations < max_iterations && residual > tolerance {
    for i in 0.. Array[Double] {
  let potential = zeros(grid.cells)
  if grid.cells >= 1 {
    potential[0] = left_value
    potential[grid.cells - 1] = right_value
  }
  if grid.cells > 2 {
    for i in 1..<(grid.cells - 1) {
      let weight = i.to_double() / (grid.cells - 1).to_double()
      potential[i] = lerp(left_value, right_value, weight) -
        source[i] * grid.dx * grid.dx * 0.5
    }
  }
  potential
}

///|
pub fn poisson_electric_field(
  grid : Grid1D,
  potential : ArrayView[Double],
) -> Array[Double] {
  let electric = zeros(grid.cells)
  for i in 0.. Double {
  let electric = poisson_electric_field(grid, potential)
  0.5 *
  vacuum_permittivity *
  grid.dx *
  sum_values(electric.map(fn(value) { value * value }))
}

///|
pub fn poisson_mean_free(potential : ArrayView[Double]) -> Array[Double] {
  let average = mean(potential)
  potential.map(fn(value) { value - average })
}

///|
pub fn poisson_iteration_history(
  grid : Grid1D,
  source : ArrayView[Double],
  steps : Int,
) -> Array[Double] {
  let potential = zeros(grid.cells)
  let history : Array[Double] = []
  let mut step = 0
  while step < steps {
    for i in 0.. Bool {
  solve_poisson_periodic(grid, source, max_iterations, tolerance).converged
}

///|
pub fn poisson_residual_series(
  grid : Grid1D,
  potentials : ArrayView[Array[Double]],
  source : ArrayView[Double],
) -> Array[Double] {
  potentials.map(fn(potential) { poisson_residual(grid, potential, source) })
}

///|
pub fn poisson_solution_range(result : PoissonResult) -> Double {
  max_value(result.potential, 0.0) - min_value(result.potential, 0.0)
}

///|
/// Compatibility wrapper with a source-density API.
pub fn solve_periodic_poisson(
  grid : Grid1D,
  source : ArrayView[Double],
  max_iterations? : Int = 1000,
  tolerance? : Double = 1.0e-8,
) -> PoissonResult {
  let potential = zeros(source.length())
  let target = source.map(fn(value) { -value })
  let mut residual = poisson_residual(grid, potential, target)
  let mut iterations = 0
  while iterations < max_iterations && residual > tolerance {
    for i in 0.. Field1D {
  let result = solve_periodic_poisson(grid, source, max_iterations~, tolerance~)
  let electric = periodic_gradient(grid, result.potential).map(fn(value) {
    -value
  })
  {
    grid,
    charge_density: source.to_owned(),
    electric,
    potential: result.potential,
  }
}