///|
pub fn bubble_pressure_bar(
  components : Array[Component],
  liquid : Array[Double],
  temperature_k : Double,
  model? : ActivityModel = ActivityModel::Ideal,
) -> EquilibriumPoint raise VleError {
  assert_same_length(components.length(), liquid.length())
  let x = normalize(liquid)
  let psat = saturation_pressures_bar(components, temperature_k)
  let gamma = activity_coefficients(model, x, temperature_k)
  let p = for i = 0, acc = 0.0; i < x.length(); {
    continue i + 1, acc + x[i] * gamma[i] * psat[i]
  } nobreak {
    acc
  }
  let y_raw = [
    for i = 0; i < x.length(); i = i + 1 => x[i] * gamma[i] * psat[i] / p
  ]
  EquilibriumPoint::new(
    temperature_k~,
    pressure_bar=p,
    liquid=x,
    vapor=normalize(y_raw),
    iterations=1,
  )
}

///|
pub fn dew_pressure_bar(
  components : Array[Component],
  vapor : Array[Double],
  temperature_k : Double,
  model? : ActivityModel = ActivityModel::Ideal,
) -> EquilibriumPoint raise VleError {
  assert_same_length(components.length(), vapor.length())
  let y = normalize(vapor)
  let psat = saturation_pressures_bar(components, temperature_k)
  let x0 = normalize([ for i = 0; i < y.length(); i = i + 1 => y[i] / psat[i] ])
  let gamma = activity_coefficients(model, x0, temperature_k)
  let denom = for i = 0, acc = 0.0; i < y.length(); {
    continue i + 1, acc + y[i] / (gamma[i] * psat[i])
  } nobreak {
    acc
  }
  let p = 1.0 / denom
  let x_raw = [
    for i = 0; i < y.length(); i = i + 1 => y[i] * p / (gamma[i] * psat[i])
  ]
  EquilibriumPoint::new(
    temperature_k~,
    pressure_bar=p,
    liquid=normalize(x_raw),
    vapor=y,
    iterations=1,
  )
}

///|
pub fn bubble_temperature_k(
  components : Array[Component],
  liquid : Array[Double],
  pressure_bar : Double,
  low_k? : Double = 250.0,
  high_k? : Double = 450.0,
  model? : ActivityModel = ActivityModel::Ideal,
) -> EquilibriumPoint raise VleError {
  assert_pressure(pressure_bar)
  let f = fn(t : Double) -> Double raise VleError {
    bubble_pressure_bar(components, liquid, t, model~).pressure_bar -
    pressure_bar
  }
  let (temperature, iterations) = bisect(
    low=low_k,
    high=high_k,
    tolerance=1.0E-7,
    max_iter=100,
    f,
  )
  let point = bubble_pressure_bar(components, liquid, temperature, model~)
  EquilibriumPoint::new(
    temperature_k=temperature,
    pressure_bar~,
    liquid=point.liquid,
    vapor=point.vapor,
    iterations~,
  )
}

///|
pub fn dew_temperature_k(
  components : Array[Component],
  vapor : Array[Double],
  pressure_bar : Double,
  low_k? : Double = 250.0,
  high_k? : Double = 450.0,
  model? : ActivityModel = ActivityModel::Ideal,
) -> EquilibriumPoint raise VleError {
  assert_pressure(pressure_bar)
  let f = fn(t : Double) -> Double raise VleError {
    dew_pressure_bar(components, vapor, t, model~).pressure_bar - pressure_bar
  }
  let (temperature, iterations) = bisect(
    low=low_k,
    high=high_k,
    tolerance=1.0E-7,
    max_iter=100,
    f,
  )
  let point = dew_pressure_bar(components, vapor, temperature, model~)
  EquilibriumPoint::new(
    temperature_k=temperature,
    pressure_bar~,
    liquid=point.liquid,
    vapor=point.vapor,
    iterations~,
  )
}

///|
pub fn relative_volatility(
  component_a : Component,
  component_b : Component,
  temperature_k : Double,
) -> Double raise VleError {
  component_a.saturation_pressure_bar(temperature_k) /
  component_b.saturation_pressure_bar(temperature_k)
}