{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fec96515",
   "metadata": {},
   "source": [
    "# VLE from fugacity calculation\n",
    "\n",
    "Use the Generalized Peng-Robinson equation of state to calculate the vapor-liquid equilibrium curve. We will need to evaluate the Peng-Robinson EOS and use the Peng-Robinson form of the fugacity equation.\n",
    "\n",
    "This notebook shows how I developed preliminary calcuations before generating the final plots. \n",
    "\n",
    "Eric Furst, December 2025"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84f5edf0",
   "metadata": {},
   "source": [
    "## Generalized Peng-Robinson Equation of State\n",
    "\n",
    "Routines to calcualte the Generalized Peng-Robinson Equation of State\n",
    "\n",
    "SIS is Stanley I. Sandler, *Chemical, Biochemcial and Engineering Thermodynamics*, 5th ed.\n",
    "\n",
    "Eric Furst, November 2025\n",
    "\n",
    "The Generalized Peng-Robinson equation of state is\n",
    "\n",
    "$$ P = \\frac{RT}{\\underline{V}-b} - \\frac{a(T)}{\\underline{V}(\\underline{V}+b) + b(\\underline{V}-b)} \\tag{Eq. 6.4-2}$$\n",
    "\n",
    "with\n",
    "\n",
    "$$ b = 0.07780 \\frac{RT_c}{P_c} \\tag{Eq. 6.7-2}$$\n",
    "$$ a(T) = a(T_c)\\alpha(T) = 0.45724 \\frac{R^2T_c^2}{P_c}\\alpha(T) \\tag{Eq. 6.7-1}$$\n",
    "$$ \\sqrt{\\alpha} = 1 + \\kappa \\left ( 1- \\sqrt{\\frac{T}{T_c}} \\right ) \\tag{Eq. 6.7-3}$$\n",
    "$$ \\kappa = 0.37464 + 1.54226\\omega − 0.26992\\omega^2 \\tag{Eq. 6.7-4}$$\n",
    "\n",
    "The acentric factor $\\omega$ and the crticial temperatures and pressures are given in SIS table 6.6-1.\n",
    "\n",
    "Calculating the pressure $P$ given $\\underline{V}$ and $T$ is straightforward, but to calculate the molar volume given $P$ and $T$, we need to solve the cubic equation of state of the form\n",
    "\n",
    "$$ Z^3 + \\alpha Z^2 + \\beta Z + \\gamma = 0 \\tag{Eq. 6.4-4}$$\n",
    "\n",
    "where $Z$ is the compressibility factor\n",
    "\n",
    "$$ Z = \\frac{P \\underbar{V{}}}{RT} $$\n",
    "\n",
    "For the Peng-Robinson EOS (see SIS Table 6.4-3),\n",
    "\n",
    "$$ \\alpha = -1 + B $$\n",
    "$$ \\beta = A - 3B^2 -2B $$\n",
    "$$ \\gamma = -AB + B^2 + B^3 $$\n",
    "\n",
    "and \n",
    "\n",
    "$$ A = \\frac{aP}{(RT)^2} $$\n",
    "$$ B = \\frac{bP}{RT} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8cd6123e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Generalized Peng-Robinson EOS \n",
    "PR_P returns the pressure given V, T, Pc, Tc, omega\n",
    "PR_Z returns the compressibility factor for all real roots of EOS given P, T, Pc, Tc, omega\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "from scipy import constants\n",
    "from numpy.polynomial import Polynomial\n",
    "\n",
    "R = constants.R # Set the gas constant to R\n",
    "\n",
    "def calc_b(Pc,Tc):\n",
    "    return 0.07780*R*Tc/Pc\n",
    "\n",
    "def calc_a(T,Pc,Tc,omega):\n",
    "    kappa = 0.37464 + 1.54226*omega - 0.26992*omega**2\n",
    "    sqrtalpha = 1 + kappa*(1-np.sqrt(T/Tc))\n",
    "    return 0.45724*R**2*Tc**2/Pc*sqrtalpha**2\n",
    "\n",
    "# Calculate the presure given V, T for PR EOS\n",
    "def PR_P(V,T,Pc,Tc,omega):\n",
    "    a = calc_a(T,Pc,Tc,omega)\n",
    "    b = calc_b(Pc,Tc)\n",
    "\n",
    "    P = R*T/(V-b) - a/(V*(V+b)+b*(V-b))\n",
    "    return P\n",
    "\n",
    "# Calculate the compressibility coefficient given P, T for PR EOS\n",
    "# Note that we can return multiple real roots (up to three)\n",
    "# The largest and smallest will be the vapor and liquid, respectively\n",
    "def PR_Z(P, T, Pc, Tc, omega):\n",
    "    # Calculate a, b, A, and B\n",
    "    a = calc_a(T, Pc, Tc, omega)\n",
    "    b = calc_b(Pc, Tc)\n",
    "    A = a*P/R**2/T**2\n",
    "    B = b*P/R/T\n",
    "\n",
    "    # Definitions of alpha, beta, gamma in SIS Table 6.4-3 for PR EOS\n",
    "    alpha = -1 + B\n",
    "    beta = A - 3*B**2 -2*B\n",
    "    gamma = -A*B + B**2 + B**3\n",
    "\n",
    "    # polynomial with coefficients in increasing order: c0 + c1 x + c2 x**2 + ...\n",
    "    p = Polynomial([ gamma, beta, alpha, 1 ])  \n",
    "\n",
    "    roots = p.roots()        # returns all (possibly complex) roots of Z\n",
    "    real_roots = roots.real[abs(roots.imag) < 1e-12]  # filter real ones\n",
    "\n",
    "    return real_roots\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c3f3877",
   "metadata": {},
   "source": [
    "## Fugacity coefficient for the Peng-Robinson equation of state\n",
    "\n",
    "The fugacity coefficient $ f/P $ is written as\n",
    "\n",
    "$$ \\ln \\frac{f}{P} = (Z - 1) -\\ln (Z-B) - \\frac{A}{2\\sqrt{2}B}\\ln \\left [ \\frac{Z + (1+\\sqrt{2})B}{Z + (1-\\sqrt{2})B} \\right ] $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ A = \\frac{aP}{(RT)^2} $$\n",
    "$$ B = \\frac{bP}{RT} $$\n",
    "\n",
    "and \n",
    "\n",
    "$$ Z = \\frac{P \\underbar{V{}}}{RT} $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "32f3e31d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Fugacity coefficient of the Peng-Robinson EOS\n",
    "Returns the natural log of the fugacity coefficient\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "\n",
    "def fugacity_PR(Z,A,B):\n",
    "    return (Z-1) - np.log(Z-B)-A/2/np.sqrt(2)/B*np.log((Z+(1+np.sqrt(2))*B)/(Z+(1-np.sqrt(2))*B))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9f11582",
   "metadata": {},
   "source": [
    "## Example calculation - VLE for nitrogen"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef603cb4",
   "metadata": {},
   "source": [
    "Start with data for nitrogen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bd8d2ff0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data for nitrogen\n",
    "# Molecular weight\n",
    "MW = 28.0134 # g/mol\n",
    "\n",
    "# Critical parameters and acentricity\n",
    "Pc = 33.978*100000 # Critical pressure in Pa\n",
    "Tc = 126.19 # Critical temp in K\n",
    "Vc = 1/(11.18*1000) # Critical volume in m^3/mol (original is 11.18 mol/l)\n",
    "omega = 0.040 # acentric factor\n",
    "\n",
    "# Triple point\n",
    "T_triple = 63.14 # temp in K\n",
    "P_triple = 0.1252*100000 # pressure in Pa"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49fa42eb",
   "metadata": {},
   "source": [
    "The next block is a minimal calculation. We guess the pressure and check to see if the equlibrium condition is satisfied. You can iterate by hand to find the correct $P^\\mathrm{sat}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "42d9332d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8275367561074513 0.8344346518047008 -0.006897895697249523\n",
      "Psat(101 K) = 840000 Pa\n"
     ]
    }
   ],
   "source": [
    "# This will be a manual guess and check\n",
    "import numpy as np\n",
    "import scipy.constants as constants\n",
    "\n",
    "R = constants.R\n",
    "\n",
    "# Specifiy the temperature\n",
    "T = 0.8*Tc\n",
    "\n",
    "# Guess a pressure (in Pa)\n",
    "P = 840000 \n",
    "\n",
    "# Calculate liquid and vapor fugacities\n",
    "A = calc_a(T,Pc,Tc,omega)*P/R**2/T**2\n",
    "B = calc_b(Pc,Tc)*P/R/T\n",
    "Z = PR_Z(P,T,Pc,Tc,omega) # This returns all roots\n",
    "\n",
    "ln_f_L = fugacity_PR(np.min(Z),A,B) # ln of fugacity coefficient for liquid\n",
    "ln_f_V = fugacity_PR(np.max(Z),A,B) # ln of fugacity coefficient for liquid\n",
    "\n",
    "# Are they equal (to within a tolerance?)\n",
    "print(np.exp(ln_f_L), np.exp(ln_f_V), np.exp(ln_f_L)-np.exp(ln_f_V))\n",
    "\n",
    "print(f\"Psat({T:.0f} K) = {P:.0f} Pa\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41f8639e",
   "metadata": {},
   "source": [
    "Next, try implementing an automatic search for $P^\\mathrm{sat}(T)$ by calculating the residual of the fugacity equation\n",
    "\n",
    "$$\\ln \\frac{f^L}{P} - \\ln \\frac{f^V}{P} = 0$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d52a581e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a function that returns the residual of the fugacity calculation\n",
    "def residual_lnphi(P):\n",
    "    \"\"\"\n",
    "    Returns ln(phi_L) - ln(phi_V) at given P (Pa).\n",
    "    At Psat, this should be zero.\n",
    "    As written, requires globals Tc, Pc, omega, and T\n",
    "    \"\"\"\n",
    "    # EOS parameters\n",
    "    A = calc_a(T, Pc, Tc, omega)*P/R**2/T**2\n",
    "    B = calc_b(Pc, Tc)*P/R/T\n",
    "\n",
    "    # All compressibility roots\n",
    "    Z_roots = np.array(PR_Z(P, T, Pc, Tc, omega))\n",
    "\n",
    "    # Keep real roots only (in case your solver returns complex with tiny imag)\n",
    "    Z_real = Z_roots[np.isreal(Z_roots)].real\n",
    "\n",
    "    ZL = np.min(Z_real)  # liquid root (smallest Z)\n",
    "    ZV = np.max(Z_real)  # vapor root (largest Z)\n",
    "\n",
    "    ln_f_L = fugacity_PR(ZL, A, B)  # ln(phi_L)\n",
    "    ln_f_V = fugacity_PR(ZV, A, B)  # ln(phi_V)\n",
    "\n",
    "    return ln_f_L - ln_f_V\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "294094c3",
   "metadata": {},
   "source": [
    "We'll use Newton's method in the scipy.optimize library to minimize the residual. We need a decent guess for the pressure. If it is too far away from the final value of $P^\\mathrm{sat}$, the routine will fail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "05077a86",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Psat(101 K) = 830968 Pa\n",
      "Residual at Psat = 1.3322676295501878e-15\n"
     ]
    }
   ],
   "source": [
    "# Example run\n",
    "from scipy.optimize import newton\n",
    "\n",
    "P_guess = 200000  # Pa\n",
    "T = 0.8*Tc\n",
    "Psat = newton(residual_lnphi, P_guess)  # or newton(residual_lnphi, P0, P1) for secant\n",
    "print(f\"Psat({T:.0f} K) = {Psat:.0f} Pa\")\n",
    "print(\"Residual at Psat =\", residual_lnphi(Psat))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b9f8dbe",
   "metadata": {},
   "source": [
    "That looks like it works. Next, let's write a routine to generate a Pandas dataframe of the temperature, pressure, molar volume of the liquid and molar volume of the vapor."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9414e99",
   "metadata": {},
   "source": [
    "### Create dataframe of VLE values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "62364367",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           T             P        vL        vV\n",
      "0  63.140000  12830.310618  0.000028  0.040580\n",
      "1  63.171351  12903.335640  0.000028  0.040369\n",
      "2  63.202703  12976.696886  0.000029  0.040159\n",
      "3  63.234054  13050.395494  0.000029  0.039950\n",
      "4  63.265406  13124.432605  0.000029  0.039743\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# Example temperature grid (pick your own range)\n",
    "T_min = T_triple\n",
    "T_max = 0.997 * Tc\n",
    "n_T   = 2000\n",
    "\n",
    "T_list = np.linspace(T_min, T_max, n_T)\n",
    "\n",
    "rows = []\n",
    "\n",
    "# Initial pressure guess: something reasonable, e.g. 0.1*Pc or from your old calc\n",
    "P_guess = 10000\n",
    "\n",
    "for T in T_list:\n",
    "    Psat = newton(residual_lnphi, P_guess)\n",
    "    \n",
    "    # All compressibility roots\n",
    "    Z_roots = np.array(PR_Z(Psat, T, Pc, Tc, omega))\n",
    "\n",
    "    # Keep real roots only (in case your solver returns complex with tiny imag)\n",
    "    Z_real = Z_roots[np.isreal(Z_roots)].real\n",
    "\n",
    "    ZL = np.min(Z_real)  # liquid root of compressibility factor (smallest Z)\n",
    "    ZV = np.max(Z_real)  # vapor root of compressibility factor (largest Z)\n",
    "\n",
    "    vL = ZL*R*T/Psat\n",
    "    vV = ZV*R*T/Psat\n",
    "\n",
    "    rows.append({\n",
    "        \"T\": T,\n",
    "        \"P\": Psat,\n",
    "        \"vL\": vL,\n",
    "        \"vV\": vV,\n",
    "    })\n",
    "\n",
    "    # Use current Psat as the next initial guess (smooth VLE curve)\n",
    "    P_guess = Psat\n",
    "\n",
    "df_eq = pd.DataFrame(rows)\n",
    "print(df_eq.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2d94860",
   "metadata": {},
   "source": [
    "Plot the P-T diagram..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "346d9640",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35364 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams[\"text.usetex\"] = True\n",
    "\n",
    "# Antoine equation given parameters A, B, and C\n",
    "# Pressure is returned in Pa\n",
    "def P_antoine(A, B, C, T):\n",
    "    return 10**(A - B/(T + C))*100000\n",
    "\n",
    "# Antoine equation parameters\n",
    "# Valid temperature range in K\n",
    "T_Antoine_min = 63.14\n",
    "T_Antoine_max = 126\n",
    "A = 3.7362\n",
    "B = 264.651\n",
    "C = -6.788\n",
    "\n",
    "# Plot the Antoine equation\n",
    "T = np.arange(T_Antoine_min, T_Antoine_max+1, 1)\n",
    "plt.plot(T, P_antoine(A, B, C, T), \"c\", label=r\"$P^\\mathrm{vap}$ (Antoine)\")\n",
    "\n",
    "sns.lineplot(data=df_eq, x=\"T\", y=\"P\", label=r\"$P^\\mathrm{vap}$ (PR EOS)\")\n",
    "plt.xlabel(\"T (K)\")\n",
    "plt.ylabel(\"P (Pa)\")\n",
    "plt.yscale(\"log\")\n",
    "plt.tight_layout()\n",
    "\n",
    "# Plot the critical point\n",
    "plt.plot(Tc, Pc, \"ro\", label=\"critical point\")\n",
    "\n",
    "# Plot the triple point\n",
    "plt.plot(T_triple, P_triple, \"bo\", label=\"triple point\")\n",
    "\n",
    "# Plot the normal boiling point (from Streng)\n",
    "plt.plot(77.4, 101325, \"mo\", label=\"normal boiling point\")\n",
    "\n",
    "plt.title(r\"$P-T$ diagram for nitrogen, Peng-Robinson EOS\")\n",
    "plt.legend(frameon=True, fancybox=True, edgecolor='none')\n",
    "\n",
    "# Line at 1 bar\n",
    "plt.plot([T_triple, Tc],[100000, 100000], linestyle = '--', color=\"gray\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c53c4fdc",
   "metadata": {},
   "source": [
    "We can also create a P-V diagram using the same dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "68379299",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/vs/zb__jwqr8xj7b90059b80rzr0000gp/T/ipykernel_26142/2492820622.py:55: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.\n",
      "  plt.legend(frameon=True, fancybox=True, edgecolor='none')\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 46276 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This line can cause problems in Google Colab\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams[\"text.usetex\"] = True\n",
    "\n",
    "# Plot (calculated) critical point\n",
    "#Vc = PR_Z(Pc,Tc,Pc,Tc,omega)*R*Tc/Pc\n",
    "#plt.plot(Vc, Pc, \"ro\", label=\"c.p.\") \n",
    "\n",
    "# Add the critical isotherm\n",
    "V = np.linspace(Vc/3,4e-2,10000)\n",
    "#plt.plot(V, PR_P(V,Tc,Pc,Tc,omega), color='gray', linestyle='-',label=f\"critical isotherm, {Tc} K\")\n",
    "plt.plot(V, PR_P(V,Tc,Pc,Tc,omega), color='gray', linestyle='-')\n",
    "\n",
    "# Plot the t.p. tie line\n",
    "row_min_T = df_eq.loc[df_eq[\"T\"].idxmin()]\n",
    "vL = row_min_T[\"vL\"]\n",
    "vV = row_min_T[\"vV\"]\n",
    "P  = row_min_T[\"P\"]\n",
    "plt.plot([vL, vV], [P, P], color='gray', linestyle='-', linewidth=1)\n",
    "\n",
    "# Pick and plot 4 evenly spaced tie lines through the dataframe\n",
    "idx = np.linspace(0, len(df_eq)/4*3, 5, dtype=int)\n",
    "for i in idx:\n",
    "    row = df_eq.iloc[i]\n",
    "    vL, vV, P = row[\"vL\"], row[\"vV\"], row[\"P\"]\n",
    "    \n",
    "    # Draw tie line\n",
    "    plt.plot([vL, vV], [P, P], '-', linewidth=1, color='gray')\n",
    "\n",
    "    # Plot the isotherm above the vapor line\n",
    "    V = np.linspace(row[\"vV\"],4e-2,1000)\n",
    "    plt.plot(V, PR_P(V,row[\"T\"],Pc,Tc,omega), color='gray', linestyle='-', linewidth=1)\n",
    "\n",
    "    # Plot the isotherm below the liquid line\n",
    "    V = np.linspace(row[\"vL\"],2.8e-5,1000)\n",
    "    plt.plot(V, PR_P(V,row[\"T\"],Pc,Tc,omega), color='gray', linestyle='-', linewidth=1)\n",
    "    \n",
    "\n",
    "# Plot the vapor coexistence line\n",
    "#sns.lineplot(data=df_eq, x=\"vV\", y=\"P\", color=\"black\", linewidth=2, label=r\"$\\underbar{V}_\\mathrm{V}$\")\n",
    "sns.lineplot(data=df_eq, x=\"vV\", y=\"P\", color=\"black\", linewidth=2)\n",
    "\n",
    "# Plot the liquid coexistence line\n",
    "#sns.lineplot(data=df_eq, x=\"vL\", y=\"P\", color=\"black\", linewidth=2, label=r\"$\\underbar{V}_\\mathrm{L}$\")\n",
    "sns.lineplot(data=df_eq, x=\"vL\", y=\"P\", color=\"black\", linewidth=2)\n",
    "\n",
    "plt.grid(color='gray', linestyle='--', linewidth=0.5, alpha=0.7)\n",
    "plt.ylim([6e3, 1e7])\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(r\"$\\underbar{V}$ (m$^3$/mol)\")\n",
    "plt.ylabel(r\"$P$ (Pa)\")\n",
    "plt.title('$P-V$ diagram for nitrogen, Peng-Robinson EOS')\n",
    "\n",
    "plt.legend(frameon=True, fancybox=True, edgecolor='none')\n",
    "#plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ad5761c",
   "metadata": {},
   "source": [
    "As we approach the critical point, it becomes more difficult to calculate the binodal. I used many points, which is inefficient, but gets us to 0.99*Pc.\n",
    "\n",
    "Since we use a Pandas dataframe, we can easily generate a table of the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1e8719c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>T</th>\n",
       "      <th>P</th>\n",
       "      <th>vL</th>\n",
       "      <th>vV</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>63.140000</td>\n",
       "      <td>1.283031e+04</td>\n",
       "      <td>0.000028</td>\n",
       "      <td>0.040580</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>63.171351</td>\n",
       "      <td>1.290334e+04</td>\n",
       "      <td>0.000028</td>\n",
       "      <td>0.040369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>63.202703</td>\n",
       "      <td>1.297670e+04</td>\n",
       "      <td>0.000029</td>\n",
       "      <td>0.040159</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>63.234054</td>\n",
       "      <td>1.305040e+04</td>\n",
       "      <td>0.000029</td>\n",
       "      <td>0.039950</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>63.265406</td>\n",
       "      <td>1.312443e+04</td>\n",
       "      <td>0.000029</td>\n",
       "      <td>0.039743</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1995</th>\n",
       "      <td>125.686024</td>\n",
       "      <td>3.319799e+06</td>\n",
       "      <td>0.000079</td>\n",
       "      <td>0.000116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996</th>\n",
       "      <td>125.717376</td>\n",
       "      <td>3.324631e+06</td>\n",
       "      <td>0.000080</td>\n",
       "      <td>0.000115</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1997</th>\n",
       "      <td>125.748727</td>\n",
       "      <td>3.329469e+06</td>\n",
       "      <td>0.000080</td>\n",
       "      <td>0.000114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1998</th>\n",
       "      <td>125.780079</td>\n",
       "      <td>3.334311e+06</td>\n",
       "      <td>0.000081</td>\n",
       "      <td>0.000114</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1999</th>\n",
       "      <td>125.811430</td>\n",
       "      <td>3.339158e+06</td>\n",
       "      <td>0.000081</td>\n",
       "      <td>0.000113</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2000 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "               T             P        vL        vV\n",
       "0      63.140000  1.283031e+04  0.000028  0.040580\n",
       "1      63.171351  1.290334e+04  0.000028  0.040369\n",
       "2      63.202703  1.297670e+04  0.000029  0.040159\n",
       "3      63.234054  1.305040e+04  0.000029  0.039950\n",
       "4      63.265406  1.312443e+04  0.000029  0.039743\n",
       "...          ...           ...       ...       ...\n",
       "1995  125.686024  3.319799e+06  0.000079  0.000116\n",
       "1996  125.717376  3.324631e+06  0.000080  0.000115\n",
       "1997  125.748727  3.329469e+06  0.000080  0.000114\n",
       "1998  125.780079  3.334311e+06  0.000081  0.000114\n",
       "1999  125.811430  3.339158e+06  0.000081  0.000113\n",
       "\n",
       "[2000 rows x 4 columns]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_eq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36bc4abe",
   "metadata": {},
   "source": [
    "Some additional things we can do include \n",
    "\n",
    "1. Plot the triple point isotherm \n",
    "1. Plotting tie lines on the VLE phase envelope and plotting isotherms around the binodal\n",
    "2. Plotting $\\ln P^\\mathrm{vap}(T)$ versus $1/T$ and fitting the Antoine equation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4076545",
   "metadata": {},
   "source": [
    "Plot $\\ln P^\\mathrm{vap}$ versus $1/T$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "373ca4aa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30312 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Add 1/T to the dataframe:\n",
    "df_eq['1_over_T'] = 1 / df_eq[\"T\"]\n",
    "\n",
    "# Plot the result\n",
    "sns.lineplot(data=df_eq, x=\"1_over_T\", y=\"P\", color=\"b\", linewidth=2, label=\"PR EOS\")\n",
    "\n",
    "# Add the Antoine equation to the plot\n",
    "# Antoine equation given parameters A, B, and C\n",
    "# Pressure is returned in Pa\n",
    "def P_antoine(A, B, C, T):\n",
    "    return 10**(A - B/(T + C))*100000\n",
    "\n",
    "# Antoine equation parameters\n",
    "# Valid temperature range in K\n",
    "T_Antoine_min = 63.14\n",
    "T_Antoine_max = 126\n",
    "A = 3.7362\n",
    "B = 264.651\n",
    "C = -6.788\n",
    "\n",
    "# Plot the Antoine equation\n",
    "T = np.arange(T_Antoine_min, T_Antoine_max+1, 1)\n",
    "plt.plot(1/T, P_antoine(A, B, C, T), \"c\", label=\"Antoine eqn\")\n",
    "\n",
    "#plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(r\"$1/T$ (K$^{-1}$)\")\n",
    "plt.ylabel(r\"$P^\\mathrm{vap}$ (Pa)\")\n",
    "plt.title(r'$P^\\mathrm{vap}$ diagram for nitrogen, Peng-Robinson EOS')\n",
    "\n",
    "plt.legend(frameon=True, fancybox=True, edgecolor='none')\n",
    "#plt.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d34c490",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
