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

///|
pub fn particle_speed(particle : Particle) -> Double {
  particle.v.abs()
}

///|
pub fn particle_speeds(particles : ArrayView[Particle]) -> Array[Double] {
  particles.map(fn(particle) { particle_speed(particle) })
}

///|
pub fn particle_position_bounds(
  particles : ArrayView[Particle],
) -> (Double, Double) {
  (
    min_value(particles.map(fn(particle) { particle.x }), 0.0),
    max_value(particles.map(fn(particle) { particle.x }), 0.0),
  )
}

///|
pub fn particle_velocity_bounds(
  particles : ArrayView[Particle],
) -> (Double, Double) {
  (
    min_value(particles.map(fn(particle) { particle.v }), 0.0),
    max_value(particles.map(fn(particle) { particle.v }), 0.0),
  )
}

///|
pub fn particles_in_cell(
  grid : Grid1D,
  particles : ArrayView[Particle],
  cell : Int,
) -> Array[Particle] {
  let output : Array[Particle] = []
  let wanted = modulo_index(cell, grid.cells)
  for particle in particles {
    if nearest_index(grid, particle.x) == wanted {
      output.push(particle)
    }
  }
  output
}

///|
pub fn particles_with_speed_limit(
  particles : ArrayView[Particle],
  limit : Double,
) -> Array[Particle] {
  particles.filter(fn(particle) { particle_speed(particle) <= limit })
}

///|
pub fn particles_with_position_range(
  particles : ArrayView[Particle],
  low : Double,
  high : Double,
) -> Array[Particle] {
  particles.filter(fn(particle) { particle.x >= low && particle.x <= high })
}

///|
pub fn rescale_particle_weights(
  particles : ArrayView[Particle],
  factor : Double,
) -> Array[Particle] {
  particles.map(fn(particle) {
    Particle::new(
      x=particle.x,
      v=particle.v,
      weight=particle.weight * factor,
      charge=particle.charge,
      mass=particle.mass,
    )
  })
}

///|
pub fn translate_particles(
  particles : ArrayView[Particle],
  shift : Double,
) -> Array[Particle] {
  particles.map(fn(particle) {
    Particle::new(
      x=particle.x + shift,
      v=particle.v,
      weight=particle.weight,
      charge=particle.charge,
      mass=particle.mass,
    )
  })
}

///|
pub fn accelerate_particles(
  particles : ArrayView[Particle],
  acceleration : Double,
  dt : Double,
) -> Array[Particle] {
  particles.map(fn(particle) {
    let velocity = particle.v + acceleration * dt
    Particle::new(
      x=particle.x + velocity * dt,
      v=velocity,
      weight=particle.weight,
      charge=particle.charge,
      mass=particle.mass,
    )
  })
}

///|
pub fn particle_energy_series(
  states : ArrayView[Array[Particle]],
) -> Array[Double] {
  states.map(fn(particles) { kinetic_energy(particles) })
}

///|
pub fn particle_charge_series(
  states : ArrayView[Array[Particle]],
) -> Array[Double] {
  states.map(fn(particles) { total_particle_charge(particles) })
}

///|
pub fn particle_momentum_series(
  states : ArrayView[Array[Particle]],
) -> Array[Double] {
  states.map(fn(particles) { particle_total_momentum(particles) })
}

///|
pub fn pic_state_to_snapshot(state : PicState) -> String {
  let output = StringBuilder()
  output.write_string("time=\{state.time}\n")
  output.write_string(particles_to_csv(state.particles))
  output.write_string(field_to_csv(state.field))
  output.to_string()
}

///|
pub fn pic_state_charge_error(state : PicState) -> Double {
  (total_charge(state.grid, state.field.charge_density) -
  total_particle_charge(state.particles)).abs()
}

///|
pub fn pic_state_energy(state : PicState) -> Double {
  kinetic_energy(state.particles) + field_energy(state.field)
}

///|
pub fn pic_state_is_finite(state : PicState) -> Bool {
  state.particles.length() >= 0 &&
  state.field.electric.length() == state.grid.cells
}