{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.3.0\n",
      "ba41b1434002b1fce6e596386861bca2a890fb6f\n"
     ]
    }
   ],
   "source": [
    "############################################################\n",
    "# Plots of meta-stable Maxwell loops\n",
    "# Inspired by https://doi.org/10.1134/S0036024406040030\n",
    "# Math tricks taken from: http://math.stackexchange.com/q/416823/92706\n",
    "# Plot also shown on page 79 of https://doi.org/10.15480/882.1207\n",
    "############################################################\n",
    "\n",
    "# load some bits and pieces\n",
    "import numpy as np\n",
    "from numpy.linalg import solve\n",
    "from numpy.linalg import lstsq\n",
    "from numpy import log\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import CoolProp as CP\n",
    "from CoolProp.CoolProp import PropsSI\n",
    "\n",
    "# Check: CoolProp version\n",
    "print(CP.__version__)\n",
    "print(CP.__gitrevision__)\n",
    "\n",
    "# Constants\n",
    "eps = 1e-3\n",
    "kilo = 1e3\n",
    "Mega = 1e6\n",
    "golden = (1 + 5 ** 0.5) / 2\n",
    "width = 12.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R = 8.31451\n",
      "MM = 0.0440098\n",
      "Rs = 188.92405782348476\n",
      "T_crt = 304.1282\n",
      "T_trp = 216.592\n"
     ]
    }
   ],
   "source": [
    "# Set FluidName\n",
    "FluidName = 'CO2'\n",
    "nPoints = 1000\n",
    "# pick any int smaller than nPoints\n",
    "myIdx = 860 \n",
    "\n",
    "# Constants, triple and critical data\n",
    "R = PropsSI('GAS_CONSTANT',FluidName)\n",
    "MM = PropsSI('MOLAR_MASS',FluidName)\n",
    "Rs = R/MM\n",
    "T_crt = PropsSI('T_CRITICAL',FluidName)\n",
    "T_trp = PropsSI('T_TRIPLE',FluidName)\n",
    "p_crt = PropsSI('P_CRITICAL',FluidName)\n",
    "p_trp = PropsSI('P_TRIPLE',FluidName)\n",
    "p_max = PropsSI('P_MAX',FluidName)\n",
    "d_crt = PropsSI('RHOMASS_CRITICAL',FluidName)\n",
    "v_crt = 1/d_crt\n",
    "d_trp_liq = PropsSI('D','T',T_trp,'Q',0,FluidName)\n",
    "d_trp_vap = PropsSI('D','T',T_trp,'Q',1,FluidName)\n",
    "print(\"R = \" + str(R))\n",
    "print(\"MM = \" + str(MM))\n",
    "print(\"Rs = \" + str(Rs))\n",
    "print(\"T_crt = \" + str(T_crt))\n",
    "print(\"T_trp = \" + str(T_trp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculation of the coefficients for the metastable region interpolation happens in this cell\n",
    "\n",
    "T_sat = np.linspace(T_trp, T_crt-eps, num=nPoints)\n",
    "# empty arrays\n",
    "# vap side\n",
    "delta_vap = np.empty(nPoints)\n",
    "tau_vap = np.empty(nPoints)\n",
    "p_vap = np.empty(nPoints)\n",
    "d_vap = np.empty(nPoints)\n",
    "v_vap = np.empty(nPoints)\n",
    "f_vap = np.empty(nPoints)\n",
    "dP_dD_T_vap = np.empty(nPoints)\n",
    "d2P_dD2_T_vap = np.empty(nPoints)\n",
    "d2P_dDdT_vap = np.empty(nPoints)\n",
    "# liq side\n",
    "delta_liq = np.empty(nPoints)\n",
    "tau_liq = np.empty(nPoints)\n",
    "p_liq = np.empty(nPoints)\n",
    "d_liq = np.empty(nPoints)\n",
    "v_liq = np.empty(nPoints)\n",
    "f_liq = np.empty(nPoints)\n",
    "dP_dD_T_liq = np.empty(nPoints)\n",
    "d2P_dD2_T_liq = np.empty(nPoints)\n",
    "d2P_dDdT_liq = np.empty(nPoints)\n",
    "# metastable coeffs: \n",
    "AShape = (8,8)\n",
    "A = np.empty(AShape)\n",
    "b = np.empty(8)\n",
    "xShape = (nPoints,8)\n",
    "x = np.empty(xShape)\n",
    "\n",
    "HEOS = CP.AbstractState(\"HEOS\", FluidName)\n",
    "# get values from CoolProp\n",
    "for idx in range(0,nPoints):\n",
    "    # AT the vap line\n",
    "    HEOS.update(CP.QT_INPUTS, 1, T_sat[idx])  \n",
    "    delta_vap[idx] = HEOS.delta() \n",
    "    tau_vap[idx] = HEOS.tau()\n",
    "    p_vap[idx] = HEOS.p()\n",
    "    d_vap[idx] = HEOS.rhomass()\n",
    "    f_vap[idx] = Rs*T_sat[idx]*( HEOS.alpha0() + HEOS.alphar() )    \n",
    "    #f_vap[idx] = HEOS.umass() - T_sat[idx]*HEOS.smass()\n",
    "    dP_dD_T_vap[idx] = HEOS.first_partial_deriv(CP.iP, CP.iDmass, CP.iT)\n",
    "    d2P_dD2_T_vap[idx] = HEOS.second_partial_deriv(CP.iP, CP.iDmass, CP.iT, CP.iDmass, CP.iT)\n",
    "    d2P_dDdT_vap[idx] = HEOS.second_partial_deriv(CP.iP, CP.iDmass, CP.iT, CP.iT, CP.iDmass)     \n",
    "    \n",
    "    # AT the liq line\n",
    "    HEOS.update(CP.QT_INPUTS, 0, T_sat[idx])  \n",
    "    delta_liq[idx] = HEOS.delta() \n",
    "    tau_liq[idx] = HEOS.tau()\n",
    "    p_liq[idx] = HEOS.p() \n",
    "    d_liq[idx] = HEOS.rhomass() \n",
    "    f_liq[idx] = Rs*T_sat[idx]*( HEOS.alpha0() + HEOS.alphar() )\n",
    "    # f_liq[idx] = HEOS.umass() - T_sat[idx]*HEOS.smass()\n",
    "    dP_dD_T_liq[idx] = HEOS.first_partial_deriv(CP.iP, CP.iDmass, CP.iT)\n",
    "    d2P_dD2_T_liq[idx] = HEOS.second_partial_deriv(CP.iP, CP.iDmass, CP.iT, CP.iDmass, CP.iT)   \n",
    "    d2P_dDdT_liq[idx] = HEOS.second_partial_deriv(CP.iP, CP.iDmass, CP.iT, CP.iT, CP.iDmass)   \n",
    "\n",
    "    # calculate metastable coeffs by solving Ax=b\n",
    "    A = np.array([  [1/tau_vap[idx], -1/delta_vap[idx]/tau_vap[idx],  log(delta_vap[idx]),          delta_vap[idx],     delta_vap[idx]**2/2,       delta_vap[idx]**3/3,         delta_vap[idx]**4/4,        delta_vap[idx]**5/5 ], \n",
    "                    [1/tau_liq[idx], -1/delta_liq[idx]/tau_liq[idx],  log(delta_liq[idx]),          delta_liq[idx],     delta_liq[idx]**2/2,       delta_liq[idx]**3/3,         delta_liq[idx]**4/4,        delta_liq[idx]**5/5 ], \n",
    "                    [             0,             d_crt/tau_vap[idx], d_crt*delta_vap[idx], d_crt*delta_vap[idx]**2, d_crt*delta_vap[idx]**3,    d_crt*delta_vap[idx]**4,    d_crt*delta_vap[idx]**5,    d_crt*delta_vap[idx]**6 ], \n",
    "                    [             0,             d_crt/tau_liq[idx], d_crt*delta_liq[idx], d_crt*delta_liq[idx]**2, d_crt*delta_liq[idx]**3,    d_crt*delta_liq[idx]**4,    d_crt*delta_liq[idx]**5,    d_crt*delta_liq[idx]**6 ], \n",
    "                    [             0,                              0,                    1,        2*delta_vap[idx],     3*delta_vap[idx]**2,        4*delta_vap[idx]**3,        5*delta_vap[idx]**4,        6*delta_vap[idx]**5 ], \n",
    "                    [             0,                              0,                    1,        2*delta_liq[idx],     3*delta_liq[idx]**2,        4*delta_liq[idx]**3,        5*delta_liq[idx]**4,        6*delta_liq[idx]**5 ], \n",
    "                    [             0,                              0,                    0,                 2/d_crt,  6*delta_vap[idx]/d_crt, 12*delta_vap[idx]**2/d_crt, 20*delta_vap[idx]**3/d_crt, 30*delta_vap[idx]**4/d_crt ], \n",
    "                    [             0,                              0,                    0,                 2/d_crt,  6*delta_liq[idx]/d_crt, 12*delta_liq[idx]**2/d_crt, 20*delta_liq[idx]**3/d_crt, 30*delta_liq[idx]**4/d_crt ]])\n",
    "    A = Rs*T_crt*A\n",
    "    b = np.array([f_vap[idx], f_liq[idx], p_vap[idx], p_liq[idx], dP_dD_T_vap[idx], dP_dD_T_liq[idx], d2P_dD2_T_vap[idx], d2P_dD2_T_liq[idx]])\n",
    "    x[idx] = solve(A,b)\n",
    "    \n",
    "    # for validation\n",
    "    if (abs(idx-myIdx)<0.9):\n",
    "        print(np.allclose(np.dot(A, x[idx]), b))\n",
    "        print(A)\n",
    "        print(b)\n",
    "        print(x[idx])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just some validation plots\n",
    "plt.figure(figsize=(width,width*4/2/golden))\n",
    "\n",
    "plt.subplot(4,2,1)\n",
    "plt.plot(T_sat, f_vap, color='red')\n",
    "plt.plot(T_sat, f_liq, color='blue')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.ylabel('Helmholtz energy')\n",
    "\n",
    "plt.subplot(4,2,2)\n",
    "plt.plot(d_vap, T_sat, color='red')\n",
    "plt.plot(d_liq, T_sat, color='blue')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.xlabel('Density in kg/m³')\n",
    "plt.ylabel('Temperature in K')\n",
    "\n",
    "plt.subplot(4,2,3)\n",
    "plt.plot(T_sat, p_vap, color='red')\n",
    "plt.plot(T_sat, p_liq, color='blue')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.ylabel('Pressure in Pa')\n",
    "\n",
    "plt.subplot(4,2,4)\n",
    "plt.plot(d_vap, p_vap, color='red')\n",
    "plt.plot(d_liq, p_liq, color='blue')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.xlabel('Density in kg/m³')\n",
    "plt.ylabel('Pressure in Pa')\n",
    "\n",
    "plt.subplot(4,2,5)\n",
    "plt.plot(T_sat, dP_dD_T_vap, color='red')\n",
    "plt.plot(T_sat, dP_dD_T_liq, color='blue')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.yscale('log')\n",
    "plt.ylabel('dP_dD_T')\n",
    "\n",
    "plt.subplot(4,2,6)\n",
    "plt.plot(d_vap, dP_dD_T_vap, color='red')\n",
    "plt.plot(d_liq, dP_dD_T_liq, color='blue')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Density in kg/m³')\n",
    "plt.ylabel('dP_dD_T')\n",
    "\n",
    "plt.subplot(4,2,7)\n",
    "plt.plot(T_sat, d2P_dD2_T_vap, color='red')\n",
    "plt.plot(T_sat, d2P_dD2_T_liq, color='blue')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.ylabel('d2P_dD2_T')\n",
    "plt.xlabel('Temperature in K')\n",
    "\n",
    "plt.subplot(4,2,8)\n",
    "plt.plot(d_vap, d2P_dD2_T_vap, color='red')\n",
    "plt.plot(d_liq, d2P_dD2_T_liq, color='blue')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.xlabel('Density in kg/m³')\n",
    "plt.ylabel('d2P_dD2_T')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tau_iso = 1.0417217720566936\n",
      "T_iso = 291.94762762762764\n",
      "p_sat(T_iso) = 5569191.928520719\n",
      "d_vap(T_iso) = 185.17587066325405\n",
      "d_liq(T_iso) = 785.8632178765519\n",
      "coeffs = [-0.38011836 -0.0331094   1.21278872 -2.00233849  1.1397745   0.285719\n",
      " -0.5838537   0.17728102]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$v/v_c$')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 80152 characters]\n",
      "text/plain": [
       "<Figure size 900x1112.46 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This cell is for plotting the Maxwell loops along one specified isotherm (chosen with myIdx)\n",
    "T_iso = T_sat[myIdx]\n",
    "tau_iso = T_crt/T_iso\n",
    "p_iso = p_vap[myIdx]\n",
    "print(\"tau_iso = \" + str(tau_iso))\n",
    "print(\"T_iso = \" + str(T_iso))\n",
    "print(\"p_sat(T_iso) = \" + str(p_iso))\n",
    "print(\"d_vap(T_iso) = \" + str(d_vap[myIdx]))\n",
    "print(\"d_liq(T_iso) = \" + str(d_liq[myIdx]))\n",
    "\n",
    "# copy coeffs of that single point\n",
    "c = x[myIdx,:]\n",
    "print(\"coeffs = \" + str(c))\n",
    "      \n",
    "# get a density range\n",
    "d_min = 0.8*d_vap[myIdx] \n",
    "#d_min = d_trp_vap\n",
    "# d_max = 1.9*d_liq[myIdx]\n",
    "d_max = PropsSI('D','T',T_iso,'P',p_max,FluidName)\n",
    "d_max = d_trp_liq\n",
    "rhos = np.linspace(d_min, d_max, num=nPoints)\n",
    "deltas = rhos/d_crt\n",
    "# for plotting, we will use volume (d_min is high v, d_max is low v)\n",
    "vs = 1/rhos\n",
    "\n",
    "# calculate Helmholtz energy and pressure for that density range, at T_iso\n",
    "# stable\n",
    "fs = np.ones(nPoints)\n",
    "ps = np.ones(nPoints)\n",
    "# equation of state, without chacking phase\n",
    "fes = np.ones(nPoints)\n",
    "pes = np.ones(nPoints)\n",
    "# meta-stable\n",
    "fms = np.ones(nPoints)\n",
    "pms = np.ones(nPoints)\n",
    "for idx in range(0,nPoints):\n",
    "    # stable\n",
    "    HEOS.unspecify_phase()\n",
    "    HEOS.update(CP.DmassT_INPUTS, rhos[idx], T_iso)  \n",
    "    #fs[idx] = Rs*T_iso*(HEOS.alpha0() + HEOS.alphar())\n",
    "    fs[idx] = HEOS.umass() - T_iso*HEOS.smass()\n",
    "    ps[idx] = HEOS.p()\n",
    "    # eos\n",
    "    HEOS.specify_phase(CP.iphase_liquid)\n",
    "    HEOS.update(CP.DmassT_INPUTS, rhos[idx], T_iso)\n",
    "    fes[idx] = HEOS.umass() - T_iso*HEOS.smass()\n",
    "    pes[idx] = HEOS.p()\n",
    "    # meta stable interpolation\n",
    "    # p = -(df/dv)_T\n",
    "    # s = -(df/dT)_v\n",
    "    # cv = T*(ds/dT)_v = -T*(dsf/dT2)_v\n",
    "    fms[idx]  = Rs*T_crt*( +c[0]/tau_iso -c[1]/tau_iso/deltas[idx] +c[2]*log(deltas[idx]) +c[3]*deltas[idx]    +c[4]*deltas[idx]**2/2 +c[5]*deltas[idx]**3/3 +c[6]*deltas[idx]**4/4 +c[7]*deltas[idx]**5/5 )\n",
    "    pms[idx]  = Rs*T_crt*d_crt*( 0       +c[1]/tau_iso             +c[2]*deltas[idx]      +c[3]*deltas[idx]**2 +c[4]*deltas[idx]**3   +c[5]*deltas[idx]**4   +c[6]*deltas[idx]**5   +c[7]*deltas[idx]**6   )\n",
    "\n",
    "# now plot \n",
    "plt.figure(figsize=(width,width*2/1/golden))\n",
    "\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(vs/v_crt, fs/Rs/T_iso, color='red', label='stable equilibrium')\n",
    "plt.plot(vs/v_crt, fes/Rs/T_iso, color='green', label='EoS')\n",
    "plt.plot(vs/v_crt, fms/Rs/T_iso, color='blue', label='metastable Maxwell loop')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.legend(loc='upper right')\n",
    "plt.ylabel(r'$\\alpha$')\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(vs/v_crt, ps/p_crt, color='red', label='stable equilibrium')\n",
    "plt.plot(vs/v_crt, pes/p_crt, color='green', label='EoS')\n",
    "plt.plot(vs/v_crt, pms/p_crt, color='blue', label='metastable Maxwell loop')\n",
    "plt.plot(d_crt/d_liq[myIdx], p_iso/p_crt, 'b|', markersize=20)\n",
    "plt.plot(d_crt/d_vap[myIdx], p_iso/p_crt, 'r|', markersize=20)\n",
    "plt.ylim(top=7)\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.legend(loc='upper right')\n",
    "plt.ylabel(r'$p/p_c$')\n",
    "plt.xlabel(r'$v/v_c$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copy coeffs for all points\n",
    "c0 = x[:,0]\n",
    "c1 = x[:,1]\n",
    "c2 = x[:,2]\n",
    "c3 = x[:,3]\n",
    "c4 = x[:,4]\n",
    "c5 = x[:,5]\n",
    "c6 = x[:,6]\n",
    "c7 = x[:,7]\n",
    "\n",
    "# plot the values of all 8 coeffs over T_sat\n",
    "plt.figure(figsize=(width,width*4/2/golden))\n",
    "\n",
    "plt.subplot(4,2,1)\n",
    "plt.plot(T_sat, c0, label='c0')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.legend(loc='upper center')\n",
    "\n",
    "plt.subplot(4,2,2)\n",
    "plt.plot(T_sat, c1, label='c1')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.legend(loc='upper center')\n",
    "\n",
    "plt.subplot(4,2,3)\n",
    "plt.plot(T_sat, c2, label='c2')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.legend(loc='upper center')\n",
    "\n",
    "plt.subplot(4,2,4)\n",
    "plt.plot(T_sat, c3, label='c3')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.legend(loc='upper center')\n",
    "\n",
    "plt.subplot(4,2,5)\n",
    "plt.plot(T_sat, c4, label='c4')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.legend(loc='upper center')\n",
    "\n",
    "plt.subplot(4,2,6)\n",
    "plt.plot(T_sat, c5, label='c5')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.legend(loc='upper center')\n",
    "\n",
    "plt.subplot(4,2,7)\n",
    "plt.plot(T_sat, c6, label='c6')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.legend(loc='upper center')\n",
    "\n",
    "plt.subplot(4,2,8)\n",
    "plt.plot(T_sat, c7, label='c7')\n",
    "plt.grid(b=True, linestyle=':')\n",
    "plt.minorticks_on()\n",
    "plt.legend(loc='upper center')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
