///|
fn properties_components(
  properties : Array[ThermoProperty],
) -> Array[Component] {
  [
    for property in properties => property.component
  ]
}

///|
pub fn mixture_liquid_enthalpy(
  properties : Array[ThermoProperty],
  composition : Array[Double],
  temperature_k : Double,
) -> Double raise VleError {
  assert_same_length(properties.length(), composition.length())
  let x = normalize(composition)
  for i = 0, value = 0.0; i < properties.length(); {
    continue i + 1,
      value + x[i] * component_liquid_enthalpy(properties[i], temperature_k)
  } nobreak {
    value
  }
}

///|
pub fn mixture_vapor_enthalpy(
  properties : Array[ThermoProperty],
  composition : Array[Double],
  temperature_k : Double,
) -> Double raise VleError {
  assert_same_length(properties.length(), composition.length())
  let y = normalize(composition)
  for i = 0, value = 0.0; i < properties.length(); {
    continue i + 1,
      value + y[i] * component_vapor_enthalpy(properties[i], temperature_k)
  } nobreak {
    value
  }
}

///|
pub fn mixture_phase_enthalpy(
  properties : Array[ThermoProperty],
  composition : Array[Double],
  temperature_k : Double,
  phase : PhaseKind,
) -> Double raise VleError {
  match phase {
    Liquid => mixture_liquid_enthalpy(properties, composition, temperature_k)
    Vapor => mixture_vapor_enthalpy(properties, composition, temperature_k)
    Solid => mixture_liquid_enthalpy(properties, composition, temperature_k)
  }
}

///|
fn overall_flash_enthalpy(
  properties : Array[ThermoProperty],
  flash : FlashResult,
) -> Double raise VleError {
  let liquid_h = mixture_liquid_enthalpy(
    properties,
    flash.liquid,
    flash.temperature_k,
  )
  let vapor_h = mixture_vapor_enthalpy(
    properties,
    flash.vapor,
    flash.temperature_k,
  )
  (1.0 - flash.vapor_fraction) * liquid_h + flash.vapor_fraction * vapor_h
}

///|
pub fn flash_enthalpy_at_temperature(
  properties : Array[ThermoProperty],
  feed : Array[Double],
  temperature_k : Double,
  pressure_bar : Double,
) -> Double raise VleError {
  assert_same_length(properties.length(), feed.length())
  let flash = flash_isothermal_ideal(
    properties_components(properties),
    feed,
    temperature_k,
    pressure_bar,
  )
  overall_flash_enthalpy(properties, flash)
}

///|
pub fn flash_isenthalpic_ideal(
  properties : Array[ThermoProperty],
  feed : Array[Double],
  feed_temperature_k~ : Double,
  pressure_bar~ : Double,
  target_enthalpy_j_per_mol~ : Double,
  low_k? : Double = 250.0,
  high_k? : Double = 450.0,
) -> EnergyFlashResult raise VleError {
  assert_same_length(properties.length(), feed.length())
  assert_temperature(feed_temperature_k)
  assert_pressure(pressure_bar)
  let options = SolverOptions::new(tolerance=1.0E-7, max_iterations=120)
  let objective = fn(temperature_k : Double) -> Double raise VleError {
    flash_enthalpy_at_temperature(properties, feed, temperature_k, pressure_bar) -
    target_enthalpy_j_per_mol
  }
  let report = solve_bisection(low=low_k, high=high_k, options~, objective)
  let flash = flash_isothermal_ideal(
    properties_components(properties),
    feed,
    report.root,
    pressure_bar,
  )
  EnergyFlashResult::new(
    temperature_k=report.root,
    pressure_bar~,
    vapor_fraction=flash.vapor_fraction,
    liquid=flash.liquid,
    vapor=flash.vapor,
    enthalpy_j_per_mol=overall_flash_enthalpy(properties, flash),
    iterations=report.iterations + flash.iterations,
  )
}

///|
pub fn sensible_energy_between(
  property : ThermoProperty,
  low_temperature_k : Double,
  high_temperature_k : Double,
  phase : PhaseKind,
) -> Double raise VleError {
  match phase {
    Liquid =>
      property.liquid_heat_capacity.integral(
        low_k=low_temperature_k,
        high_k=high_temperature_k,
      )
    Vapor =>
      property.vapor_heat_capacity.integral(
        low_k=low_temperature_k,
        high_k=high_temperature_k,
      )
    Solid =>
      property.liquid_heat_capacity.integral(
        low_k=low_temperature_k,
        high_k=high_temperature_k,
      )
  }
}

///|
pub fn latent_energy_fraction(
  property : ThermoProperty,
  temperature_k : Double,
  critical_temperature_k : Double,
) -> Double raise VleError {
  latent_heat_at(property, temperature_k, critical_temperature_k) /
  property.latent_heat_j_per_mol
}

///|
pub fn phase_energy_gap(
  property : ThermoProperty,
  temperature_k : Double,
) -> Double raise VleError {
  component_vapor_enthalpy(property, temperature_k) -
  component_liquid_enthalpy(property, temperature_k)
}

///|
pub fn mixture_energy_gap(
  properties : Array[ThermoProperty],
  composition : Array[Double],
  temperature_k : Double,
) -> Double raise VleError {
  mixture_vapor_enthalpy(properties, composition, temperature_k) -
  mixture_liquid_enthalpy(properties, composition, temperature_k)
}

///|
pub fn weighted_phase_enthalpy(
  liquid_enthalpy : Double,
  vapor_enthalpy : Double,
  vapor_fraction : Double,
) -> Double raise VleError {
  if vapor_fraction < 0.0 || vapor_fraction > 1.0 {
    raise VleError::InvalidParameter(
      "vapor fraction must be between zero and one",
    )
  }
  (1.0 - vapor_fraction) * liquid_enthalpy + vapor_fraction * vapor_enthalpy
}

///|
pub fn estimate_vapor_fraction_from_enthalpy(
  liquid_enthalpy : Double,
  vapor_enthalpy : Double,
  target_enthalpy : Double,
) -> Double raise VleError {
  if vapor_enthalpy <= liquid_enthalpy {
    raise VleError::InvalidParameter(
      "vapor enthalpy must exceed liquid enthalpy",
    )
  }
  let fraction = (target_enthalpy - liquid_enthalpy) /
    (vapor_enthalpy - liquid_enthalpy)
  if fraction < 0.0 || fraction > 1.0 {
    raise VleError::InvalidParameter(
      "target enthalpy is outside the phase envelope",
    )
  }
  fraction
}