///|
pub(all) struct Vector3 {
  x : Double
  y : Double
  z : Double
} derive(Eq, Debug, ToJson)

///|
pub fn Vector3::new(x : Double, y : Double, z : Double) -> Vector3 {
  { x, y, z }
}

///|
pub fn Vector3::zero() -> Vector3 {
  { x: 0.0, y: 0.0, z: 0.0 }
}

///|
pub fn Vector3::add(a : Vector3, b : Vector3) -> Vector3 {
  { x: a.x + b.x, y: a.y + b.y, z: a.z + b.z }
}

///|
pub fn Vector3::sub(a : Vector3, b : Vector3) -> Vector3 {
  { x: a.x - b.x, y: a.y - b.y, z: a.z - b.z }
}

///|
pub fn Vector3::scale(a : Vector3, factor : Double) -> Vector3 {
  { x: a.x * factor, y: a.y * factor, z: a.z * factor }
}

///|
pub fn Vector3::dot(a : Vector3, b : Vector3) -> Double {
  a.x * b.x + a.y * b.y + a.z * b.z
}

///|
pub fn Vector3::cross(a : Vector3, b : Vector3) -> Vector3 {
  {
    x: a.y * b.z - a.z * b.y,
    y: a.z * b.x - a.x * b.z,
    z: a.x * b.y - a.y * b.x,
  }
}

///|
pub fn Vector3::norm(a : Vector3) -> Double {
  a.dot(a).sqrt()
}

///|
pub fn Vector3::distance(a : Vector3, b : Vector3) -> Double {
  a.sub(b).norm()
}

///|
pub fn Vector3::normalized(a : Vector3) -> Vector3 {
  let magnitude = a.norm()
  if magnitude == 0.0 {
    Vector3::zero()
  } else {
    a.scale(1.0 / magnitude)
  }
}

///|
pub fn Vector3::component_min(a : Vector3, b : Vector3) -> Vector3 {
  { x: a.x.min(b.x), y: a.y.min(b.y), z: a.z.min(b.z) }
}

///|
pub fn Vector3::component_max(a : Vector3, b : Vector3) -> Vector3 {
  { x: a.x.max(b.x), y: a.y.max(b.y), z: a.z.max(b.z) }
}

///|
pub fn boris_push(
  velocity : Vector3,
  electric : Vector3,
  magnetic : Vector3,
  dt : Double,
  charge : Double,
  mass : Double,
) -> Vector3 {
  if mass == 0.0 {
    velocity
  } else {
    let half = charge * dt / (2.0 * mass)
    let v_minus = velocity.add(electric.scale(half))
    let t = magnetic.scale(half)
    let s = t.scale(2.0 / (1.0 + t.dot(t)))
    let v_prime = v_minus.add(v_minus.cross(t))
    let v_plus = v_minus.add(v_prime.cross(s))
    v_plus.add(electric.scale(half))
  }
}

///|
pub fn boris_position(
  position : Vector3,
  velocity : Vector3,
  dt : Double,
) -> Vector3 {
  position.add(velocity.scale(dt))
}

///|
pub(all) struct VerletState {
  position : Double
  velocity : Double
  acceleration : Double
} derive(Debug, ToJson)

///|
pub fn VerletState::new(
  position : Double,
  velocity : Double,
  acceleration : Double,
) -> VerletState {
  { position, velocity, acceleration }
}

///|
pub fn advance_verlet(
  state : VerletState,
  dt : Double,
  next_acceleration : Double,
) -> VerletState {
  let next_position = state.position +
    state.velocity * dt +
    0.5 * state.acceleration * dt * dt
  let next_velocity = state.velocity +
    0.5 * (state.acceleration + next_acceleration) * dt
  {
    position: next_position,
    velocity: next_velocity,
    acceleration: next_acceleration,
  }
}

///|
pub fn advance_leapfrog(
  position : Double,
  velocity : Double,
  acceleration : Double,
  dt : Double,
) -> (Double, Double) {
  let half_velocity = velocity + 0.5 * acceleration * dt
  (position + half_velocity * dt, half_velocity + 0.5 * acceleration * dt)
}

///|
pub fn acceleration_from_field(
  particle : Particle,
  electric : Double,
) -> Double {
  if particle.mass == 0.0 {
    0.0
  } else {
    particle.charge * electric / particle.mass
  }
}

///|
pub fn push_particle_unbounded(
  particle : Particle,
  electric : Double,
  dt : Double,
) -> Particle {
  let acceleration = acceleration_from_field(particle, electric)
  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 push_particles(
  particles : ArrayView[Particle],
  electric : ArrayView[Double],
  dt : Double,
) -> Array[Particle] {
  particles.map(fn(particle) {
    let index = if electric.length() == 0 {
      0
    } else {
      modulo_index(particles.length(), electric.length())
    }
    push_particle_unbounded(
      particle,
      if electric.length() == 0 {
        0.0
      } else {
        electric[index]
      },
      dt,
    )
  })
}

///|
pub fn integrate_position(
  position : Double,
  velocity : Double,
  acceleration : Double,
  dt : Double,
) -> Double {
  position + velocity * dt + 0.5 * acceleration * dt * dt
}

///|
pub fn integrate_velocity(
  velocity : Double,
  acceleration : Double,
  dt : Double,
) -> Double {
  velocity + acceleration * dt
}

///|
pub fn kinetic_energy_vector(mass : Double, velocity : Vector3) -> Double {
  0.5 * mass * velocity.dot(velocity)
}

///|
pub fn momentum_vector(mass : Double, velocity : Vector3) -> Vector3 {
  velocity.scale(mass)
}

///|
pub fn angular_frequency(
  charge : Double,
  magnetic : Vector3,
  mass : Double,
) -> Double {
  if mass == 0.0 {
    0.0
  } else {
    charge * magnetic.norm() / mass
  }
}

///|
pub fn cyclotron_period(
  charge : Double,
  magnetic : Vector3,
  mass : Double,
) -> Double {
  let frequency = angular_frequency(charge, magnetic, mass).abs()
  if frequency == 0.0 {
    0.0
  } else {
    2.0 * pi / frequency
  }
}