{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2ab7df80",
   "metadata": {},
   "source": [
    "# LKP (Lee-Kesler-Plöcker)\n",
    "\n",
    "The LKP model is a sort of hybrid between corresponding states and multiparameter EOS, simple EOS are developed for a reference fluid, and a simple fluid, and the acentric factor of the mixture is used to weight the two.\n",
    "\n",
    "The reduced residual Helmholtz energy for the mixture is evaluated from\n",
    "\n",
    "$$ \\alpha^{\\rm r} = \\left(1-\\frac{\\omega_{\\rm mix}}{\\omega_{\\rm ref}}\\right)\\alpha^{\\rm r}_{\\rm simple} + \\frac{\\omega_{\\rm mix}}{\\omega_{\\rm ref}}\\alpha^{\\rm r}_{\\rm ref} $$\n",
    "\n",
    "where the contributions are each of the form\n",
    "\n",
    "$$\\alpha^{\\rm r}_{X}(\\tau, \\delta) = B\\left(\\frac{\\delta}{Z_c}\\right) + \\frac{C}{2}\\left(\\frac{\\delta}{Z_c}\\right)^2 + \\frac{D}{5}\\left(\\frac{\\delta}{Z_c}\\right)^5 - \\frac{c_4\\tau^3}{2\\gamma}\\left(\\gamma\\left(\\frac{\\delta}{Z_c}\\right)^2+\\beta+1\\right)\\exp\\left(-\\gamma\\left(\\frac{\\delta}{Z_c}\\right)^2\\right) + \\frac{c_4\\tau^3}{2\\gamma}(\\beta+1) $$\n",
    "\n",
    "where $X$ is one of simple or reference (abbreviation: ref) with the matching sets of coefficients taken from this table:\n",
    "\n",
    "|   var     |    simple   |   reference  |\n",
    "|-----------|-------------|--------------|\n",
    "|  $b_1$    | 0.1181193   |  0.2026579   |\n",
    "|  $b_2$    | 0.265728    |  0.331511    |\n",
    "|  $b_3$    | 0.154790    |  0.276550e-1 |\n",
    "|  $b_4$    | 0.303230e-1 |  0.203488    |\n",
    "|  $c_1$    | 0.236744e-1 |  0.313385e-1 |\n",
    "|  $c_2$    | 0.186984e-1 |  0.503618e-1 |\n",
    "|  $c_3$    | 0           |  0.169010e-1 |\n",
    "|  $c_4$    | 0.427240e-1 |  0.41577e-1  |\n",
    "|  $d_1$    | 0.155428e-4 |  0.487360e-4 |\n",
    "|  $d_2$    | 0.623689e-4 |  0.740336e-5 |\n",
    "|  $\\beta$  | 0.653920    |  1.226       |\n",
    "|  $\\gamma$ | 0.601670e-1 |  0.03754     |\n",
    "|  $\\omega$ | 0.0         |  0.3978      |\n",
    "\n",
    "The terms in the contributions are given by:\n",
    "\n",
    "$$ B = b_1 - b_2\\tau - b_3\\tau^2 - b_4\\tau^3$$\n",
    "$$ C = c_1 - c_2\\tau + c_3\\tau^3$$\n",
    "$$ D = d_1 + d_2\\tau$$\n",
    "\n",
    "For density, the reduced density $\\delta$ is defined by \n",
    "$$\\delta = \\frac{\\rho}{\\rho_{\\rm red}} = v_{\\rm c,mix}\\rho$$\n",
    "in which the reducing density is the reciprocal of the pseudo-critical volume obtained from\n",
    "$$ v_{\\rm c, mix} = \\sum_{i=1}^{N-1}\\sum_{j=i+1}^N x_ix_jv_{ij} $$\n",
    "$$ v_{c,ij} = \\frac{1}{8}(v_{c,i}^{1/3} + v_{c,j}^{1/3})^3$$\n",
    "and the critical volumes are estimated from\n",
    "$$ v_{c,i} = (0.2905-0.085\\omega_i)\\frac{RT_{c,i}}{p_{c,i}} $$\n",
    "\n",
    "For temperature, the reciprocal reduced density is defined by\n",
    "$$ \\tau = \\frac{T_{\\rm c,mix}}{T}$$\n",
    "with\n",
    "$$ T_{\\rm c,mix} = \\frac{1}{v_{c,mix}^{\\eta}}\\sum_{i=1}^{N-1}\\sum_{j=i+1}^{N}x_ix_jv_{c,ij}^{\\eta}T_{c,ij}$$\n",
    "with $\\eta=0.25$ and \n",
    "$$T_{c,ij} = k_{ij}\\sqrt{T_{c,i}T_{c,j}}$$\n",
    "\n",
    "Note: the default interaction parameter $k_{ij}$ is therefore 1, rather than 0 in the case of SAFT and cubic models.\n",
    "\n",
    "Finally the parameter $Z_c$ is defined by\n",
    "$$ Z_c = 0.2905-0.085\\omega_{\\rm mix} $$\n",
    "with the mixture acentric factor defined by \n",
    "$$ \\omega_{\\rm mix} = \\sum_i x_i\\omega_i$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a582cd08",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:41:23.493040Z",
     "iopub.status.busy": "2024-03-15T22:41:23.492861Z",
     "iopub.status.idle": "2024-03-15T22:41:23.562372Z",
     "shell.execute_reply": "2024-03-15T22:41:23.561844Z"
    }
   },
   "outputs": [],
   "source": [
    "import teqp, numpy as np\n",
    "spec = {\n",
    "    \"Tcrit / K\": [190.564, 126.192],\n",
    "    \"pcrit / Pa\": [4.5992e6, 3.3958e6],\n",
    "    \"acentric\": [0.011, 0.037],\n",
    "    \"R / J/mol/K\": 8.3144598,\n",
    "    \"kmat\": [[1.0, 0.977],[0.977, 1.0]]\n",
    "}\n",
    "model = teqp.make_model({'kind': 'LKP', 'model': spec}, validate=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "89c910ff",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-03-15T22:41:23.566309Z",
     "iopub.status.busy": "2024-03-15T22:41:23.565375Z",
     "iopub.status.idle": "2024-03-15T22:41:23.569613Z",
     "shell.execute_reply": "2024-03-15T22:41:23.569171Z"
    }
   },
   "outputs": [],
   "source": [
    "# A little sanity check, with the check value from TREND\n",
    "expected = -0.18568096994998817\n",
    "diff = abs(model.get_Ar00(300, 8000.1, np.array([0.8, 0.2])) - expected)\n",
    "assert(diff < 1e-13)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}