///|
pub fn flash_isothermal_ideal(
  components : Array[Component],
  feed : Array[Double],
  temperature_k : Double,
  pressure_bar : Double,
) -> FlashResult raise VleError {
  assert_pressure(pressure_bar)
  assert_same_length(components.length(), feed.length())
  let z = normalize(feed)
  let ks = k_values(components, temperature_k, pressure_bar)
  let f0 = rachford_rice(ks, z, 0.0)
  let f1 = rachford_rice(ks, z, 1.0)
  if f0 <= 0.0 {
    FlashResult::new(
      temperature_k~,
      pressure_bar~,
      vapor_fraction=0.0,
      liquid=z,
      vapor=z,
      iterations=0,
    )
  } else if f1 >= 0.0 {
    FlashResult::new(
      temperature_k~,
      pressure_bar~,
      vapor_fraction=1.0,
      liquid=z,
      vapor=z,
      iterations=0,
    )
  } else {
    let (vapor_fraction, iterations) = bisect(
      low=0.0,
      high=1.0,
      tolerance=1.0E-10,
      max_iter=150,
      fn(v : Double) -> Double { rachford_rice(ks, z, v) },
    )
    let x = normalize(
      [
        for i in 0.. z[i] / (1.0 + vapor_fraction * (ks[i] - 1.0))
      ],
    )
    let y = normalize([ for i in 0.. ks[i] * x[i] ])
    FlashResult::new(
      temperature_k~,
      pressure_bar~,
      vapor_fraction~,
      liquid=x,
      vapor=y,
      iterations~,
    )
  }
}

///|
fn rachford_rice(
  ks : Array[Double],
  z : Array[Double],
  vapor_fraction : Double,
) -> Double {
  for i = 0, acc = 0.0; i < ks.length(); {
    continue i + 1,
      acc + z[i] * (ks[i] - 1.0) / (1.0 + vapor_fraction * (ks[i] - 1.0))
  } nobreak {
    acc
  }
}

///|
pub fn benchmark_pressure_error(
  point : BenchmarkPoint,
  predicted_pressure_bar : Double,
) -> Double {
  abs_double(predicted_pressure_bar - point.pressure_bar)
}

///|
pub fn benchmark_is_within_pressure_tolerance(
  point : BenchmarkPoint,
  predicted_pressure_bar : Double,
) -> Bool {
  benchmark_pressure_error(point, predicted_pressure_bar) <= point.tolerance
}