///|
pub(all) struct ConservationReport {
  initial_charge : Double
  final_charge : Double
  charge_drift : Double
  initial_energy : Double
  final_energy : Double
  energy_drift : Double
  particle_count_preserved : Bool
} derive(Debug, ToJson)

///|
pub fn relative_change(initial : Double, final_value : Double) -> Double {
  safe_ratio((final_value - initial).abs(), initial.abs().max(1.0e-30), 0.0)
}

///|
pub fn conservation_report(
  charges : ArrayView[Double],
  energies : ArrayView[Double],
) -> ConservationReport {
  let initial_charge = if charges.length() == 0 { 0.0 } else { charges[0] }
  let final_charge = if charges.length() == 0 {
    0.0
  } else {
    charges[charges.length() - 1]
  }
  let initial_energy = if energies.length() == 0 { 0.0 } else { energies[0] }
  let final_energy = if energies.length() == 0 {
    0.0
  } else {
    energies[energies.length() - 1]
  }
  let mut maximum_energy_drift = 0.0
  for value in energies {
    let drift = relative_change(initial_energy, value)
    if drift > maximum_energy_drift {
      maximum_energy_drift = drift
    }
  }
  {
    initial_charge,
    final_charge,
    charge_drift: relative_change(initial_charge, final_charge),
    initial_energy,
    final_energy,
    energy_drift: maximum_energy_drift,
    particle_count_preserved: true,
  }
}

///|
pub fn charge_conservation_error(
  grid : Grid1D,
  densities : ArrayView[Array[Double]],
  expected : Double,
) -> Double {
  let mut maximum = 0.0
  for density in densities {
    let error = (total_charge(grid, density) - expected).abs()
    if error > maximum {
      maximum = error
    }
  }
  maximum
}

///|
pub fn energy_series(
  particles : ArrayView[Array[Particle]],
  fields : ArrayView[Field1D],
) -> Array[Double] {
  let count = if particles.length() < fields.length() {
    particles.length()
  } else {
    fields.length()
  }
  Array::makei(count, fn(i) {
    kinetic_energy(particles[i]) + field_energy(fields[i])
  })
}

///|
pub fn charge_series(
  grid : Grid1D,
  densities : ArrayView[Array[Double]],
) -> Array[Double] {
  densities.map(fn(density) { total_charge(grid, density) })
}

///|
pub fn momentum_series(particles : ArrayView[Array[Particle]]) -> Array[Double] {
  particles.map(fn(values) {
    values.fold(init=0.0, fn(acc, particle) {
      acc + particle.mass * particle.v * particle.weight
    })
  })
}

///|
pub fn mass_series(particles : ArrayView[Array[Particle]]) -> Array[Double] {
  particles.map(fn(values) {
    values.fold(init=0.0, fn(acc, particle) {
      acc + particle.mass * particle.weight
    })
  })
}

///|
pub fn series_range(values : ArrayView[Double]) -> Double {
  max_value(values, 0.0) - min_value(values, 0.0)
}

///|
pub fn series_mean(values : ArrayView[Double]) -> Double {
  mean(values)
}

///|
pub fn series_rms(values : ArrayView[Double]) -> Double {
  if values.length() == 0 {
    0.0
  } else {
    (values.fold(init=0.0, fn(acc, value) { acc + value * value }) /
    values.length().to_double()).sqrt()
  }
}

///|
pub fn monotonic_time(values : ArrayView[Double]) -> Bool {
  let mut valid = true
  for i in 1.. Array[Double] {
  if values.length() == 0 {
    []
  } else {
    values.map(fn(value) { relative_change(values[0], value) })
  }
}