///|
/// Lightweight activity-coefficient relations for screening non-ideal data.
pub(all) enum ActivityModel {
  ConstantK(k~ : Double)
  Margules(two_sided~ : Double, three_sided~ : Double)
} derive(Debug, Eq, ToJson)

///|
pub fn ActivityModel::equilibrium(
  self : ActivityModel,
  x : Double,
) -> Double raise @domain.ExtractionError {
  if x < 0.0 || x > 1.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="composition must be in [0,1]",
    )
  }
  match self {
    ConstantK(k~) => {
      @domain.require_non_negative("distribution coefficient", k)
      k * x
    }
    Margules(two_sided~, three_sided~) => {
      @domain.require_finite("Margules two-sided parameter", two_sided)
      @domain.require_finite("Margules three-sided parameter", three_sided)
      let complement = 1.0 - x
      let activity = exp_approx(
        two_sided * complement * complement +
        three_sided * complement * complement * complement,
      )
      x * activity
    }
  }
}

///|
fn exp_approx(value : Double) -> Double {
  let mut term = 1.0
  let mut sum = 1.0
  for index in 1..<=12 {
    term *= value / index.to_double()
    sum += term
  }
  sum
}

///|
pub(all) struct EquilibriumCurve {
  model : ActivityModel
  points : Array[Double]
  values : Array[Double]
  monotone : Bool
  maximum : Double
} derive(Debug, Eq, ToJson)

///|
pub fn sample_curve(
  model : ActivityModel,
  count : Int,
) -> EquilibriumCurve raise @domain.ExtractionError {
  if count < 2 {
    raise @domain.ExtractionError::InsufficientData(required=2, actual=count)
  }
  let points = [
    for index in 0.. index.to_double() / (count - 1).to_double()
  ]
  let values = [ for x in points => model.equilibrium(x) ]
  let mut monotone = true
  for index in 1.. Double raise @domain.ExtractionError {
  if x == 0.0 {
    1.0
  } else {
    model.equilibrium(x) / x
  }
}

///|
pub fn validate_curve(
  curve : EquilibriumCurve,
  tolerance~ : Double,
) -> Array[String] raise @domain.ExtractionError {
  @domain.require_positive("curve tolerance", tolerance)
  let errors = []
  for index, value in curve.values {
    if value < -tolerance {
      errors.push("negative equilibrium value at index {index}")
    }
    if index > 0 && curve.points[index] <= curve.points[index - 1] {
      errors.push("curve x coordinates are not increasing")
    }
  }
  errors
}