///|
pub struct WilsonBinary {
  lambda12 : Double
  lambda21 : Double
  source : String
} derive(Debug)

///|
pub fn WilsonBinary::new(
  lambda12~ : Double,
  lambda21~ : Double,
  source~ : String,
) -> WilsonBinary {
  { lambda12, lambda21, source }
}

///|
pub struct NrtlBinary {
  tau12 : Double
  tau21 : Double
  alpha : Double
  source : String
} derive(Debug)

///|
pub fn NrtlBinary::new(
  tau12~ : Double,
  tau21~ : Double,
  alpha~ : Double,
  source~ : String,
) -> NrtlBinary {
  { tau12, tau21, alpha, source }
}

///|
pub struct UniquacComponent {
  r : Double
  q : Double
} derive(Debug)

///|
pub fn UniquacComponent::new(r~ : Double, q~ : Double) -> UniquacComponent {
  { r, q }
}

///|
pub struct UniquacBinary {
  comp1 : UniquacComponent
  comp2 : UniquacComponent
  tau12 : Double
  tau21 : Double
  source : String
} derive(Debug)

///|
pub fn UniquacBinary::new(
  comp1~ : UniquacComponent,
  comp2~ : UniquacComponent,
  tau12~ : Double,
  tau21~ : Double,
  source~ : String,
) -> UniquacBinary {
  { comp1, comp2, tau12, tau21, source }
}

///|
pub(all) enum ActivityModel {
  Ideal
  Wilson(WilsonBinary)
  Nrtl(NrtlBinary)
  Uniquac(UniquacBinary)
} derive(Debug)

///|
pub fn ActivityModel::ideal() -> ActivityModel {
  ActivityModel::Ideal
}

///|
pub struct WilsonMatrix {
  lambda : Array[Array[Double]]
  source : String
} derive(Debug)

///|
pub fn WilsonMatrix::new(
  lambda~ : Array[Array[Double]],
  source~ : String,
) -> WilsonMatrix {
  { lambda, source }
}

///|
pub struct NrtlMatrix {
  tau : Array[Array[Double]]
  alpha : Array[Array[Double]]
  source : String
} derive(Debug)

///|
pub fn NrtlMatrix::new(
  tau~ : Array[Array[Double]],
  alpha~ : Array[Array[Double]],
  source~ : String,
) -> NrtlMatrix {
  { tau, alpha, source }
}

///|
pub struct UniquacMatrix {
  r : Array[Double]
  q : Array[Double]
  tau : Array[Array[Double]]
  source : String
} derive(Debug)

///|
pub fn UniquacMatrix::new(
  r~ : Array[Double],
  q~ : Array[Double],
  tau~ : Array[Array[Double]],
  source~ : String,
) -> UniquacMatrix {
  { r, q, tau, source }
}

///|
pub(all) enum MulticomponentActivityModel {
  Ideal
  Wilson(WilsonMatrix)
  Nrtl(NrtlMatrix)
  Uniquac(UniquacMatrix)
} derive(Debug)

///|
pub fn multicomponent_activity_coefficients(
  model : MulticomponentActivityModel,
  liquid : Array[Double],
  temperature_k : Double,
) -> Array[Double] raise VleError {
  let x = normalize(liquid)
  assert_temperature(temperature_k)
  match model {
    Ideal => [ for _ in x => 1.0 ]
    Wilson(params) => wilson_matrix_gamma(params, x)
    Nrtl(params) => nrtl_matrix_gamma(params, x)
    Uniquac(params) => uniquac_matrix_gamma(params, x)
  }
}

///|
fn assert_square(matrix : Array[Array[Double]], n : Int) -> Unit raise VleError {
  if matrix.length() != n {
    raise VleError::LengthMismatch(expected=n, actual=matrix.length())
  }
  for row in matrix {
    if row.length() != n {
      raise VleError::LengthMismatch(expected=n, actual=row.length())
    }
  }
}

///|
fn wilson_matrix_gamma(
  params : WilsonMatrix,
  x : Array[Double],
) -> Array[Double] raise VleError {
  let n = x.length()
  assert_square(params.lambda, n)
  [
    for i in 0.. {
      let sum_i = for j = 0, acc = 0.0; j < n; {
        continue j + 1, acc + x[j] * clamp_positive(params.lambda[i][j])
      } nobreak {
        acc
      }
      let correction = for j = 0, acc = 0.0; j < n; {
        let sum_j = for k = 0, inner = 0.0; k < n; {
          continue k + 1, inner + x[k] * clamp_positive(params.lambda[j][k])
        } nobreak {
          inner
        }
        continue j + 1, acc + x[j] * clamp_positive(params.lambda[j][i]) / sum_j
      } nobreak {
        acc
      }
      @math.exp(1.0 - @math.ln(sum_i) - correction)
    }
  ]
}

///|
fn nrtl_matrix_gamma(
  params : NrtlMatrix,
  x : Array[Double],
) -> Array[Double] raise VleError {
  let n = x.length()
  assert_square(params.tau, n)
  assert_square(params.alpha, n)
  let g = [
    for i in 0.. {
      [
        for j in 0.. @math.exp(0.0 - params.alpha[i][j] * params.tau[i][j])
      ]
    }
  ]
  [
    for i in 0.. {
      let ln_gamma = for j = 0, acc = 0.0; j < n; {
        let denominator_j = for k = 0, inner = 0.0; k < n; {
          continue k + 1, inner + x[k] * g[k][j]
        } nobreak {
          inner
        }
        let weighted_tau = for k = 0, inner = 0.0; k < n; {
          continue k + 1, inner + x[k] * params.tau[k][j] * g[k][j]
        } nobreak {
          inner
        }
        let curvature = params.tau[i][j] - weighted_tau / denominator_j
        continue j + 1, acc + x[j] * g[i][j] / denominator_j * curvature
      } nobreak {
        acc
      }
      @math.exp(ln_gamma)
    }
  ]
}

///|
fn uniquac_matrix_gamma(
  params : UniquacMatrix,
  x : Array[Double],
) -> Array[Double] raise VleError {
  let n = x.length()
  assert_same_length(n, params.r.length())
  assert_same_length(n, params.q.length())
  assert_square(params.tau, n)
  let z = 10.0
  let r_total = for i = 0, acc = 0.0; i < n; {
    continue i + 1, acc + x[i] * clamp_positive(params.r[i])
  } nobreak {
    acc
  }
  let q_total = for i = 0, acc = 0.0; i < n; {
    continue i + 1, acc + x[i] * clamp_positive(params.q[i])
  } nobreak {
    acc
  }
  let phi = [ for i in 0.. x[i] * params.r[i] / r_total ]
  let theta = [ for i in 0.. x[i] * params.q[i] / q_total ]
  let l = [
    for i in 0.. {
      z / 2.0 * (params.r[i] - params.q[i]) - (params.r[i] - 1.0)
    }
  ]
  let l_mix = for i = 0, acc = 0.0; i < n; {
    continue i + 1, acc + x[i] * l[i]
  } nobreak {
    acc
  }
  [
    for i in 0.. {
      let combinatorial = @math.ln(phi[i] / x[i]) +
        z / 2.0 * params.q[i] * @math.ln(theta[i] / phi[i]) +
        l[i] -
        phi[i] / x[i] * l_mix
      let residual = for j = 0, acc = 0.0; j < n; {
        let denominator_j = for k = 0, inner = 0.0; k < n; {
          continue k + 1, inner + theta[k] * params.tau[k][j]
        } nobreak {
          inner
        }
        let second = for k = 0, inner = 0.0; k < n; {
          continue k + 1,
            inner +
            theta[k] *
            params.tau[j][k] /
            (for m = 0, nested = 0.0; m < n; {
              continue m + 1, nested + theta[m] * params.tau[m][k]
            } nobreak {
              nested
            })
        } nobreak {
          inner
        }
        continue j + 1,
          acc +
          params.q[i] * (theta[j] * params.tau[j][i] / denominator_j - second)
      } nobreak {
        acc
      }
      @math.exp(combinatorial + residual)
    }
  ]
}

///|
pub fn activity_coefficients(
  model : ActivityModel,
  liquid : Array[Double],
  temperature_k : Double,
) -> Array[Double] raise VleError {
  let x = normalize(liquid)
  assert_temperature(temperature_k)
  match model {
    Ideal => [ for _ in x => 1.0 ]
    Wilson(params) => wilson_binary_gamma(params, x)
    Nrtl(params) => nrtl_binary_gamma(params, x)
    Uniquac(params) => uniquac_binary_gamma(params, x)
  }
}

///|
fn ensure_binary(x : Array[Double]) -> Unit raise VleError {
  assert_same_length(2, x.length())
}

///|
fn wilson_binary_gamma(
  params : WilsonBinary,
  x : Array[Double],
) -> Array[Double] raise VleError {
  ensure_binary(x)
  let x1 = clamp_positive(x[0])
  let x2 = clamp_positive(x[1])
  let l12 = clamp_positive(params.lambda12)
  let l21 = clamp_positive(params.lambda21)
  let a = x1 + l12 * x2
  let b = x2 + l21 * x1
  let ln1 = 0.0 - @math.ln(a) + x2 * (l12 / a - l21 / b)
  let ln2 = 0.0 - @math.ln(b) - x1 * (l12 / a - l21 / b)
  [@math.exp(ln1), @math.exp(ln2)]
}

///|
fn nrtl_binary_gamma(
  params : NrtlBinary,
  x : Array[Double],
) -> Array[Double] raise VleError {
  ensure_binary(x)
  let x1 = clamp_positive(x[0])
  let x2 = clamp_positive(x[1])
  let g12 = @math.exp(0.0 - params.alpha * params.tau12)
  let g21 = @math.exp(0.0 - params.alpha * params.tau21)
  let d1 = x1 + x2 * g21
  let d2 = x2 + x1 * g12
  let ln1 = x2 *
    x2 *
    (params.tau21 * g21 * g21 / (d1 * d1) + params.tau12 * g12 / (d2 * d2))
  let ln2 = x1 *
    x1 *
    (params.tau12 * g12 * g12 / (d2 * d2) + params.tau21 * g21 / (d1 * d1))
  [@math.exp(ln1), @math.exp(ln2)]
}

///|
fn uniquac_binary_gamma(
  params : UniquacBinary,
  x : Array[Double],
) -> Array[Double] raise VleError {
  ensure_binary(x)
  let z = 10.0
  let x1 = clamp_positive(x[0])
  let x2 = clamp_positive(x[1])
  let r1 = clamp_positive(params.comp1.r)
  let r2 = clamp_positive(params.comp2.r)
  let q1 = clamp_positive(params.comp1.q)
  let q2 = clamp_positive(params.comp2.q)
  let phi1 = r1 * x1 / (r1 * x1 + r2 * x2)
  let phi2 = r2 * x2 / (r1 * x1 + r2 * x2)
  let theta1 = q1 * x1 / (q1 * x1 + q2 * x2)
  let theta2 = q2 * x2 / (q1 * x1 + q2 * x2)
  let l1 = z / 2.0 * (r1 - q1) - (r1 - 1.0)
  let l2 = z / 2.0 * (r2 - q2) - (r2 - 1.0)
  let combin1 = @math.ln(phi1 / x1) +
    z / 2.0 * q1 * @math.ln(theta1 / phi1) +
    l1 -
    phi1 / x1 * (x1 * l1 + x2 * l2)
  let combin2 = @math.ln(phi2 / x2) +
    z / 2.0 * q2 * @math.ln(theta2 / phi2) +
    l2 -
    phi2 / x2 * (x1 * l1 + x2 * l2)
  let res1 = q1 *
    (
      1.0 -
      @math.ln(theta1 + theta2 * params.tau21) -
      (
        theta1 / (theta1 + theta2 * params.tau21) +
        theta2 * params.tau12 / (theta2 + theta1 * params.tau12)
      )
    )
  let res2 = q2 *
    (
      1.0 -
      @math.ln(theta2 + theta1 * params.tau12) -
      (
        theta2 / (theta2 + theta1 * params.tau12) +
        theta1 * params.tau21 / (theta1 + theta2 * params.tau21)
      )
    )
  [@math.exp(combin1 + res1), @math.exp(combin2 + res2)]
}