///|
/// Continuous vector fields for particles, brush strokes, and procedural
/// camera offsets. A field owns its evaluator so integration stays backend
/// independent.
pub enum VectorFieldMode2 {
  UniformField(Point2)
  RadialField(Point2, Double)
  VortexField(Point2, Double)
  NoiseField(NoiseSource2, Double)
} derive(Debug)

///|
pub fn vector_field_modes2() -> Array[VectorFieldMode2] {
  [
    UniformField(point2(x=1.0, y=0.0)),
    RadialField(point2(x=0.0, y=0.0), 1.0),
    VortexField(point2(x=0.0, y=0.0), 1.0),
    NoiseField(noise_source2(1), 1.0),
  ]
}

///|
pub fn vector_field_value(
  mode : VectorFieldMode2,
  position : Point2,
  time : Double,
) -> Point2 {
  match mode {
    UniformField(value) => value
    RadialField(center, strength) => {
      let delta = point2(
        x=position.x() - center.x(),
        y=position.y() - center.y(),
      )
      let distance = point2_length(delta).max(0.000000000001)
      point2(
        x=delta.x() * strength / distance,
        y=delta.y() * strength / distance,
      )
    }
    VortexField(center, strength) => {
      let delta = point2(
        x=position.x() - center.x(),
        y=position.y() - center.y(),
      )
      let distance = point2_length(delta).max(0.000000000001)
      point2(
        x=-delta.y() * strength / distance,
        y=delta.x() * strength / distance,
      )
    }
    NoiseField(source, strength) => {
      let gradient = source.gradient(position.x() + time, position.y())
      point2(x=gradient.y() * strength, y=-gradient.x() * strength)
    }
  }
}

///|
pub struct VectorField2 {
  mode : VectorFieldMode2
  scale : Double
  falloff : Double
  bounds : Bounds2D?
} derive(Debug)

///|
pub fn vector_field2(mode : VectorFieldMode2) -> VectorField2 {
  { mode, scale: 1.0, falloff: 0.0, bounds: None }
}

///|
pub fn VectorField2::with_scale(
  self : VectorField2,
  scale : Double,
) -> VectorField2 {
  { mode: self.mode, scale, falloff: self.falloff, bounds: self.bounds }
}

///|
pub fn VectorField2::with_falloff(
  self : VectorField2,
  falloff : Double,
) -> VectorField2 {
  {
    mode: self.mode,
    scale: self.scale,
    falloff: falloff.max(0.0),
    bounds: self.bounds,
  }
}

///|
pub fn VectorField2::with_bounds(
  self : VectorField2,
  bounds : Bounds2D,
) -> VectorField2 {
  {
    mode: self.mode,
    scale: self.scale,
    falloff: self.falloff,
    bounds: Some(bounds),
  }
}

///|
pub fn VectorField2::mode(self : VectorField2) -> VectorFieldMode2 {
  self.mode
}

///|
pub fn VectorField2::sample(
  self : VectorField2,
  position : Point2,
  time : Double,
) -> Point2 {
  let value = vector_field_value(self.mode, position, time)
  let scale = if self.falloff == 0.0 {
    self.scale
  } else {
    self.scale / (1.0 + self.falloff * point2_length(value))
  }
  point2(x=value.x() * scale, y=value.y() * scale)
}

///|
pub struct FieldParticle2 {
  position : Point2
  velocity : Point2
  age : Double
  lifetime : Double
} derive(Debug)

///|
pub fn field_particle2(
  position : Point2,
  velocity : Point2,
  lifetime : Double,
) -> FieldParticle2 {
  { position, velocity, age: 0.0, lifetime: lifetime.max(0.000001) }
}

///|
pub fn FieldParticle2::position(self : FieldParticle2) -> Point2 {
  self.position
}

///|
pub fn FieldParticle2::velocity(self : FieldParticle2) -> Point2 {
  self.velocity
}

///|
pub fn FieldParticle2::age(self : FieldParticle2) -> Double {
  self.age
}

///|
pub fn FieldParticle2::lifetime(self : FieldParticle2) -> Double {
  self.lifetime
}

///|
pub fn FieldParticle2::alive(self : FieldParticle2) -> Bool {
  self.age < self.lifetime
}

///|
pub fn FieldParticle2::progress(self : FieldParticle2) -> Double {
  clamp01(self.age / self.lifetime)
}

///|
pub fn VectorField2::step(
  self : VectorField2,
  particle : FieldParticle2,
  time : Double,
  dt : Double,
  drag : Double,
) -> FieldParticle2 {
  let delta = dt.max(0.0)
  let field = self.sample(particle.position, time)
  let velocity = point2(
    x=particle.velocity.x() + field.x() * delta,
    y=particle.velocity.y() + field.y() * delta,
  )
  let damp = @math.pow(1.0 - clamp01(drag), delta)
  let damped = point2(x=velocity.x() * damp, y=velocity.y() * damp)
  {
    position: point2(
      x=particle.position.x() + damped.x() * delta,
      y=particle.position.y() + damped.y() * delta,
    ),
    velocity: damped,
    age: particle.age + delta,
    lifetime: particle.lifetime,
  }
}

///|
pub fn VectorField2::integrate(
  self : VectorField2,
  particle : FieldParticle2,
  start : Double,
  end : Double,
  steps : Int,
  drag : Double,
) -> Array[Point2] {
  let result : Array[Point2] = []
  let actual = steps.max(1)
  let dt = (end - start) / actual.to_double()
  let mut current = particle
  result.push(current.position)
  for index in 0.. Array[Point2] {
  let result : Array[Point2] = []
  for point in points {
    let velocity = field.sample(point, time)
    result.push(
      point2(x=point.x() + velocity.x() * dt, y=point.y() + velocity.y() * dt),
    )
  }
  result
}

///|
pub fn field_divergence2(
  field : VectorField2,
  position : Point2,
  time : Double,
  epsilon? : Double = 0.001,
) -> Double {
  let step = epsilon.max(0.000000000001)
  let right = field.sample(point2(x=position.x() + step, y=position.y()), time)
  let left = field.sample(point2(x=position.x() - step, y=position.y()), time)
  let top = field.sample(point2(x=position.x(), y=position.y() + step), time)
  let bottom = field.sample(point2(x=position.x(), y=position.y() - step), time)
  (right.x() - left.x() + (top.y() - bottom.y())) / (2.0 * step)
}

///|
pub fn field_curl2(
  field : VectorField2,
  position : Point2,
  time : Double,
  epsilon? : Double = 0.001,
) -> Double {
  let step = epsilon.max(0.000000000001)
  let right = field.sample(point2(x=position.x() + step, y=position.y()), time)
  let left = field.sample(point2(x=position.x() - step, y=position.y()), time)
  let top = field.sample(point2(x=position.x(), y=position.y() + step), time)
  let bottom = field.sample(point2(x=position.x(), y=position.y() - step), time)
  (top.x() - bottom.x() - (right.y() - left.y())) / (2.0 * step)
}