///|
/// Trace a streamline with first-order Euler integration.
pub fn trace_streamline(
  velocity : VectorField2D,
  start~ : Point,
  step_size~ : Double,
  max_steps~ : Int,
) -> Array[Point] {
  let points : Array[Point] = [start]
  let mut current = start
  for _ in 0.. (velocity.size.width - 1).to_double() ||
      current.y > (velocity.size.height - 1).to_double() {
      break
    }
  }
  points
}

///|
/// Sample a scalar field along a polyline.
pub fn sample_polyline(
  field : Field2D,
  points : ArrayView[Point],
) -> Array[Double] {
  let values = Array::new()
  for point in points {
    values.push(field.sample(x=point.x, y=point.y))
  }
  values
}

///|
/// Return a normalized histogram with a fixed number of bins.
pub fn histogram(
  values : ArrayView[Double],
  bins~ : Int,
  low~ : Double,
  high~ : Double,
) -> Array[Int] {
  let result = Array::make(bins.max(0), 0)
  if high > low && bins > 0 {
    for value in values {
      let t = clamp_double(
        (value - low) / (high - low),
        low=0.0,
        high=0.999999999,
      )
      let index = (t * bins.to_double()).to_int().clamp(min=0, max=bins - 1)
      result[index] += 1
    }
  }
  result
}

///|
/// Compute a central-difference gradient at one scalar cell.
pub fn gradient_at(field : Field2D, x~ : Int, y~ : Int) -> Point {
  Point::new(
    x=(field.get(x + 1, y) - field.get(x - 1, y)) * 0.5,
    y=(field.get(x, y + 1) - field.get(x, y - 1)) * 0.5,
  )
}

///|
/// Compute the scalar Laplacian at one cell.
pub fn laplacian_at(field : Field2D, x~ : Int, y~ : Int) -> Double {
  field.get(x - 1, y) +
  field.get(x + 1, y) +
  field.get(x, y - 1) +
  field.get(x, y + 1) -
  4.0 * field.get(x, y)
}

///|
/// Return the maximum speed in a vector field.
pub fn max_speed(velocity : VectorField2D) -> Double {
  let mut maximum = 0.0
  for y in 0..