///|
pub fn generate_equilibrium_profile(
  relative_volatility~ : Double,
  points~ : Int,
) -> Array[EquilibriumProfilePoint] raise VleError {
  if relative_volatility <= 0.0 {
    raise VleError::InvalidParameter("relative volatility must be positive")
  }
  if points < 2 {
    raise VleError::InvalidParameter("profile needs at least two points")
  }
  [
    for i in 0.. {
      let x = i.to_double() / (points - 1).to_double()
      EquilibriumProfilePoint::{
        liquid_fraction: x,
        vapor_fraction: binary_equilibrium_y(x, relative_volatility),
        relative_volatility,
      }
    }
  ]
}

///|
pub fn profile_area(
  profile : Array[EquilibriumProfilePoint],
) -> Double raise VleError {
  if profile.length() < 2 {
    raise VleError::InvalidParameter("profile area needs two points")
  }
  for i = 0, area = 0.0; i + 1 < profile.length(); {
    let width = profile[i + 1].liquid_fraction - profile[i].liquid_fraction
    if width < 0.0 {
      raise VleError::InvalidParameter("profile must be ordered")
    }
    let gap0 = profile[i].vapor_fraction - profile[i].liquid_fraction
    let gap1 = profile[i + 1].vapor_fraction - profile[i + 1].liquid_fraction
    continue i + 1, area + width * (gap0 + gap1) / 2.0
  } nobreak {
    area
  }
}

///|
pub fn profile_is_ordered(profile : Array[EquilibriumProfilePoint]) -> Bool {
  for i = 0; i + 1 < profile.length(); i = i + 1 {
    if profile[i + 1].liquid_fraction < profile[i].liquid_fraction ||
      profile[i + 1].vapor_fraction < profile[i].vapor_fraction {
      return false
    }
  }
  true
}

///|
pub fn profile_maximum_gap(
  profile : Array[EquilibriumProfilePoint],
) -> Double raise VleError {
  if profile.length() == 0 {
    raise VleError::EmptyMixture
  }
  for i = 0, value = 0.0; i < profile.length(); {
    let gap = profile[i].vapor_fraction - profile[i].liquid_fraction
    continue i + 1, if gap > value { gap } else { value }
  } nobreak {
    value
  }
}

///|
pub fn profile_minimum_gap(
  profile : Array[EquilibriumProfilePoint],
) -> Double raise VleError {
  if profile.length() == 0 {
    raise VleError::EmptyMixture
  }
  for i = 0, value = profile[0].vapor_fraction - profile[0].liquid_fraction; i <
     profile.length(); {
    let gap = profile[i].vapor_fraction - profile[i].liquid_fraction
    continue i + 1, if gap < value { gap } else { value }
  } nobreak {
    value
  }
}

///|
pub fn profile_average_gap(
  profile : Array[EquilibriumProfilePoint],
) -> Double raise VleError {
  if profile.length() == 0 {
    raise VleError::EmptyMixture
  }
  for i = 0, value = 0.0; i < profile.length(); {
    continue i + 1,
      value + profile[i].vapor_fraction - profile[i].liquid_fraction
  } nobreak {
    value / profile.length().to_double()
  }
}

///|
pub fn profile_interpolate(
  profile : Array[EquilibriumProfilePoint],
  liquid_fraction : Double,
) -> Double raise VleError {
  if profile.length() < 2 {
    raise VleError::InvalidParameter("profile interpolation needs two points")
  }
  if liquid_fraction < profile[0].liquid_fraction ||
    liquid_fraction > profile[profile.length() - 1].liquid_fraction {
    raise VleError::InvalidParameter("liquid fraction is outside profile")
  }
  for i = 0; i + 1 < profile.length(); i = i + 1 {
    let left = profile[i]
    let right = profile[i + 1]
    if liquid_fraction >= left.liquid_fraction &&
      liquid_fraction <= right.liquid_fraction {
      return linear_interpolate(
        left.liquid_fraction,
        left.vapor_fraction,
        right.liquid_fraction,
        right.vapor_fraction,
        liquid_fraction,
      )
    }
  }
  profile[profile.length() - 1].vapor_fraction
}

///|
pub fn equilibrium_curve_slope(
  relative_volatility : Double,
  liquid_fraction : Double,
) -> Double raise VleError {
  if relative_volatility <= 0.0 ||
    liquid_fraction < 0.0 ||
    liquid_fraction > 1.0 {
    raise VleError::InvalidParameter("equilibrium slope inputs are invalid")
  }
  let denominator = 1.0 + (relative_volatility - 1.0) * liquid_fraction
  relative_volatility / (denominator * denominator)
}

///|
pub fn profile_gap_integral(
  profile : Array[EquilibriumProfilePoint],
  low_fraction : Double,
  high_fraction : Double,
) -> Double raise VleError {
  if high_fraction <= low_fraction {
    raise VleError::InvalidRange(low=low_fraction, high=high_fraction)
  }
  let left = profile_interpolate(profile, low_fraction)
  let right = profile_interpolate(profile, high_fraction)
  (high_fraction - low_fraction) *
  (left - low_fraction + (right - high_fraction)) /
  2.0
}

///|
pub fn profile_separation_index(
  profile : Array[EquilibriumProfilePoint],
) -> Double raise VleError {
  let area = profile_area(profile)
  let width = profile[profile.length() - 1].liquid_fraction -
    profile[0].liquid_fraction
  if width <= 0.0 {
    raise VleError::InvalidRange(
      low=profile[0].liquid_fraction,
      high=profile[profile.length() - 1].liquid_fraction,
    )
  }
  area / width
}

///|
pub fn operating_line_y(
  liquid_fraction : Double,
  reflux_ratio : Double,
  distillate_fraction : Double,
) -> Double raise VleError {
  if reflux_ratio <= 0.0 ||
    distillate_fraction < 0.0 ||
    distillate_fraction > 1.0 {
    raise VleError::InvalidParameter("operating line inputs are invalid")
  }
  let slope = reflux_ratio / (reflux_ratio + 1.0)
  let intercept = distillate_fraction / (reflux_ratio + 1.0)
  slope * liquid_fraction + intercept
}

///|
pub fn operating_line_gap(
  profile_point : EquilibriumProfilePoint,
  reflux_ratio : Double,
  distillate_fraction : Double,
) -> Double raise VleError {
  operating_line_y(
    profile_point.liquid_fraction,
    reflux_ratio,
    distillate_fraction,
  ) -
  profile_point.vapor_fraction
}

///|
pub fn profile_operating_line_crossings(
  profile : Array[EquilibriumProfilePoint],
  reflux_ratio : Double,
  distillate_fraction : Double,
) -> Int raise VleError {
  if profile.length() < 2 {
    raise VleError::InvalidParameter("profile needs two points")
  }
  for i = 0, crossings = 0; i + 1 < profile.length(); {
    let left = operating_line_gap(profile[i], reflux_ratio, distillate_fraction)
    let right = operating_line_gap(
      profile[i + 1],
      reflux_ratio,
      distillate_fraction,
    )
    let crosses = if left == 0.0 || left * right < 0.0 { 1 } else { 0 }
    continue i + 1, crossings + crosses
  } nobreak {
    crossings
  }
}

///|
pub fn profile_points_above_operating_line(
  profile : Array[EquilibriumProfilePoint],
  reflux_ratio : Double,
  distillate_fraction : Double,
) -> Int raise VleError {
  for i = 0, count = 0; i < profile.length(); {
    let gap = operating_line_gap(profile[i], reflux_ratio, distillate_fraction)
    continue i + 1, count + (if gap < 0.0 { 1 } else { 0 })
  } nobreak {
    count
  }
}