///|
/// Material-balance result for one equilibrium contact.
pub(all) struct SingleStageResult {
  raffinate_concentration : Double
  extract_concentration : Double
  raffinate_solute : Double
  extract_solute : Double
  recovery : Double
  balance_residual : Double
} derive(Debug, Eq, ToJson)

///|
fn make_single_stage_result(
  feed_flow : Double,
  feed_concentration : Double,
  solvent_flow : Double,
  solvent_concentration : Double,
  raffinate_concentration : Double,
  extract_concentration : Double,
) -> SingleStageResult {
  let inlet = feed_flow * feed_concentration +
    solvent_flow * solvent_concentration
  let raffinate_solute = feed_flow * raffinate_concentration
  let extract_solute = solvent_flow * extract_concentration
  let residual = inlet - raffinate_solute - extract_solute
  let recovery = if inlet > 0.0 { extract_solute / inlet } else { 0.0 }
  {
    raffinate_concentration,
    extract_concentration,
    raffinate_solute,
    extract_solute,
    recovery,
    balance_residual: residual,
  }
}

///|
/// Solve one ideal equilibrium stage with constant distribution coefficient.
/// Concentrations and flow rates must use a consistent solute-free basis.
pub fn single_stage_constant_k(
  feed_flow~ : Double,
  feed_concentration~ : Double,
  solvent_flow~ : Double,
  solvent_concentration? : Double = 0.0,
  distribution~ : Double,
) -> SingleStageResult raise @domain.ExtractionError {
  @domain.require_positive("feed flow", feed_flow)
  @domain.require_positive("solvent flow", solvent_flow)
  @domain.require_non_negative("feed concentration", feed_concentration)
  @domain.require_non_negative("solvent concentration", solvent_concentration)
  @domain.require_non_negative("distribution coefficient", distribution)
  let denominator = feed_flow + distribution * solvent_flow
  if denominator <= 0.0 {
    raise @domain.ExtractionError::SingularBalance(
      context="single-stage denominator is zero",
    )
  }
  let x = (
      feed_flow * feed_concentration + solvent_flow * solvent_concentration
    ) /
    denominator
  make_single_stage_result(
    feed_flow,
    feed_concentration,
    solvent_flow,
    solvent_concentration,
    x,
    distribution * x,
  )
}

///|
/// Solve one stage against a tabulated equilibrium relation using bisection.
pub fn single_stage_table(
  feed_flow~ : Double,
  feed_concentration~ : Double,
  solvent_flow~ : Double,
  solvent_concentration? : Double = 0.0,
  equilibrium~ : @interpolation.EquilibriumTable,
  tolerance? : Double = 1.0e-10,
  max_iterations? : Int = 200,
) -> SingleStageResult raise @domain.ExtractionError {
  @domain.require_positive("feed flow", feed_flow)
  @domain.require_positive("solvent flow", solvent_flow)
  @domain.require_non_negative("feed concentration", feed_concentration)
  @domain.require_non_negative("solvent concentration", solvent_concentration)
  @domain.require_positive("tolerance", tolerance)
  if max_iterations <= 0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="max_iterations must be positive",
    )
  }
  let inlet = feed_flow * feed_concentration +
    solvent_flow * solvent_concentration
  let residual = x => feed_flow * x + solvent_flow * equilibrium.y_at(x) - inlet
  let mut lower = equilibrium.relation.lower()
  let mut upper = equilibrium.relation.upper()
  let lower_residual = residual(lower)
  let upper_residual = residual(upper)
  if lower_residual > 0.0 || upper_residual < 0.0 {
    raise @domain.ExtractionError::OutOfRange(
      name="stage solution",
      value=inlet,
      lower=feed_flow * lower + solvent_flow * equilibrium.y_at(lower),
      upper=feed_flow * upper + solvent_flow * equilibrium.y_at(upper),
    )
  }
  let mut final_residual = upper_residual
  let mut iterations = 0
  for iteration in 0.. tolerance * 10.0 && upper - lower > tolerance {
    raise @domain.ExtractionError::ConvergenceFailure(
      iterations~,
      residual=final_residual,
    )
  }
  let x = lower + (upper - lower) / 2.0
  let y = equilibrium.y_at(x)
  make_single_stage_result(
    feed_flow, feed_concentration, solvent_flow, solvent_concentration, x, y,
  )
}