///|
pub(all) struct BandMetrics {
  area : Double
  positive_area : Double
  positive_centroid : Double?
  height : Double
  height_position : Double
  sampled_points : Int32
}

///|
fn segment_moment(x0 : Double, x1 : Double, y0 : Double, y1 : Double) -> Double {
  // Use local coordinates to avoid subtracting large squared coordinates.
  let width = x1 - x0
  let area = width * (y0 + y1) / 2.0
  x0 * area + width * width * (y0 + 2.0 * y1) / 6.0
}

///|
/// Exact integral of the piecewise-linear interpolant, including partial
/// endpoint segments. Optional chord removal is local to this band.
/// Positive centroid uses only signal above zero, independently of signed area.
pub fn integrate_band(
  input : Spectrum,
  left : Double,
  right : Double,
  remove_chord~ : Bool = false,
) -> BandMetrics raise SpectralError {
  finite(left, "band left")
  finite(right, "band right")
  if left >= right {
    raise Invalid("band endpoints must increase")
  }
  let left_value = interpolate(input, left)
  let right_value = interpolate(input, right)
  let axis = [left]
  let values = [left_value]
  for i in 0.. left && input.axis[i] < right {
      axis.push(input.axis[i])
      values.push(input.values[i])
    }
  }
  axis.push(right)
  values.push(right_value)
  if remove_chord {
    for i in 0.. height {
      height = values[i]
      height_position = axis[i]
    }
    if i == 0 {
      continue
    }
    let x0 = axis[i - 1]
    let x1 = axis[i]
    let y0 = values[i - 1]
    let y1 = values[i]
    let piece = (x1 - x0) * (y0 + y1) / 2.0
    let adjusted = piece - compensation
    let next = area + adjusted
    compensation = next - area - adjusted
    area = next
    if y0 >= 0.0 && y1 >= 0.0 {
      positive_area = positive_area + piece
      moment = moment + segment_moment(x0, x1, y0, y1)
    } else if y0 > 0.0 {
      let zero = x0 + (x1 - x0) * y0 / (y0 - y1)
      positive_area = positive_area + (zero - x0) * y0 / 2.0
      moment = moment + segment_moment(x0, zero, y0, 0.0)
    } else if y1 > 0.0 {
      let zero = x0 + (x1 - x0) * -y0 / (y1 - y0)
      positive_area = positive_area + (x1 - zero) * y1 / 2.0
      moment = moment + segment_moment(zero, x1, 0.0, y1)
    }
  }
  ignore(checked_result(area, "band area"))
  ignore(checked_result(positive_area, "positive band area"))
  let positive_centroid = if positive_area > 0.0 {
    Some(checked_result(moment / positive_area, "band centroid"))
  } else {
    None
  }
  {
    area,
    positive_area,
    positive_centroid,
    height,
    height_position,
    sampled_points: axis.length(),
  }
}

///|
/// Divide signed areas; reject a small denominator using a caller threshold.
pub fn band_ratio(
  numerator : BandMetrics,
  denominator : BandMetrics,
  minimum_area~ : Double = 1.0e-12,
) -> Double raise SpectralError {
  finite(minimum_area, "minimum denominator area")
  finite(numerator.area, "numerator area")
  finite(denominator.area, "denominator area")
  if minimum_area <= 0.0 || denominator.area.abs() < minimum_area {
    raise Numerical("band ratio denominator is too small")
  }
  checked_result(numerator.area / denominator.area, "band ratio")
}