///|
pub(all) struct PicSeriesSummary {
  samples : Int
  particles_initial : Int
  particles_final : Int
  charge_drift : Double
  energy_drift : Double
  maximum_field : Double
} derive(Debug, ToJson)

///|
pub fn trace_pic_diagnostics(
  state : PicState,
  dt : Double,
  steps : Int,
) -> Array[PicDiagnostics] {
  let output : Array[PicDiagnostics] = []
  let mut current = state
  output.push(PicDiagnostics::from_state(current))
  for _ in 0.. PicSeriesSummary {
  let first_particles = if values.length() == 0 {
    0
  } else {
    values[0].particles
  }
  let last_particles = if values.length() == 0 {
    0
  } else {
    values[values.length() - 1].particles
  }
  let charges = values.map(fn(value) { value.total_charge })
  let energies = values.map(fn(value) { value.total_energy() })
  let fields = values.map(fn(value) { value.field_energy })
  {
    samples: values.length(),
    particles_initial: first_particles,
    particles_final: last_particles,
    charge_drift: if charges.length() == 0 {
      0.0
    } else {
      relative_change(charges[0], charges[charges.length() - 1])
    },
    energy_drift: if energies.length() == 0 {
      0.0
    } else {
      relative_change(energies[0], energies[energies.length() - 1])
    },
    maximum_field: max_value(fields, 0.0).sqrt(),
  }
}

///|
pub fn pic_series_summary_to_csv(summary : PicSeriesSummary) -> String {
  "samples,particles_initial,particles_final,charge_drift,energy_drift,maximum_field\n\{summary.samples},\{summary.particles_initial},\{summary.particles_final},\{summary.charge_drift},\{summary.energy_drift},\{summary.maximum_field}\n"
}

///|
pub fn pic_diagnostics_times(
  values : ArrayView[PicDiagnostics],
) -> Array[Double] {
  values.map(fn(value) { value.time })
}

///|
pub fn pic_diagnostics_charges(
  values : ArrayView[PicDiagnostics],
) -> Array[Double] {
  values.map(fn(value) { value.total_charge })
}

///|
pub fn pic_diagnostics_energies(
  values : ArrayView[PicDiagnostics],
) -> Array[Double] {
  values.map(fn(value) { value.total_energy() })
}

///|
pub fn pic_diagnostics_to_table(
  values : ArrayView[PicDiagnostics],
) -> NumericTable {
  let table = NumericTable::new(["time", "particles", "charge", "energy"])
  for value in values {
    table.add_row([
      value.time,
      value.particles.to_double(),
      value.total_charge,
      value.total_energy(),
    ])
  }
  table
}

///|
pub fn pic_energy_drift(values : ArrayView[PicDiagnostics]) -> Double {
  let energies = pic_diagnostics_energies(values)
  if energies.length() < 2 {
    0.0
  } else {
    relative_change(energies[0], energies[energies.length() - 1])
  }
}

///|
pub fn pic_charge_drift(values : ArrayView[PicDiagnostics]) -> Double {
  let charges = pic_diagnostics_charges(values)
  if charges.length() < 2 {
    0.0
  } else {
    relative_change(charges[0], charges[charges.length() - 1])
  }
}

///|
pub fn pic_diagnostics_pass(
  values : ArrayView[PicDiagnostics],
  charge_tolerance : Double,
  energy_tolerance : Double,
) -> Bool {
  pic_charge_drift(values) <= charge_tolerance &&
  pic_energy_drift(values) <= energy_tolerance
}