///|
/// Make an inclusive regular grid; caller explicitly chooses axis units.
pub fn regular_axis(
  start : Double,
  stop : Double,
  count : Int32,
) -> Array[Double] raise SpectralError {
  finite(start, "grid start")
  finite(stop, "grid stop")
  if count < 3 || count > 65536 {
    raise Budget("grid count must lie in [3,65536]")
  }
  if stop <= start {
    raise Invalid("grid endpoints must increase")
  }
  if start.abs() > 1.0e100 || stop.abs() > 1.0e100 {
    raise Invalid("grid magnitude exceeds 1e100")
  }
  let step = (stop - start) / (count - 1).to_double()
  let result = Array::makei(count, fn(i) { start + step * i.to_double() })
  result[count - 1] = stop
  for i in 1.. Int32 {
  let mut low = 0
  let mut high = axis.length() - 1
  while high - low > 1 {
    let mid = low + (high - low) / 2
    if axis[mid] <= x {
      low = mid
    } else {
      high = mid
    }
  }
  low
}

///|
/// Linear interpolation within the measured domain; extrapolation is refused.
pub fn interpolate(input : Spectrum, x : Double) -> Double raise SpectralError {
  finite(x, "interpolation coordinate")
  if x < input.axis[0] || x > input.axis[input.length() - 1] {
    raise Invalid("interpolation would extrapolate")
  }
  let i = lower_segment(input.axis, x)
  let fraction = (x - input.axis[i]) / (input.axis[i + 1] - input.axis[i])
  checked_result(
    input.values[i] * (1.0 - fraction) + input.values[i + 1] * fraction,
    "interpolation",
  )
}

///|
/// Two-pointer resampling in O(n+m), with owned axis and intensity arrays.
pub fn resample(
  input : Spectrum,
  axis : Array[Double],
) -> Spectrum raise SpectralError {
  if axis.length() < 3 || axis.length() > 65536 {
    raise Budget("resampling needs 3..65536 coordinates")
  }
  for i in 0.. 0 && axis[i] <= axis[i - 1] {
      raise Invalid("resampling axis must increase")
    }
  }
  if axis[0] < input.axis[0] ||
    axis[axis.length() - 1] > input.axis[input.length() - 1] {
    raise Invalid("resampling would extrapolate")
  }
  let output = Array::make(axis.length(), 0.0)
  let mut segment = 0
  for i in 0.. Spectrum raise SpectralError {
  finite(left, "crop left")
  finite(right, "crop right")
  if left >= right {
    raise Invalid("crop endpoints must increase")
  }
  let axis = [left]
  let values = [interpolate(input, left)]
  for i in 0.. left && input.axis[i] < right {
      axis.push(input.axis[i])
      values.push(input.values[i])
    }
  }
  axis.push(right)
  values.push(interpolate(input, right))
  spectrum(axis, values)
}

///|
/// Relative spacing tolerance is measured against mean spacing.
pub fn is_uniform_axis(
  input : Spectrum,
  tolerance~ : Double = 1.0e-6,
) -> Bool raise SpectralError {
  finite(tolerance, "spacing tolerance")
  if tolerance < 0.0 || tolerance > 1.0 {
    raise Invalid("spacing tolerance outside [0,1]")
  }
  let step = (input.axis[input.length() - 1] - input.axis[0]) /
    (input.length() - 1).to_double()
  for i in 1.. tolerance * step {
      return false
    }
  }
  true
}