///|
pub fn binary_composition_grid(
  points : Int,
) -> Array[Array[Double]] raise VleError {
  if points < 2 {
    raise VleError::InvalidParameter("binary composition grid needs two points")
  }
  [
    for i in 0.. {
      let first = 1.0 - i.to_double() / (points - 1).to_double()
      [first, 1.0 - first]
    }
  ]
}

///|
pub fn bubble_pressure_multicomponent(
  components : Array[Component],
  liquid : Array[Double],
  temperature_k : Double,
  model : MulticomponentActivityModel,
) -> EquilibriumPoint raise VleError {
  assert_same_length(components.length(), liquid.length())
  let x = normalize(liquid)
  let psat = saturation_pressures_bar(components, temperature_k)
  let gamma = multicomponent_activity_coefficients(model, x, temperature_k)
  let p = for i = 0, value = 0.0; i < x.length(); {
    continue i + 1, value + x[i] * gamma[i] * psat[i]
  } nobreak {
    value
  }
  let y = [ for i in 0.. x[i] * gamma[i] * psat[i] / p ]
  EquilibriumPoint::new(
    temperature_k~,
    pressure_bar=p,
    liquid=x,
    vapor=normalize(y),
    iterations=1,
  )
}

///|
pub fn dew_pressure_multicomponent(
  components : Array[Component],
  vapor : Array[Double],
  temperature_k : Double,
  model : MulticomponentActivityModel,
) -> EquilibriumPoint raise VleError {
  assert_same_length(components.length(), vapor.length())
  let y = normalize(vapor)
  let psat = saturation_pressures_bar(components, temperature_k)
  let x_guess = normalize([ for i in 0.. y[i] / psat[i] ])
  let gamma = multicomponent_activity_coefficients(
    model, x_guess, temperature_k,
  )
  let denominator = for i = 0, value = 0.0; i < y.length(); {
    continue i + 1, value + y[i] / (gamma[i] * psat[i])
  } nobreak {
    value
  }
  let pressure = 1.0 / denominator
  let x = normalize(
    [
      for i in 0.. y[i] * pressure / (gamma[i] * psat[i])
    ],
  )
  EquilibriumPoint::new(
    temperature_k~,
    pressure_bar=pressure,
    liquid=x,
    vapor=y,
    iterations=1,
  )
}

///|
pub fn bubble_temperature_multicomponent(
  components : Array[Component],
  liquid : Array[Double],
  pressure_bar : Double,
  low_k? : Double = 250.0,
  high_k? : Double = 450.0,
  model : MulticomponentActivityModel,
) -> EquilibriumPoint raise VleError {
  assert_pressure(pressure_bar)
  let options = SolverOptions::new(tolerance=1.0E-7, max_iterations=120)
  let f = fn(t : Double) -> Double raise VleError {
    bubble_pressure_multicomponent(components, liquid, t, model).pressure_bar -
    pressure_bar
  }
  let report = solve_bisection(low=low_k, high=high_k, options~, f)
  let result = bubble_pressure_multicomponent(
    components,
    liquid,
    report.root,
    model,
  )
  EquilibriumPoint::new(
    temperature_k=report.root,
    pressure_bar~,
    liquid=result.liquid,
    vapor=result.vapor,
    iterations=report.iterations,
  )
}

///|
pub fn dew_temperature_multicomponent(
  components : Array[Component],
  vapor : Array[Double],
  pressure_bar : Double,
  low_k? : Double = 250.0,
  high_k? : Double = 450.0,
  model : MulticomponentActivityModel,
) -> EquilibriumPoint raise VleError {
  assert_pressure(pressure_bar)
  let options = SolverOptions::new(tolerance=1.0E-7, max_iterations=120)
  let f = fn(t : Double) -> Double raise VleError {
    dew_pressure_multicomponent(components, vapor, t, model).pressure_bar -
    pressure_bar
  }
  let report = solve_bisection(low=low_k, high=high_k, options~, f)
  let result = dew_pressure_multicomponent(
    components,
    vapor,
    report.root,
    model,
  )
  EquilibriumPoint::new(
    temperature_k=report.root,
    pressure_bar~,
    liquid=result.liquid,
    vapor=result.vapor,
    iterations=report.iterations,
  )
}

///|
pub fn bubble_curve_binary(
  components : Array[Component],
  model : MulticomponentActivityModel,
  pressure_bar~ : Double,
  points~ : Int,
  low_k~ : Double,
  high_k~ : Double,
) -> Array[EquilibriumPoint] raise VleError {
  assert_same_length(2, components.length())
  let grid = binary_composition_grid(points)
  [
    for composition in grid => {
      bubble_temperature_multicomponent(
        components,
        composition,
        pressure_bar,
        low_k~,
        high_k~,
        model,
      )
    }
  ]
}

///|
pub fn dew_curve_binary(
  components : Array[Component],
  model : MulticomponentActivityModel,
  pressure_bar~ : Double,
  points~ : Int,
  low_k~ : Double,
  high_k~ : Double,
) -> Array[EquilibriumPoint] raise VleError {
  assert_same_length(2, components.length())
  let grid = binary_composition_grid(points)
  [
    for composition in grid => {
      let bubble = bubble_temperature_multicomponent(
        components,
        composition,
        pressure_bar,
        low_k~,
        high_k~,
        model,
      )
      dew_temperature_multicomponent(
        components,
        bubble.vapor,
        pressure_bar,
        low_k~,
        high_k~,
        model,
      )
    }
  ]
}

///|
pub fn curve_temperature_span(
  curve : Array[EquilibriumPoint],
) -> Double raise VleError {
  if curve.length() == 0 {
    raise VleError::EmptyMixture
  }
  let minimum = for i = 1, value = curve[0].temperature_k; i < curve.length(); {
    continue i + 1,
      if curve[i].temperature_k < value {
        curve[i].temperature_k
      } else {
        value
      }
  } nobreak {
    value
  }
  let maximum = for i = 1, value = curve[0].temperature_k; i < curve.length(); {
    continue i + 1,
      if curve[i].temperature_k > value {
        curve[i].temperature_k
      } else {
        value
      }
  } nobreak {
    value
  }
  maximum - minimum
}

///|
pub fn curve_pressure_span(
  curve : Array[EquilibriumPoint],
) -> Double raise VleError {
  if curve.length() == 0 {
    raise VleError::EmptyMixture
  }
  let minimum = for i = 1, value = curve[0].pressure_bar; i < curve.length(); {
    continue i + 1,
      if curve[i].pressure_bar < value {
        curve[i].pressure_bar
      } else {
        value
      }
  } nobreak {
    value
  }
  let maximum = for i = 1, value = curve[0].pressure_bar; i < curve.length(); {
    continue i + 1,
      if curve[i].pressure_bar > value {
        curve[i].pressure_bar
      } else {
        value
      }
  } nobreak {
    value
  }
  maximum - minimum
}

///|
pub fn curve_compositions_sum_to_one(
  curve : Array[EquilibriumPoint],
  tolerance : Double,
) -> Bool raise VleError {
  if tolerance <= 0.0 {
    raise VleError::InvalidParameter("curve tolerance must be positive")
  }
  for point in curve {
    if abs_double(sum(point.liquid) - 1.0) > tolerance ||
      abs_double(sum(point.vapor) - 1.0) > tolerance {
      return false
    }
  }
  true
}