///|
/// A binary liquid activity model for non-ideal equilibrium screening.
pub(all) enum NonIdealActivityModel {
  SymmetricMargules(parameter~ : Double)
  TwoSuffixMargules(first_to_second~ : Double, second_to_first~ : Double)
  PolynomialCorrection(first~ : Double, second~ : Double, third~ : Double)
} derive(Debug, Eq, ToJson)

///|
/// Activity coefficients and excess Gibbs-energy diagnostics for a binary pair.
pub(all) struct ActivityPair {
  first_fraction : Double
  second_fraction : Double
  first_coefficient : Double
  second_coefficient : Double
  first_ln_coefficient : Double
  second_ln_coefficient : Double
  excess_gibbs_over_rt : Double
} derive(Debug, Eq, ToJson)

///|
fn clamp_fraction(value : Double) -> Double {
  value.max(1.0e-12).min(1.0 - 1.0e-12)
}

///|
fn approximate_exp(value : Double) -> Double {
  let clipped = value.max(-40.0).min(40.0)
  let mut term = 1.0
  let mut sum = 1.0
  for index in 1..<28 {
    term = term * clipped / index.to_double()
    sum += term
  }
  sum.max(1.0e-12)
}

///|
fn model_ln_gamma(
  model : NonIdealActivityModel,
  first_fraction : Double,
) -> (Double, Double) raise @domain.ExtractionError {
  @domain.require_fraction("first fraction", first_fraction)
  let x = clamp_fraction(first_fraction)
  let y = 1.0 - x
  match model {
    SymmetricMargules(parameter~) => {
      @domain.require_finite("Margules parameter", parameter)
      let first = parameter * y * y
      let second = parameter * x * x
      (first, second)
    }
    TwoSuffixMargules(first_to_second~, second_to_first~) => {
      @domain.require_finite("first Margules parameter", first_to_second)
      @domain.require_finite("second Margules parameter", second_to_first)
      let first = y *
        y *
        (first_to_second + 2.0 * (second_to_first - first_to_second) * x)
      let second = x *
        x *
        (second_to_first + 2.0 * (first_to_second - second_to_first) * y)
      (first, second)
    }
    PolynomialCorrection(first~, second~, third~) => {
      @domain.require_finite("first correction", first)
      @domain.require_finite("second correction", second)
      @domain.require_finite("third correction", third)
      let first_power = x * x
      let second_power = y * y
      let first_value = first * second_power +
        second * x * second_power +
        third * first_power * second_power
      let second_value = first * first_power +
        second * y * first_power +
        third * second_power * first_power
      (first_value, second_value)
    }
  }
}

///|
/// Evaluate activity coefficients for a binary liquid composition.
pub fn activity_pair(
  model : NonIdealActivityModel,
  first_fraction~ : Double,
) -> ActivityPair raise @domain.ExtractionError {
  let second_fraction = 1.0 - first_fraction
  let (first_ln, second_ln) = model_ln_gamma(model, first_fraction)
  {
    first_fraction,
    second_fraction,
    first_coefficient: approximate_exp(first_ln),
    second_coefficient: approximate_exp(second_ln),
    first_ln_coefficient: first_ln,
    second_ln_coefficient: second_ln,
    excess_gibbs_over_rt: first_fraction * first_ln +
    second_fraction * second_ln,
  }
}

///|
pub fn activity_coefficient(
  model : NonIdealActivityModel,
  component_fraction~ : Double,
  component_index~ : Int,
) -> Double raise @domain.ExtractionError {
  let pair = activity_pair(model, first_fraction=component_fraction)
  match component_index {
    1 => pair.first_coefficient
    2 => pair.second_coefficient
    _ =>
      raise @domain.ExtractionError::InvalidArgument(
        message="binary component index must be 1 or 2",
      )
  }
}

///|
pub fn excess_gibbs_energy(
  model : NonIdealActivityModel,
  first_fraction~ : Double,
  temperature~ : @units.Temperature,
) -> Double raise @domain.ExtractionError {
  @domain.require_positive("temperature", temperature.kelvin())
  let gas_constant = 8.314462618
  gas_constant *
  temperature.kelvin() *
  activity_pair(model, first_fraction~).excess_gibbs_over_rt
}

///|
/// Polynomial vapor-pressure correlation in a bounded temperature range.
pub(all) struct VaporPressureCorrelation {
  reference_temperature_kelvin : Double
  coefficient_zero : Double
  coefficient_one : Double
  coefficient_two : Double
  coefficient_three : Double
  minimum_temperature_kelvin : Double
  maximum_temperature_kelvin : Double
  pressure_unit : @units.PressureUnit
} derive(Debug, Eq, ToJson)

///|
pub fn VaporPressureCorrelation::new(
  reference_temperature~ : @units.Temperature,
  coefficient_zero~ : Double,
  coefficient_one~ : Double,
  coefficient_two~ : Double,
  coefficient_three~ : Double,
  minimum_temperature~ : @units.Temperature,
  maximum_temperature~ : @units.Temperature,
  pressure_unit~ : @units.PressureUnit,
) -> VaporPressureCorrelation raise @domain.ExtractionError {
  @domain.require_finite("vapor pressure coefficient zero", coefficient_zero)
  @domain.require_finite("vapor pressure coefficient one", coefficient_one)
  @domain.require_finite("vapor pressure coefficient two", coefficient_two)
  @domain.require_finite("vapor pressure coefficient three", coefficient_three)
  let lower = minimum_temperature.kelvin()
  let upper = maximum_temperature.kelvin()
  if upper <= lower {
    raise @domain.ExtractionError::InvalidArgument(
      message="vapor pressure temperature range must be increasing",
    )
  }
  if reference_temperature.kelvin() < lower ||
    reference_temperature.kelvin() > upper {
    raise @domain.ExtractionError::OutOfRange(
      name="reference temperature",
      value=reference_temperature.kelvin(),
      lower~,
      upper~,
    )
  }
  {
    reference_temperature_kelvin: reference_temperature.kelvin(),
    coefficient_zero,
    coefficient_one,
    coefficient_two,
    coefficient_three,
    minimum_temperature_kelvin: lower,
    maximum_temperature_kelvin: upper,
    pressure_unit,
  }
}

///|
fn polynomial_pressure(
  correlation : VaporPressureCorrelation,
  temperature_kelvin : Double,
) -> Double {
  let delta = temperature_kelvin - correlation.reference_temperature_kelvin
  correlation.coefficient_zero +
  correlation.coefficient_one * delta +
  correlation.coefficient_two * delta * delta +
  correlation.coefficient_three * delta * delta * delta
}

///|
pub fn VaporPressureCorrelation::pressure(
  self : VaporPressureCorrelation,
  temperature : @units.Temperature,
) -> @units.Pressure raise @domain.ExtractionError {
  let kelvin = temperature.kelvin()
  if kelvin < self.minimum_temperature_kelvin ||
    kelvin > self.maximum_temperature_kelvin {
    raise @domain.ExtractionError::OutOfRange(
      name="vapor pressure temperature",
      value=kelvin,
      lower=self.minimum_temperature_kelvin,
      upper=self.maximum_temperature_kelvin,
    )
  }
  let pressure = polynomial_pressure(self, kelvin)
  @units.Pressure::new(value=pressure, unit=self.pressure_unit)
}

///|
pub fn VaporPressureCorrelation::slope(
  self : VaporPressureCorrelation,
  temperature : @units.Temperature,
) -> Double raise @domain.ExtractionError {
  let kelvin = temperature.kelvin()
  if kelvin < self.minimum_temperature_kelvin ||
    kelvin > self.maximum_temperature_kelvin {
    raise @domain.ExtractionError::OutOfRange(
      name="vapor pressure temperature",
      value=kelvin,
      lower=self.minimum_temperature_kelvin,
      upper=self.maximum_temperature_kelvin,
    )
  }
  let delta = kelvin - self.reference_temperature_kelvin
  self.coefficient_one +
  2.0 * self.coefficient_two * delta +
  3.0 * self.coefficient_three * delta * delta
}

///|
pub fn VaporPressureCorrelation::relative_pressure_change(
  self : VaporPressureCorrelation,
  first : @units.Temperature,
  second : @units.Temperature,
) -> Double raise @domain.ExtractionError {
  let first_pressure = self.pressure(first).pascal()
  let second_pressure = self.pressure(second).pascal()
  if first_pressure <= 0.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="reference vapor pressure must be positive",
    )
  }
  (second_pressure - first_pressure) / first_pressure
}

///|
/// A binary equilibrium estimate from activity and vapor-pressure data.
pub(all) struct BinaryEquilibriumPoint {
  liquid_fraction : Double
  vapor_fraction : Double
  first_k_value : Double
  second_k_value : Double
  bubble_pressure : @units.Pressure
  excess_gibbs_energy : Double
} derive(Debug, Eq, ToJson)

///|
pub fn binary_equilibrium_point(
  model : NonIdealActivityModel,
  first_fraction~ : Double,
  temperature~ : @units.Temperature,
  pressure~ : @units.Pressure,
  first_vapor_pressure~ : VaporPressureCorrelation,
  second_vapor_pressure~ : VaporPressureCorrelation,
) -> BinaryEquilibriumPoint raise @domain.ExtractionError {
  @domain.require_positive("system pressure", pressure.pascal())
  let activities = activity_pair(model, first_fraction~)
  let first_psat = first_vapor_pressure.pressure(temperature).pascal()
  let second_psat = second_vapor_pressure.pressure(temperature).pascal()
  let first_k = activities.first_coefficient * first_psat / pressure.pascal()
  let second_k = activities.second_coefficient * second_psat / pressure.pascal()
  let bubble_pa = first_fraction * activities.first_coefficient * first_psat +
    (1.0 - first_fraction) * activities.second_coefficient * second_psat
  let bubble_pressure = @units.Pressure::new(
    value=bubble_pa,
    unit=@units.PressureUnit::Pascal,
  )
  let denominator = first_fraction * first_k + (1.0 - first_fraction) * second_k
  let vapor_fraction = if denominator <= 0.0 {
    0.0
  } else {
    first_fraction * first_k / denominator
  }
  {
    liquid_fraction: first_fraction,
    vapor_fraction: vapor_fraction.max(0.0).min(1.0),
    first_k_value: first_k,
    second_k_value: second_k,
    bubble_pressure,
    excess_gibbs_energy: excess_gibbs_energy(
      model,
      first_fraction~,
      temperature~,
    ),
  }
}

///|
/// Flash calculation result for a multicomponent feed.
pub(all) struct FlashResult {
  vapor_fraction : Double
  liquid_composition : Array[Double]
  vapor_composition : Array[Double]
  feed_composition : Array[Double]
  rachford_rice_residual : Double
  iterations : Int
  converged : Bool
} derive(Debug, Eq, ToJson)

///|
fn rachford_rice(
  feed : ArrayView[Double],
  k_values : ArrayView[Double],
  vapor_fraction : Double,
) -> Double {
  let mut result = 0.0
  for index, feed_value in feed {
    let denominator = 1.0 + vapor_fraction * (k_values[index] - 1.0)
    result += feed_value * (k_values[index] - 1.0) / denominator
  }
  result
}

///|
fn normalize_composition(
  values : ArrayView[Double],
) -> Array[Double] raise @domain.ExtractionError {
  if values.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  let mut total = 0.0
  for value in values {
    @domain.require_non_negative("composition", value)
    total += value
  }
  if total <= 0.0 {
    raise @domain.ExtractionError::InvalidComposition(total~)
  }
  [
    for value in values => value / total
  ]
}

///|
/// Solve the isothermal Rachford–Rice flash equation for supplied K-values.
pub fn isothermal_flash(
  feed_composition~ : ArrayView[Double],
  k_values~ : ArrayView[Double],
  tolerance? : Double = 1.0e-10,
  maximum_iterations? : Int = 200,
) -> FlashResult raise @domain.ExtractionError {
  if feed_composition.length() != k_values.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="feed and K-value dimensions differ",
    )
  }
  @domain.require_positive("flash tolerance", tolerance)
  if maximum_iterations < 1 {
    raise @domain.ExtractionError::InvalidArgument(
      message="flash iteration limit must be positive",
    )
  }
  let feed = normalize_composition(feed_composition)
  for value in k_values {
    @domain.require_positive("K-value", value)
  }
  let mut lower = 0.0
  let mut upper = 1.0
  let lower_value = rachford_rice(feed, k_values, lower)
  let upper_value = rachford_rice(feed, k_values, upper)
  if lower_value < 0.0 {
    let liquid = feed.copy()
    let vapor = [ for value in feed => value ]
    {
      vapor_fraction: 0.0,
      liquid_composition: liquid,
      vapor_composition: vapor,
      feed_composition: feed,
      rachford_rice_residual: lower_value,
      iterations: 0,
      converged: true,
    }
  } else if upper_value > 0.0 {
    let liquid = feed.copy()
    let vapor = feed.copy()
    {
      vapor_fraction: 1.0,
      liquid_composition: liquid,
      vapor_composition: vapor,
      feed_composition: feed,
      rachford_rice_residual: upper_value,
      iterations: 0,
      converged: true,
    }
  } else {
    let mut fraction = 0.5
    let mut residual = rachford_rice(feed, k_values, fraction)
    let mut iterations = 0
    for index in 0.. 0.0 {
        lower = fraction
      } else {
        upper = fraction
      }
    }
    let liquid = [
      for index, feed_value in feed => {
        feed_value / (1.0 + fraction * (k_values[index] - 1.0))
      }
    ]
    let vapor = [
      for index, liquid_value in liquid => k_values[index] * liquid_value
    ]
    let liquid_total = liquid.fold(init=0.0, (sum, value) => sum + value)
    let vapor_total = vapor.fold(init=0.0, (sum, value) => sum + value)
    {
      vapor_fraction: fraction,
      liquid_composition: [
        for value in liquid => value / liquid_total
      ],
      vapor_composition: [
        for value in vapor => value / vapor_total
      ],
      feed_composition: feed,
      rachford_rice_residual: residual,
      iterations,
      converged: residual.abs() <= tolerance,
    }
  }
}

///|
/// Estimate the K-value of a component from activity and vapor pressure.
pub fn k_value(
  activity_coefficient~ : Double,
  vapor_pressure~ : @units.Pressure,
  system_pressure~ : @units.Pressure,
) -> Double raise @domain.ExtractionError {
  @domain.require_positive("activity coefficient", activity_coefficient)
  @domain.require_positive("system pressure", system_pressure.pascal())
  activity_coefficient * vapor_pressure.pascal() / system_pressure.pascal()
}

///|
/// Convert K-values into relative volatility ratios against a reference component.
pub fn relative_volatility(
  k_values : ArrayView[Double],
  reference_index : Int,
) -> Array[Double] raise @domain.ExtractionError {
  if k_values.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  if reference_index < 0 || reference_index >= k_values.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="reference component index is outside K-value array",
    )
  }
  let reference = k_values[reference_index]
  @domain.require_positive("reference K-value", reference)
  [
    for value in k_values => {
      @domain.require_positive("K-value", value)
      value / reference
    }
  ]
}

///|
/// A pressure correction for a simple vapor-phase fugacity factor.
pub fn pressure_correction(
  pressure : @units.Pressure,
  second_virial_coefficient~ : Double,
  temperature : @units.Temperature,
) -> Double raise @domain.ExtractionError {
  @domain.require_finite("second virial coefficient", second_virial_coefficient)
  let gas_constant = 8.314462618
  let z = 1.0 +
    pressure.pascal() *
    second_virial_coefficient /
    (gas_constant * temperature.kelvin())
  z.max(0.05)
}

///|
/// Apply a pressure correction to a K-value while preserving its sign.
pub fn corrected_k_value(
  k_value~ : Double,
  correction~ : Double,
) -> Double raise @domain.ExtractionError {
  @domain.require_positive("K-value", k_value)
  @domain.require_positive("fugacity correction", correction)
  k_value / correction
}