///|
pub(all) struct Particle {
  mut x : Double
  mut v : Double
  weight : Double
  charge : Double
  mass : Double
} derive(Debug, ToJson)

///|
pub fn Particle::new(
  x~ : Double,
  v~ : Double,
  weight? : Double = 1.0,
  charge? : Double = -elementary_charge,
  mass? : Double = electron_mass,
) -> Particle {
  { x, v, weight, charge, mass }
}

///|
pub(all) struct PicState {
  grid : Grid1D
  particles : Array[Particle]
  field : Field1D
  time : Double
} derive(Debug, ToJson)

///|
pub fn PicState::new(grid : Grid1D, particles : Array[Particle]) -> PicState {
  let rho = deposit_charge(grid, particles)
  { grid, particles, field: solve_periodic_field(grid, rho), time: 0.0 }
}

///|
pub fn deposit_charge(
  grid : Grid1D,
  particles : ArrayView[Particle],
) -> Array[Double] {
  let rho = 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 q = p.charge * p.weight / grid.dx
    rho[left] = rho[left] + q * (1.0 - frac)
    rho[right] = rho[right] + q * frac
  }
  rho
}

///|
pub fn PicState::step(state : PicState, dt : Double) -> PicState {
  let grid = state.grid
  let rho = deposit_charge(grid, state.particles)
  let field = solve_periodic_field(grid, rho)
  let next_particles = state.particles.map(fn(p) {
    let e = sample_linear(grid, field.electric, p.x)
    let accel = p.charge * e / p.mass
    let v = p.v + accel * dt
    Particle::new(
      x=grid.wrap(p.x + v * dt),
      v~,
      weight=p.weight,
      charge=p.charge,
      mass=p.mass,
    )
  })
  { grid, particles: next_particles, field, time: state.time + dt }
}

///|
pub fn PicState::run(state : PicState, dt : Double, steps : Int) -> PicState {
  for current = state, _i = 0; _i < steps; {
    continue current.step(dt), _i + 1
  } nobreak {
    current
  }
}

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

///|
pub fn field_energy(field : Field1D) -> Double {
  0.5 *
  vacuum_permittivity *
  field.grid.dx *
  field.electric.fold(init=0.0, fn(acc, e) { acc + e * e })
}