///|
pub(all) enum ContactPattern {
  SingleStage
  Crosscurrent(Int)
  Countercurrent(Int)
} derive(Debug, Eq, ToJson)

///|
pub(all) struct SolventDesign {
  pattern : ContactPattern
  solvent_flow : Double
  solvent_to_feed_ratio : Double
  achieved_raffinate : Double
  achieved_recovery : Double
  iterations : Int
  target_raffinate : Double
  target_residual : Double
} derive(Debug, Eq, ToJson)

///|
fn validate_design_inputs(
  carrier_flow : Double,
  feed_concentration : Double,
  distribution : Double,
  target_raffinate : Double,
) -> Unit raise @domain.ExtractionError {
  @domain.require_positive("carrier flow", carrier_flow)
  @domain.require_positive("feed concentration", feed_concentration)
  @domain.require_positive("distribution coefficient", distribution)
  @domain.require_non_negative("target raffinate", target_raffinate)
  if target_raffinate >= feed_concentration {
    raise @domain.ExtractionError::InvalidArgument(
      message="target raffinate must be below feed concentration",
    )
  }
}

///|
/// Closed-form fresh-solvent requirement for one ideal stage.
pub fn solvent_for_single_stage(
  carrier_flow~ : Double,
  feed_concentration~ : Double,
  distribution~ : Double,
  target_raffinate~ : Double,
) -> SolventDesign raise @domain.ExtractionError {
  validate_design_inputs(
    carrier_flow, feed_concentration, distribution, target_raffinate,
  )
  if target_raffinate <= 0.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="a finite single-stage solvent flow cannot reach zero concentration",
    )
  }
  let solvent_flow = carrier_flow *
    (feed_concentration / target_raffinate - 1.0) /
    distribution
  let result = @liquid.single_stage_constant_k(
    feed_flow=carrier_flow,
    feed_concentration~,
    solvent_flow~,
    distribution~,
  )
  {
    pattern: SingleStage,
    solvent_flow,
    solvent_to_feed_ratio: solvent_flow / carrier_flow,
    achieved_raffinate: result.raffinate_concentration,
    achieved_recovery: result.recovery,
    iterations: 0,
    target_raffinate,
    target_residual: result.raffinate_concentration - target_raffinate,
  }
}

///|
fn bisection_upper_bound(
  initial? : Double = 1.0,
  maximum? : Double = 1.0e12,
  meets_target : (Double) -> Bool raise @domain.ExtractionError,
) -> Double raise @domain.ExtractionError {
  @domain.require_positive("initial solvent bound", initial)
  @domain.require_positive("maximum solvent bound", maximum)
  let mut upper = initial
  for expansion in 0..<100 {
    if meets_target(upper) {
      return upper
    }
    upper *= 2.0
    if upper > maximum {
      raise @domain.ExtractionError::ConvergenceFailure(
        iterations=expansion + 1,
        residual=upper,
      )
    }
  }
  raise @domain.ExtractionError::ConvergenceFailure(
    iterations=100,
    residual=upper,
  )
}

///|
pub fn solvent_for_crosscurrent(
  carrier_flow~ : Double,
  feed_concentration~ : Double,
  distribution~ : Double,
  target_raffinate~ : Double,
  stage_count~ : Int,
  tolerance? : Double = 1.0e-9,
) -> SolventDesign raise @domain.ExtractionError {
  validate_design_inputs(
    carrier_flow, feed_concentration, distribution, target_raffinate,
  )
  @domain.require_positive("tolerance", tolerance)
  if stage_count <= 0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="stage_count must be positive",
    )
  }
  let calculate = total_solvent => {
    @liquid.crosscurrent_equal_solvent(
      feed_flow=carrier_flow,
      feed_concentration~,
      total_solvent~,
      stage_count~,
      distribution~,
    )
  }
  let upper_bound = bisection_upper_bound(solvent => {
    calculate(solvent).final_raffinate_concentration <= target_raffinate
  })
  let mut lower = 0.0
  let mut upper = upper_bound
  let mut iterations = 0
  let mut best = calculate(upper)
  for iteration in 0..<120 {
    iterations = iteration + 1
    let middle = lower + (upper - lower) / 2.0
    let candidate = calculate(middle)
    if candidate.final_raffinate_concentration <= target_raffinate {
      upper = middle
      best = candidate
    } else {
      lower = middle
    }
    if upper - lower <= tolerance * (1.0 + upper) {
      break
    }
  }
  {
    pattern: Crosscurrent(stage_count),
    solvent_flow: upper,
    solvent_to_feed_ratio: upper / carrier_flow,
    achieved_raffinate: best.final_raffinate_concentration,
    achieved_recovery: best.overall_recovery,
    iterations,
    target_raffinate,
    target_residual: best.final_raffinate_concentration - target_raffinate,
  }
}

///|
pub fn solvent_for_countercurrent(
  carrier_flow~ : Double,
  feed_concentration~ : Double,
  distribution~ : Double,
  target_raffinate~ : Double,
  stage_count~ : Int,
  solvent_inlet_concentration? : Double = 0.0,
  tolerance? : Double = 1.0e-9,
) -> SolventDesign raise @domain.ExtractionError {
  validate_design_inputs(
    carrier_flow, feed_concentration, distribution, target_raffinate,
  )
  @domain.require_non_negative(
    "solvent inlet concentration", solvent_inlet_concentration,
  )
  @domain.require_positive("tolerance", tolerance)
  if stage_count <= 0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="stage_count must be positive",
    )
  }
  let calculate = solvent_flow => {
    @liquid.countercurrent_constant_k(
      carrier_flow~,
      feed_concentration~,
      solvent_flow~,
      solvent_inlet_concentration~,
      distribution~,
      stage_count~,
    )
  }
  let upper_bound = bisection_upper_bound(solvent => {
    calculate(solvent).raffinate_outlet_concentration <= target_raffinate
  })
  let mut lower = 0.0
  let mut upper = upper_bound
  let mut iterations = 0
  let mut best = calculate(upper)
  for iteration in 0..<120 {
    iterations = iteration + 1
    let middle = lower + (upper - lower) / 2.0
    if middle <= 0.0 {
      lower = middle
      continue
    }
    let candidate = calculate(middle)
    if candidate.raffinate_outlet_concentration <= target_raffinate {
      upper = middle
      best = candidate
    } else {
      lower = middle
    }
    if upper - lower <= tolerance * (1.0 + upper) {
      break
    }
  }
  {
    pattern: Countercurrent(stage_count),
    solvent_flow: upper,
    solvent_to_feed_ratio: upper / carrier_flow,
    achieved_raffinate: best.raffinate_outlet_concentration,
    achieved_recovery: best.recovery,
    iterations,
    target_raffinate,
    target_residual: best.raffinate_outlet_concentration - target_raffinate,
  }
}

///|
/// Compare feasible ideal contact patterns at the same target.
pub fn compare_contact_patterns(
  carrier_flow~ : Double,
  feed_concentration~ : Double,
  distribution~ : Double,
  target_raffinate~ : Double,
  maximum_stages? : Int = 8,
) -> Array[SolventDesign] raise @domain.ExtractionError {
  if maximum_stages <= 0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="maximum_stages must be positive",
    )
  }
  let designs = [
    solvent_for_single_stage(
      carrier_flow~,
      feed_concentration~,
      distribution~,
      target_raffinate~,
    ),
  ]
  for stages in 2..<=maximum_stages {
    designs.push(
      solvent_for_crosscurrent(
        carrier_flow~,
        feed_concentration~,
        distribution~,
        target_raffinate~,
        stage_count=stages,
      ),
    )
    designs.push(
      solvent_for_countercurrent(
        carrier_flow~,
        feed_concentration~,
        distribution~,
        target_raffinate~,
        stage_count=stages,
      ),
    )
  }
  designs.sort_by((left, right) => left.solvent_flow.compare(right.solvent_flow))
  designs
}