///|
fn activity_flash_k_values(
  components : Array[Component],
  liquid : Array[Double],
  temperature_k : Double,
  pressure_bar : Double,
  model : MulticomponentActivityModel,
) -> Array[Double] raise VleError {
  let psat = saturation_pressures_bar(components, temperature_k)
  let gamma = multicomponent_activity_coefficients(model, liquid, temperature_k)
  [
    for i in 0.. gamma[i] * psat[i] / pressure_bar
  ]
}

///|
fn make_activity_flash_result(
  temperature_k : Double,
  pressure_bar : Double,
  vapor_fraction : Double,
  composition : Array[Double],
  k_values : Array[Double],
  iterations : Int,
  residual : Double,
) -> ActivityFlashResult raise VleError {
  let z = normalize(composition)
  if vapor_fraction == 0.0 || vapor_fraction == 1.0 {
    return ActivityFlashResult::new(
      temperature_k~,
      pressure_bar~,
      vapor_fraction~,
      liquid=z,
      vapor=z,
      k_values~,
      iterations~,
      converged=true,
      residual~,
    )
  }
  let x = normalize(
    [
      for i in 0.. {
        z[i] / (1.0 + vapor_fraction * (k_values[i] - 1.0))
      }
    ],
  )
  let y = normalize([ for i in 0.. k_values[i] * x[i] ])
  ActivityFlashResult::new(
    temperature_k~,
    pressure_bar~,
    vapor_fraction~,
    liquid=x,
    vapor=y,
    k_values~,
    iterations~,
    converged=true,
    residual~,
  )
}

///|
pub fn flash_isothermal_activity(
  components : Array[Component],
  feed : Array[Double],
  temperature_k : Double,
  pressure_bar : Double,
  model : MulticomponentActivityModel,
  damping? : Double = 0.5,
  tolerance? : Double = 1.0E-8,
  max_iterations? : Int = 100,
) -> ActivityFlashResult raise VleError {
  assert_pressure(pressure_bar)
  assert_same_length(components.length(), feed.length())
  let z = normalize(feed)
  if damping <= 0.0 || damping > 1.0 {
    raise VleError::InvalidParameter("flash damping must be in (0, 1]")
  }
  if tolerance <= 0.0 || max_iterations <= 0 {
    raise VleError::InvalidParameter("flash iteration options must be positive")
  }
  let initial_k = activity_flash_k_values(
    components, z, temperature_k, pressure_bar, model,
  )
  for iteration = 0, current_k = initial_k; iteration < max_iterations; {
    let f0 = rachford_rice(current_k, z, 0.0)
    let f1 = rachford_rice(current_k, z, 1.0)
    if f0 <= 0.0 {
      break make_activity_flash_result(
        temperature_k, pressure_bar, 0.0, z, current_k, iteration, 0.0,
      )
    }
    if f1 >= 0.0 {
      break make_activity_flash_result(
        temperature_k, pressure_bar, 1.0, z, current_k, iteration, 0.0,
      )
    }
    let rr_options = SolverOptions::new(tolerance=1.0E-10, max_iterations=160)
    let rr = solve_bisection(low=0.0, high=1.0, options=rr_options, fn(
      v : Double,
    ) -> Double {
      rachford_rice(current_k, z, v)
    })
    let x = normalize(
      [
        for i in 0.. z[i] / (1.0 + rr.root * (current_k[i] - 1.0))
      ],
    )
    let raw_k = activity_flash_k_values(
      components, x, temperature_k, pressure_bar, model,
    )
    let next_k = [
      for i in 0.. {
        current_k[i] + damping * (raw_k[i] - current_k[i])
      }
    ]
    let residual = weighted_rms(current_k, next_k)
    if residual <= tolerance {
      break make_activity_flash_result(
        temperature_k,
        pressure_bar,
        rr.root,
        z,
        next_k,
        iteration + 1,
        residual,
      )
    }
    continue iteration + 1, next_k
  } nobreak {
    raise VleError::SolverDidNotConverge(iterations=max_iterations)
  }
}

///|
pub fn activity_flash_residual(result : ActivityFlashResult) -> Double {
  result.residual
}

///|
pub fn activity_flash_is_two_phase(result : ActivityFlashResult) -> Bool {
  result.vapor_fraction > 0.0 && result.vapor_fraction < 1.0
}

///|
pub fn activity_flash_material_balance_error(
  result : ActivityFlashResult,
  feed : Array[Double],
) -> Double raise VleError {
  let z = normalize(feed)
  assert_same_length(z.length(), result.liquid.length())
  let reconstructed = [
    for i in 0.. {
      (1.0 - result.vapor_fraction) * result.liquid[i] +
      result.vapor_fraction * result.vapor[i]
    }
  ]
  vector_distance(z, reconstructed)
}

///|
pub fn activity_flash_k_error(
  result : ActivityFlashResult,
  components : Array[Component],
  temperature_k : Double,
  pressure_bar : Double,
  model : MulticomponentActivityModel,
) -> Double raise VleError {
  let expected = activity_flash_k_values(
    components,
    result.liquid,
    temperature_k,
    pressure_bar,
    model,
  )
  weighted_rms(expected, result.k_values)
}

///|
pub fn activity_flash_phase_fractions(
  result : ActivityFlashResult,
) -> Array[Double] {
  [1.0 - result.vapor_fraction, result.vapor_fraction]
}