///|
pub(all) struct ExtractionStep {
  stage : Int
  raffinate_out : Double
  extract_out : Double
  raffinate_in : Double
  extract_in : Double
  horizontal_span : Double
} derive(Debug, Eq, ToJson)

///|
pub(all) struct StageCountResult {
  steps : Array[ExtractionStep]
  theoretical_stages : Double
  complete_stages : Int
  feed_concentration_reached : Double
  target_raffinate : Double
  balance_residual : Double
} derive(Debug, Eq, ToJson)

///|
/// McCabe-Thiele style stepping for countercurrent liquid extraction.
///
/// Starting at the raffinate product end, each step first reaches the
/// equilibrium curve and then follows a stage material-balance line toward
/// the feed. A fractional final stage is reported by linear interpolation of
/// the last horizontal concentration span.
pub fn mccabe_thiele_stages(
  equilibrium~ : @interpolation.EquilibriumTable,
  carrier_flow~ : Double,
  solvent_flow~ : Double,
  feed_concentration~ : Double,
  target_raffinate~ : Double,
  solvent_inlet_concentration? : Double = 0.0,
  maximum_stages? : Int = 200,
  tolerance? : Double = 1.0e-10,
) -> StageCountResult raise @domain.ExtractionError {
  @domain.require_positive("carrier flow", carrier_flow)
  @domain.require_positive("solvent flow", solvent_flow)
  @domain.require_non_negative("feed concentration", feed_concentration)
  @domain.require_non_negative("target raffinate", target_raffinate)
  @domain.require_non_negative(
    "solvent inlet concentration", solvent_inlet_concentration,
  )
  @domain.require_positive("tolerance", tolerance)
  if target_raffinate >= feed_concentration {
    raise @domain.ExtractionError::InvalidArgument(
      message="target raffinate must be below feed concentration",
    )
  }
  if maximum_stages <= 0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="maximum_stages must be positive",
    )
  }
  let steps : Array[ExtractionStep] = []
  let mut current_x = target_raffinate
  let mut entering_y = solvent_inlet_concentration
  let mut theoretical_stages = 0.0
  let mut complete_stages = 0
  let mut residual = 0.0
  for stage_index in 0..= feed_concentration - tolerance {
      let fraction = ((feed_concentration - current_x) / span).clamp(
        min=0.0,
        max=1.0,
      )
      theoretical_stages += fraction
      if fraction >= 1.0 - tolerance {
        complete_stages += 1
      }
      current_x = feed_concentration
      let predicted_inlet = carrier_flow * feed_concentration +
        solvent_flow * solvent_inlet_concentration
      let predicted_outlet = carrier_flow * target_raffinate +
        solvent_flow * leaving_y
      residual = predicted_inlet - predicted_outlet
      return {
        steps,
        theoretical_stages,
        complete_stages,
        feed_concentration_reached: current_x,
        target_raffinate,
        balance_residual: residual,
      }
    }
    theoretical_stages += 1.0
    complete_stages += 1
    current_x = next_x
    entering_y = leaving_y
  }
  raise @domain.ExtractionError::ConvergenceFailure(
    iterations=maximum_stages,
    residual=feed_concentration - current_x,
  )
}

///|
/// Build points on the overall operating line for plotting or tabulation.
pub fn operating_line(
  carrier_flow~ : Double,
  solvent_flow~ : Double,
  raffinate_product~ : Double,
  solvent_inlet~ : Double,
  xs : ArrayView[Double],
) -> Array[@domain.EquilibriumPoint] raise @domain.ExtractionError {
  @domain.require_positive("carrier flow", carrier_flow)
  @domain.require_positive("solvent flow", solvent_flow)
  @domain.require_non_negative("raffinate product", raffinate_product)
  @domain.require_non_negative("solvent inlet", solvent_inlet)
  [
    for x in xs => {
      @domain.require_non_negative("operating-line x", x)
      @domain.EquilibriumPoint::new(
        x~,
        y=solvent_inlet + carrier_flow * (x - raffinate_product) / solvent_flow,
      )
    }
  ]
}