///|
/// Screening calculations for binary equilibrium envelopes.
pub(all) struct EnvelopePoint {
  temperature : Double
  pressure : Double
  liquid_fraction : Double
  vapor_fraction : Double
  liquid_composition : Double
  vapor_composition : Double
  k_value : Double
  phase : String
} derive(Debug, Eq, ToJson)

///|
pub(all) struct EnvelopeCurve {
  component : String
  points : Array[EnvelopePoint]
  bubble_pressure : Double
  dew_pressure : Double
  minimum_temperature : Double
  maximum_temperature : Double
} derive(Debug, Eq, ToJson)

///|
pub(all) struct TieLine {
  liquid_composition : Double
  vapor_composition : Double
  temperature : Double
  lever_liquid : Double
  lever_vapor : Double
  separation : Double
} derive(Debug, Eq, ToJson)

///|
pub fn equilibrium_k_value(
  relative_volatility : Double,
  composition : Double,
) -> Double raise @domain.ExtractionError {
  if relative_volatility <= 0.0 {
    raise @domain.ExtractionError::InvalidCoefficient(
      name="relative volatility",
      value=relative_volatility,
    )
  }
  if composition < 0.0 || composition > 1.0 {
    raise @domain.ExtractionError::InvalidFraction(
      name="composition",
      value=composition,
    )
  }
  relative_volatility / (1.0 + (relative_volatility - 1.0) * composition)
}

///|
pub fn vapor_composition_from_k(
  liquid_composition : Double,
  k_value : Double,
) -> Double raise @domain.ExtractionError {
  if liquid_composition < 0.0 || liquid_composition > 1.0 || k_value <= 0.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="equilibrium composition is invalid",
    )
  }
  let numerator = k_value * liquid_composition
  let denominator = 1.0 + (k_value - 1.0) * liquid_composition
  numerator / denominator
}

///|
pub fn liquid_composition_from_k(
  vapor_composition : Double,
  k_value : Double,
) -> Double raise @domain.ExtractionError {
  if vapor_composition < 0.0 || vapor_composition > 1.0 || k_value <= 0.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="equilibrium composition is invalid",
    )
  }
  let denominator = k_value - (k_value - 1.0) * vapor_composition
  vapor_composition / denominator
}

///|
pub fn bubble_pressure(
  liquid_composition : Double,
  saturation_pressure_one : Double,
  saturation_pressure_two : Double,
) -> Double raise @domain.ExtractionError {
  if liquid_composition < 0.0 ||
    liquid_composition > 1.0 ||
    saturation_pressure_one <= 0.0 ||
    saturation_pressure_two <= 0.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="bubble pressure inputs are invalid",
    )
  }
  liquid_composition * saturation_pressure_one +
  (1.0 - liquid_composition) * saturation_pressure_two
}

///|
pub fn dew_pressure(
  vapor_composition : Double,
  saturation_pressure_one : Double,
  saturation_pressure_two : Double,
) -> Double raise @domain.ExtractionError {
  if vapor_composition < 0.0 ||
    vapor_composition > 1.0 ||
    saturation_pressure_one <= 0.0 ||
    saturation_pressure_two <= 0.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="dew pressure inputs are invalid",
    )
  }
  1.0 /
  (
    vapor_composition / saturation_pressure_one +
    (1.0 - vapor_composition) / saturation_pressure_two
  )
}

///|
pub fn bubble_temperature_estimate(
  pressure : Double,
  normal_boiling_one : Double,
  normal_boiling_two : Double,
  liquid_composition : Double,
) -> Double raise @domain.ExtractionError {
  if pressure <= 0.0 ||
    normal_boiling_one <= 0.0 ||
    normal_boiling_two <= 0.0 ||
    liquid_composition < 0.0 ||
    liquid_composition > 1.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="bubble temperature inputs are invalid",
    )
  }
  let weighted = liquid_composition * normal_boiling_one +
    (1.0 - liquid_composition) * normal_boiling_two
  weighted * (1.0 - 0.08 * (pressure - 1.0) / pressure)
}

///|
pub fn dew_temperature_estimate(
  pressure : Double,
  normal_boiling_one : Double,
  normal_boiling_two : Double,
  vapor_composition : Double,
) -> Double raise @domain.ExtractionError {
  if pressure <= 0.0 ||
    normal_boiling_one <= 0.0 ||
    normal_boiling_two <= 0.0 ||
    vapor_composition < 0.0 ||
    vapor_composition > 1.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="dew temperature inputs are invalid",
    )
  }
  let weighted = vapor_composition * normal_boiling_one +
    (1.0 - vapor_composition) * normal_boiling_two
  weighted * (1.0 - 0.08 * (pressure - 1.0) / pressure)
}

///|
pub fn make_envelope_point(
  temperature : Double,
  pressure : Double,
  liquid_composition : Double,
  relative_volatility : Double,
) -> EnvelopePoint raise @domain.ExtractionError {
  let k = equilibrium_k_value(relative_volatility, liquid_composition)
  let vapor = vapor_composition_from_k(liquid_composition, k)
  {
    temperature,
    pressure,
    liquid_fraction: 1.0 - vapor,
    vapor_fraction: vapor,
    liquid_composition,
    vapor_composition: vapor,
    k_value: k,
    phase: if vapor > 0.999 {
      "vapor"
    } else if vapor < 0.001 {
      "liquid"
    } else {
      "two-phase"
    },
  }
}

///|
pub fn build_envelope(
  component~ : String,
  temperatures : ArrayView[Double],
  pressures : ArrayView[Double],
  relative_volatility : Double,
) -> EnvelopeCurve raise @domain.ExtractionError {
  if temperatures.length() == 0 || pressures.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  let points = []
  for temperature in temperatures {
    for pressure in pressures {
      for index in 0..<11 {
        let composition = index.to_double() / 10.0
        points.push(
          make_envelope_point(
            temperature, pressure, composition, relative_volatility,
          ),
        )
      }
    }
  }
  let mut minimum = temperatures[0]
  let mut maximum = temperatures[0]
  for temperature in temperatures {
    minimum = minimum.min(temperature)
    maximum = maximum.max(temperature)
  }
  {
    component,
    points,
    bubble_pressure: pressures[0],
    dew_pressure: pressures[pressures.length() - 1],
    minimum_temperature: minimum,
    maximum_temperature: maximum,
  }
}

///|
pub fn tie_line(
  liquid_composition : Double,
  vapor_composition : Double,
  overall_composition : Double,
  temperature : Double,
) -> TieLine raise @domain.ExtractionError {
  if liquid_composition < 0.0 ||
    liquid_composition > 1.0 ||
    vapor_composition < 0.0 ||
    vapor_composition > 1.0 ||
    overall_composition < 0.0 ||
    overall_composition > 1.0 ||
    vapor_composition <= liquid_composition {
    raise @domain.ExtractionError::InvalidArgument(
      message="tie line compositions are invalid",
    )
  }
  let lever_vapor = (overall_composition - liquid_composition) /
    (vapor_composition - liquid_composition)
  let lever_liquid = 1.0 - lever_vapor
  {
    liquid_composition,
    vapor_composition,
    temperature,
    lever_liquid,
    lever_vapor,
    separation: vapor_composition - liquid_composition,
  }
}

///|
pub fn tie_line_reconstruct(line : TieLine) -> Double {
  line.lever_liquid * line.liquid_composition +
  line.lever_vapor * line.vapor_composition
}

///|
pub fn envelope_points_at_temperature(
  curve : EnvelopeCurve,
  temperature : Double,
  tolerance : Double,
) -> Array[EnvelopePoint] {
  let selected = []
  for point in curve.points {
    if (point.temperature - temperature).abs() <= tolerance {
      selected.push(point)
    }
  }
  selected
}

///|
pub fn envelope_points_at_pressure(
  curve : EnvelopeCurve,
  pressure : Double,
  tolerance : Double,
) -> Array[EnvelopePoint] {
  let selected = []
  for point in curve.points {
    if (point.pressure - pressure).abs() <= tolerance {
      selected.push(point)
    }
  }
  selected
}

///|
pub fn envelope_minimum_separation(curve : EnvelopeCurve) -> Double {
  if curve.points.length() == 0 {
    return 0.0
  }
  let mut result = 1.0e30
  for point in curve.points {
    let separation = point.vapor_composition - point.liquid_composition
    if separation < result {
      result = separation
    }
  }
  result
}

///|
pub fn envelope_maximum_k(curve : EnvelopeCurve) -> Double {
  let mut result = 0.0
  for point in curve.points {
    result = result.max(point.k_value)
  }
  result
}

///|
pub fn envelope_average_temperature(curve : EnvelopeCurve) -> Double {
  if curve.points.length() == 0 {
    return 0.0
  }
  let mut total = 0.0
  for point in curve.points {
    total += point.temperature
  }
  total / curve.points.length().to_double()
}

///|
pub fn envelope_filter_phase(
  curve : EnvelopeCurve,
  phase : String,
) -> Array[EnvelopePoint] {
  let selected = []
  for point in curve.points {
    if point.phase == phase {
      selected.push(point)
    }
  }
  selected
}

///|
pub fn envelope_csv(curve : EnvelopeCurve) -> String {
  let lines = [
    "temperature,pressure,liquid_fraction,vapor_fraction,liquid_x,vapor_y,k,phase",
  ]
  for point in curve.points {
    lines.push(
      point.temperature.to_string() +
      "," +
      point.pressure.to_string() +
      "," +
      point.liquid_fraction.to_string() +
      "," +
      point.vapor_fraction.to_string() +
      "," +
      point.liquid_composition.to_string() +
      "," +
      point.vapor_composition.to_string() +
      "," +
      point.k_value.to_string() +
      "," +
      point.phase,
    )
  }
  lines.join("\n")
}

///|
pub fn envelope_markdown(curve : EnvelopeCurve) -> String {
  let lines = [
    "# Equilibrium envelope: " + curve.component,
    "",
    "| T | P | x | y | K | Phase |",
    "|---:|---:|---:|---:|---:|---|",
  ]
  for point in curve.points {
    lines.push(
      "| " +
      point.temperature.to_string() +
      " | " +
      point.pressure.to_string() +
      " | " +
      point.liquid_composition.to_string() +
      " | " +
      point.vapor_composition.to_string() +
      " | " +
      point.k_value.to_string() +
      " | " +
      point.phase +
      " |",
    )
  }
  lines.join("\n")
}

///|
pub fn interpolate_envelope_point(
  first : EnvelopePoint,
  second : EnvelopePoint,
  temperature : Double,
) -> EnvelopePoint raise @domain.ExtractionError {
  let span = second.temperature - first.temperature
  if span == 0.0 {
    raise @domain.ExtractionError::DuplicateCoordinate(value=temperature)
  }
  let fraction = (temperature - first.temperature) / span
  {
    temperature,
    pressure: first.pressure + fraction * (second.pressure - first.pressure),
    liquid_fraction: first.liquid_fraction +
    fraction * (second.liquid_fraction - first.liquid_fraction),
    vapor_fraction: first.vapor_fraction +
    fraction * (second.vapor_fraction - first.vapor_fraction),
    liquid_composition: first.liquid_composition +
    fraction * (second.liquid_composition - first.liquid_composition),
    vapor_composition: first.vapor_composition +
    fraction * (second.vapor_composition - first.vapor_composition),
    k_value: first.k_value + fraction * (second.k_value - first.k_value),
    phase: if fraction < 0.5 {
      first.phase
    } else {
      second.phase
    },
  }
}

///|
pub fn envelope_recovery_factor(
  curve : EnvelopeCurve,
  target_separation : Double,
) -> Double {
  if target_separation <= 0.0 {
    return 0.0
  }
  (envelope_minimum_separation(curve) / target_separation).min(1.0).max(0.0)
}

///|
pub fn classify_relative_volatility(relative_volatility : Double) -> String {
  if relative_volatility < 1.05 {
    "difficult"
  } else if relative_volatility < 1.5 {
    "moderate"
  } else if relative_volatility < 3.0 {
    "favorable"
  } else {
    "easy"
  }
}