///|
fn polynomial_value(
  coefficients : Array[Double],
  temperature_k : Double,
) -> Double {
  for i = coefficients.length() - 1, value = 0.0; i >= 0; {
    continue i - 1, value * temperature_k + coefficients[i]
  } nobreak {
    value
  }
}

///|
pub fn HeatCapacityPolynomial::is_in_range(
  self : HeatCapacityPolynomial,
  temperature_k : Double,
) -> Bool {
  temperature_k >= self.t_min_k && temperature_k <= self.t_max_k
}

///|
pub fn HeatCapacityPolynomial::evaluate(
  self : HeatCapacityPolynomial,
  temperature_k : Double,
) -> Double raise VleError {
  assert_temperature(temperature_k)
  if !self.is_in_range(temperature_k) {
    raise VleError::InvalidParameter(
      "temperature is outside the heat-capacity correlation range",
    )
  }
  let value = polynomial_value(self.coefficients, temperature_k)
  if value <= 0.0 {
    raise VleError::InvalidParameter("heat capacity must be positive")
  }
  value
}

///|
pub fn HeatCapacityPolynomial::derivative(
  self : HeatCapacityPolynomial,
  temperature_k : Double,
) -> Double raise VleError {
  assert_temperature(temperature_k)
  if !self.is_in_range(temperature_k) {
    raise VleError::InvalidParameter(
      "temperature is outside the heat-capacity correlation range",
    )
  }
  if self.coefficients.length() <= 1 {
    return 0.0
  }
  for i = self.coefficients.length() - 1, value = 0.0; i >= 1; {
    continue i - 1, value * temperature_k + i.to_double() * self.coefficients[i]
  } nobreak {
    value
  }
}

///|
pub fn HeatCapacityPolynomial::integral(
  self : HeatCapacityPolynomial,
  low_k~ : Double,
  high_k~ : Double,
) -> Double raise VleError {
  assert_temperature(low_k)
  assert_temperature(high_k)
  if !self.is_in_range(low_k) || !self.is_in_range(high_k) {
    raise VleError::InvalidParameter(
      "integration bounds are outside the heat-capacity correlation range",
    )
  }
  for i = 0, value = 0.0; i < self.coefficients.length(); {
    let exponent = (i + 1).to_double()
    continue i + 1,
      value +
      self.coefficients[i] /
      exponent *
      (@math.pow(high_k, exponent) - @math.pow(low_k, exponent))
  } nobreak {
    value
  }
}

///|
pub fn watts_per_mol_to_joules_per_mol(value : Double) -> Double {
  value
}

///|
pub fn watson_latent_heat(
  reference_latent_heat_j_per_mol~ : Double,
  reference_temperature_k~ : Double,
  critical_temperature_k~ : Double,
  temperature_k~ : Double,
) -> Double raise VleError {
  assert_temperature(reference_temperature_k)
  assert_temperature(critical_temperature_k)
  assert_temperature(temperature_k)
  if reference_latent_heat_j_per_mol <= 0.0 {
    raise VleError::InvalidParameter("reference latent heat must be positive")
  }
  if reference_temperature_k >= critical_temperature_k ||
    temperature_k >= critical_temperature_k {
    raise VleError::InvalidParameter(
      "Watson correlation requires temperatures below the critical point",
    )
  }
  let reference_reduced = 1.0 - reference_temperature_k / critical_temperature_k
  let reduced = 1.0 - temperature_k / critical_temperature_k
  reference_latent_heat_j_per_mol * @math.pow(reduced / reference_reduced, 0.38)
}

///|
pub fn linear_interpolate(
  x0 : Double,
  y0 : Double,
  x1 : Double,
  y1 : Double,
  x : Double,
) -> Double raise VleError {
  if x1 == x0 {
    raise VleError::InvalidParameter("interpolation abscissas must differ")
  }
  y0 + (y1 - y0) * (x - x0) / (x1 - x0)
}

///|
pub fn integrate_trapezoid(
  xs : Array[Double],
  ys : Array[Double],
) -> Double raise VleError {
  assert_same_length(xs.length(), ys.length())
  if xs.length() < 2 {
    raise VleError::InvalidParameter("trapezoid integration needs two points")
  }
  for i = 0, value = 0.0; i + 1 < xs.length(); {
    let width = xs[i + 1] - xs[i]
    if width < 0.0 {
      raise VleError::InvalidParameter("integration abscissas must be ordered")
    }
    continue i + 1, value + width * (ys[i] + ys[i + 1]) / 2.0
  } nobreak {
    value
  }
}

///|
pub fn sample_polynomial(
  correlation : HeatCapacityPolynomial,
  low_k~ : Double,
  high_k~ : Double,
  count~ : Int,
) -> Array[Double] raise VleError {
  if count < 2 {
    raise VleError::InvalidParameter("a sample needs at least two points")
  }
  if high_k <= low_k {
    raise VleError::InvalidRange(low=low_k, high=high_k)
  }
  let step = (high_k - low_k) / (count - 1).to_double()
  [
    for i in 0.. correlation.evaluate(low_k + i.to_double() * step)
  ]
}

///|
pub fn average_heat_capacity(
  correlation : HeatCapacityPolynomial,
  low_k~ : Double,
  high_k~ : Double,
) -> Double raise VleError {
  if high_k <= low_k {
    raise VleError::InvalidRange(low=low_k, high=high_k)
  }
  correlation.integral(low_k~, high_k~) / (high_k - low_k)
}

///|
pub fn heat_capacity_change(
  correlation : HeatCapacityPolynomial,
  low_k~ : Double,
  high_k~ : Double,
) -> Double raise VleError {
  correlation.evaluate(high_k) - correlation.evaluate(low_k)
}

///|
pub fn component_liquid_enthalpy(
  property : ThermoProperty,
  temperature_k : Double,
) -> Double raise VleError {
  property.liquid_heat_capacity.integral(
    low_k=property.reference_temperature_k,
    high_k=temperature_k,
  )
}

///|
pub fn component_vapor_enthalpy(
  property : ThermoProperty,
  temperature_k : Double,
) -> Double raise VleError {
  let sensible = property.vapor_heat_capacity.integral(
    low_k=property.normal_boiling_temperature_k,
    high_k=temperature_k,
  )
  let liquid_to_boiling = component_liquid_enthalpy(
    property,
    property.normal_boiling_temperature_k,
  )
  liquid_to_boiling + property.latent_heat_j_per_mol + sensible
}

///|
pub fn component_phase_enthalpy(
  property : ThermoProperty,
  temperature_k : Double,
  phase : PhaseKind,
) -> Double raise VleError {
  match phase {
    Liquid => component_liquid_enthalpy(property, temperature_k)
    Vapor => component_vapor_enthalpy(property, temperature_k)
    Solid => component_liquid_enthalpy(property, temperature_k)
  }
}

///|
pub fn liquid_sensible_heat_capacity(
  property : ThermoProperty,
  temperature_k : Double,
) -> Double raise VleError {
  property.liquid_heat_capacity.evaluate(temperature_k)
}

///|
pub fn vapor_sensible_heat_capacity(
  property : ThermoProperty,
  temperature_k : Double,
) -> Double raise VleError {
  property.vapor_heat_capacity.evaluate(temperature_k)
}

///|
pub fn latent_heat_at(
  property : ThermoProperty,
  temperature_k : Double,
  critical_temperature_k : Double,
) -> Double raise VleError {
  watson_latent_heat(
    reference_latent_heat_j_per_mol=property.latent_heat_j_per_mol,
    reference_temperature_k=property.normal_boiling_temperature_k,
    critical_temperature_k~,
    temperature_k~,
  )
}