///|
fn validate_recovery(value : Double) -> Unit raise VleError {
  if value < 0.0 || value > 1.0 {
    raise VleError::InvalidParameter("recovery must be between zero and one")
  }
}

///|
pub fn split_binary_stream(
  feed : Array[Double],
  light_key_recovery~ : Double,
  heavy_key_recovery~ : Double,
) -> BinarySplitResult raise VleError {
  assert_same_length(2, feed.length())
  validate_recovery(light_key_recovery)
  validate_recovery(heavy_key_recovery)
  let z = normalize(feed)
  let distillate_flows = [z[0] * light_key_recovery, z[1] * heavy_key_recovery]
  let bottoms_flows = [z[0] - distillate_flows[0], z[1] - distillate_flows[1]]
  let distillate_total = sum(distillate_flows)
  let bottoms_total = sum(bottoms_flows)
  if distillate_total <= 0.0 || bottoms_total <= 0.0 {
    raise VleError::InvalidParameter("split must produce two non-empty streams")
  }
  let distillate = normalize(distillate_flows)
  let bottoms = normalize(bottoms_flows)
  let error = for i = 0, value = 0.0; i < 2; {
    continue i + 1,
      value +
      abs_double(
        z[i] - distillate_total * distillate[i] - bottoms_total * bottoms[i],
      )
  } nobreak {
    value
  }
  BinarySplitResult::new(
    feed=z,
    distillate~,
    bottoms~,
    distillate_total~,
    bottoms_total~,
    component_balance_error=error,
  )
}

///|
pub fn fenske_minimum_stages(
  relative_volatility~ : Double,
  light_key_distillate_fraction~ : Double,
  light_key_bottoms_fraction~ : Double,
) -> Double raise VleError {
  if relative_volatility <= 1.0 {
    raise VleError::InvalidParameter("relative volatility must exceed one")
  }
  if light_key_distillate_fraction <= 0.0 ||
    light_key_distillate_fraction >= 1.0 ||
    light_key_bottoms_fraction <= 0.0 ||
    light_key_bottoms_fraction >= 1.0 {
    raise VleError::InvalidParameter("Fenske fractions must be inside (0, 1)")
  }
  let distillate_ratio = light_key_distillate_fraction /
    (1.0 - light_key_distillate_fraction)
  let bottoms_ratio = light_key_bottoms_fraction /
    (1.0 - light_key_bottoms_fraction)
  @math.ln(distillate_ratio / bottoms_ratio) / @math.ln(relative_volatility)
}

///|
pub fn minimum_reflux_ratio_binary(
  feed_light_fraction~ : Double,
  distillate_light_fraction~ : Double,
  relative_volatility~ : Double,
) -> Double raise VleError {
  if feed_light_fraction <= 0.0 ||
    feed_light_fraction >= 1.0 ||
    distillate_light_fraction <= feed_light_fraction ||
    distillate_light_fraction >= 1.0 {
    raise VleError::InvalidParameter(
      "feed and distillate fractions are inconsistent",
    )
  }
  if relative_volatility <= 1.0 {
    raise VleError::InvalidParameter("relative volatility must exceed one")
  }
  let y_intersection = relative_volatility *
    feed_light_fraction /
    (1.0 + (relative_volatility - 1.0) * feed_light_fraction)
  if y_intersection <= feed_light_fraction ||
    distillate_light_fraction <= y_intersection {
    raise VleError::InvalidParameter("minimum reflux pinch is not separable")
  }
  (distillate_light_fraction - y_intersection) /
  (y_intersection - feed_light_fraction)
}

///|
pub fn shortcut_distillation_binary(
  feed : Array[Double],
  relative_volatility~ : Double,
  light_key_distillate_fraction~ : Double,
  light_key_bottoms_fraction~ : Double,
  reflux_multiplier? : Double = 1.5,
) -> ShortcutColumnResult raise VleError {
  if reflux_multiplier <= 1.0 {
    raise VleError::InvalidParameter("reflux multiplier must exceed one")
  }
  let split = split_binary_stream(
    feed,
    light_key_recovery=light_key_distillate_fraction,
    heavy_key_recovery=1.0 - light_key_bottoms_fraction,
  )
  let minimum_stages = fenske_minimum_stages(
    relative_volatility~,
    light_key_distillate_fraction~,
    light_key_bottoms_fraction~,
  )
  let minimum_reflux = minimum_reflux_ratio_binary(
    feed_light_fraction=split.feed[0],
    distillate_light_fraction=split.distillate[0],
    relative_volatility~,
  )
  let operating_reflux = minimum_reflux * reflux_multiplier
  let stage_factor = 1.0 + 1.0 / (operating_reflux - minimum_reflux + 1.0)
  let estimated_stages = @math.ceil(minimum_stages * stage_factor).to_int()
  ShortcutColumnResult::{
    minimum_stages,
    minimum_reflux,
    operating_reflux,
    estimated_stages,
    feed: split.feed,
    distillate: split.distillate,
    bottoms: split.bottoms,
  }
}

///|
pub fn stream_component_flows(
  composition : Array[Double],
  total_moles : Double,
) -> Array[Double] raise VleError {
  if total_moles < 0.0 {
    raise VleError::InvalidParameter("total moles cannot be negative")
  }
  let x = normalize(composition)
  [
    for value in x => value * total_moles
  ]
}

///|
pub fn stream_total_moles(flows : Array[Double]) -> Double raise VleError {
  for value in flows {
    if value < 0.0 {
      raise VleError::InvalidParameter("component flow cannot be negative")
    }
  }
  sum(flows)
}

///|
pub fn stream_composition_from_flows(
  flows : Array[Double],
) -> Array[Double] raise VleError {
  ignore(stream_total_moles(flows))
  normalize(flows)
}

///|
pub fn combine_streams(
  first_flows : Array[Double],
  second_flows : Array[Double],
) -> Array[Double] raise VleError {
  assert_same_length(first_flows.length(), second_flows.length())
  for value in first_flows {
    if value < 0.0 {
      raise VleError::InvalidParameter("first stream contains negative flow")
    }
  }
  for value in second_flows {
    if value < 0.0 {
      raise VleError::InvalidParameter("second stream contains negative flow")
    }
  }
  [
    for i in 0.. first_flows[i] + second_flows[i]
  ]
}

///|
pub fn split_stream_by_fraction(
  flows : Array[Double],
  first_fraction : Double,
) -> (Array[Double], Array[Double]) raise VleError {
  validate_recovery(first_fraction)
  for value in flows {
    if value < 0.0 {
      raise VleError::InvalidParameter("stream flow cannot be negative")
    }
  }
  (
    [
      for value in flows => value * first_fraction
    ],
    [
      for value in flows => value * (1.0 - first_fraction)
    ],
  )
}

///|
pub fn component_recovery(
  feed_flow : Double,
  product_flow : Double,
) -> Double raise VleError {
  if feed_flow <= 0.0 || product_flow < 0.0 {
    raise VleError::InvalidParameter("flows are inconsistent")
  }
  let recovery = product_flow / feed_flow
  validate_recovery(recovery)
  recovery
}

///|
pub fn binary_equilibrium_y(
  liquid_light_fraction : Double,
  relative_volatility : Double,
) -> Double raise VleError {
  if liquid_light_fraction < 0.0 ||
    liquid_light_fraction > 1.0 ||
    relative_volatility <= 0.0 {
    raise VleError::InvalidParameter("binary equilibrium inputs are invalid")
  }
  relative_volatility *
  liquid_light_fraction /
  (1.0 + (relative_volatility - 1.0) * liquid_light_fraction)
}

///|
pub fn binary_equilibrium_x(
  vapor_light_fraction : Double,
  relative_volatility : Double,
) -> Double raise VleError {
  if vapor_light_fraction < 0.0 ||
    vapor_light_fraction > 1.0 ||
    relative_volatility <= 0.0 {
    raise VleError::InvalidParameter("binary equilibrium inputs are invalid")
  }
  vapor_light_fraction /
  (relative_volatility - (relative_volatility - 1.0) * vapor_light_fraction)
}

///|
pub fn q_line_intersection(
  feed_light_fraction : Double,
  q_value : Double,
  equilibrium_relative_volatility : Double,
) -> (Double, Double) raise VleError {
  if feed_light_fraction < 0.0 ||
    feed_light_fraction > 1.0 ||
    equilibrium_relative_volatility <= 0.0 {
    raise VleError::InvalidParameter("q-line inputs are invalid")
  }
  if q_value == 1.0 {
    let y = binary_equilibrium_y(
      feed_light_fraction, equilibrium_relative_volatility,
    )
    return (feed_light_fraction, y)
  }
  let slope = q_value / (q_value - 1.0)
  let intercept = 0.0 - feed_light_fraction / (q_value - 1.0)
  let x = feed_light_fraction
  (x, slope * x + intercept)
}

///|
pub fn material_balance_error(
  feed : Array[Double],
  distillate : Array[Double],
  bottoms : Array[Double],
  distillate_total : Double,
  bottoms_total : Double,
) -> Double raise VleError {
  assert_same_length(feed.length(), distillate.length())
  assert_same_length(feed.length(), bottoms.length())
  let z = normalize(feed)
  let x = normalize(distillate)
  let b = normalize(bottoms)
  for i = 0, value = 0.0; i < z.length(); {
    continue i + 1,
      value + abs_double(z[i] - distillate_total * x[i] - bottoms_total * b[i])
  } nobreak {
    value
  }
}