///|
pub(all) struct EfficiencyStep {
  stage : Int
  efficiency : Double
  raffinate_out : Double
  equilibrium_extract : Double
  actual_extract_out : Double
  raffinate_in : Double
  extract_in : Double
  equilibrium_approach : Double
} derive(Debug, Eq, ToJson)

///|
pub(all) struct EfficiencyStageResult {
  steps : Array[EfficiencyStep]
  actual_stages : Double
  complete_stages : Int
  target_raffinate : Double
  feed_concentration_reached : Double
  average_efficiency : Double
} derive(Debug, Eq, ToJson)

///|
fn validate_efficiency(
  value : Double,
  stage : Int,
) -> Unit raise @domain.ExtractionError {
  @domain.require_fraction("stage \{stage} efficiency", value)
  if value <= 0.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="stage efficiency must be greater than zero",
    )
  }
}

///|
/// Step through a countercurrent extractor with per-stage Murphree-style
/// extract-phase efficiencies.
///
/// For each contact, `y_actual = y_in + eta (y_equilibrium - y_in)`. The
/// resulting solute transfer determines the raffinate concentration entering
/// the next contact toward the feed end.
pub fn stages_with_efficiency_profile(
  equilibrium~ : @interpolation.EquilibriumTable,
  carrier_flow~ : Double,
  solvent_flow~ : Double,
  feed_concentration~ : Double,
  target_raffinate~ : Double,
  efficiencies : ArrayView[Double],
  solvent_inlet_concentration? : Double = 0.0,
  tolerance? : Double = 1.0e-10,
) -> EfficiencyStageResult raise @domain.ExtractionError {
  @domain.require_positive("carrier flow", carrier_flow)
  @domain.require_positive("solvent flow", solvent_flow)
  @domain.require_positive("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 efficiencies.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  let steps : Array[EfficiencyStep] = []
  let mut current_x = target_raffinate
  let mut entering_y = solvent_inlet_concentration
  let mut actual_stages = 0.0
  let mut complete_stages = 0
  let mut efficiency_sum = 0.0
  for index, efficiency in efficiencies {
    validate_efficiency(efficiency, index + 1)
    efficiency_sum += efficiency
    let equilibrium_y = equilibrium.y_at(current_x)
    let driving_force = equilibrium_y - entering_y
    if driving_force <= tolerance {
      raise @domain.ExtractionError::ConvergenceFailure(
        iterations=index,
        residual=feed_concentration - current_x,
      )
    }
    let actual_y = entering_y + efficiency * driving_force
    let next_x = current_x +
      solvent_flow * (actual_y - entering_y) / carrier_flow
    let span = next_x - current_x
    if span <= tolerance {
      raise @domain.ExtractionError::ConvergenceFailure(
        iterations=index,
        residual=feed_concentration - current_x,
      )
    }
    steps.push({
      stage: index + 1,
      efficiency,
      raffinate_out: current_x,
      equilibrium_extract: equilibrium_y,
      actual_extract_out: actual_y,
      raffinate_in: next_x,
      extract_in: entering_y,
      equilibrium_approach: if driving_force > 0.0 {
        (actual_y - entering_y) / driving_force
      } else {
        0.0
      },
    })
    if next_x >= feed_concentration - tolerance {
      let fraction = ((feed_concentration - current_x) / span).clamp(
        min=0.0,
        max=1.0,
      )
      actual_stages += fraction
      if fraction >= 1.0 - tolerance {
        complete_stages += 1
      }
      return {
        steps,
        actual_stages,
        complete_stages,
        target_raffinate,
        feed_concentration_reached: feed_concentration,
        average_efficiency: efficiency_sum / steps.length().to_double(),
      }
    }
    actual_stages += 1.0
    complete_stages += 1
    current_x = next_x
    entering_y = actual_y
  }
  raise @domain.ExtractionError::ConvergenceFailure(
    iterations=efficiencies.length(),
    residual=feed_concentration - current_x,
  )
}

///|
/// Uniform-efficiency wrapper that allocates at most `maximum_stages` steps.
pub fn stages_with_uniform_efficiency(
  equilibrium~ : @interpolation.EquilibriumTable,
  carrier_flow~ : Double,
  solvent_flow~ : Double,
  feed_concentration~ : Double,
  target_raffinate~ : Double,
  efficiency~ : Double,
  solvent_inlet_concentration? : Double = 0.0,
  maximum_stages? : Int = 500,
  tolerance? : Double = 1.0e-10,
) -> EfficiencyStageResult raise @domain.ExtractionError {
  validate_efficiency(efficiency, 1)
  if maximum_stages <= 0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="maximum_stages must be positive",
    )
  }
  stages_with_efficiency_profile(
    equilibrium~,
    carrier_flow~,
    solvent_flow~,
    feed_concentration~,
    target_raffinate~,
    Array::make(maximum_stages, efficiency),
    solvent_inlet_concentration~,
    tolerance~,
  )
}

///|
/// Convert ideal theoretical stages to a shortcut actual-stage estimate.
pub fn shortcut_actual_stages(
  theoretical_stages : Double,
  overall_efficiency : Double,
) -> Double raise @domain.ExtractionError {
  @domain.require_non_negative("theoretical stages", theoretical_stages)
  validate_efficiency(overall_efficiency, 1)
  theoretical_stages / overall_efficiency
}