{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helmholtz energy conversion of cubics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in Michelsen's book, a generalized volume-translated cubic EOS can be given the common form:\n",
    "$$ p = \\frac{RT}{v-b}-\\frac{a(T)}{(v+\\Delta_1b)(v+\\Delta_2b)} $$\n",
    "\n",
    "$$ p = \\frac{RT}{v+c-b}-\\frac{a(T)}{(v+c+\\Delta_1b)(v+c+\\Delta_2b)} $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPython console for SymPy 0.7.6.1 (Python 2.7.10-32-bit) (ground types: python)\n"
     ]
    }
   ],
   "source": [
    "from __future__ import division\n",
    "from sympy import *\n",
    "from IPython.display import display, Math, Latex\n",
    "from IPython.core.display import display_html\n",
    "init_session(quiet=True)\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rho_r,rho,tau,delta,T_r = symbols('rho_r,rho,tau,delta,T_r', positive = True, real = True)\n",
    "R,Delta_1,Delta_2,T,v,c = symbols('R,Delta_1,Delta_2,T,v,c', positive = True, real = True)\n",
    "x_i,x_j,x_k = symbols('x_i,x_j,x_k', positive=True, real = True)\n",
    "a = symbols('a', cls=Function, positive = True, real = True)(tau)\n",
    "b = symbols('b', positive = True, real = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "normal_cubic = True # True for no volume translation\n",
    "cubic = 'SRK'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$p = \\frac{R T}{- b + c + v} - \\frac{a{\\left (\\tau \\right )}}{\\left(c + v\\right) \\left(b + c + v\\right)}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 215,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = R*T/(v+c-b)-a/((v+c+Delta_1*b)*(v+c+Delta_2*b))\n",
    "if cubic == 'SRK':\n",
    "    p = p.subs(Delta_1,1).subs(Delta_2,0)\n",
    "elif cubic == 'VdW':\n",
    "    p = p.subs(Delta_1,0).subs(Delta_2,0)\n",
    "elif cubic == 'PR':\n",
    "    p = p.subs(Delta_1,1+sqrt(2)).subs(Delta_2,1-sqrt(2))\n",
    "Math('p = '+latex(p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$Z = \\frac{v}{R T} \\left(\\frac{R T}{- b + c + v} - \\frac{a{\\left (\\tau \\right )}}{\\left(c + v\\right) \\left(b + c + v\\right)}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 216,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = p*v/(R*T)\n",
    "Math('Z = '+latex(Z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$Z = \\frac{1}{R T \\rho} \\left(\\frac{R T}{- b + c + \\frac{1}{\\rho}} - \\frac{a{\\left (\\tau \\right )}}{\\left(c + \\frac{1}{\\rho}\\right) \\left(b + c + \\frac{1}{\\rho}\\right)}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$Z = \\frac{\\tau}{R T_{r} \\delta \\rho_{r}} \\left(\\frac{R T_{r}}{\\tau \\left(- b + c + \\frac{1}{\\delta \\rho_{r}}\\right)} - \\frac{a{\\left (\\tau \\right )}}{\\left(c + \\frac{1}{\\delta \\rho_{r}}\\right) \\left(b + c + \\frac{1}{\\delta \\rho_{r}}\\right)}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{d\\alpha^r}{d\\delta}|\\tau = - \\frac{\\tau a{\\left (\\tau \\right )}}{R T_{r} b c \\delta^{2} \\rho_{r} + R T_{r} b \\delta + R T_{r} c^{2} \\delta^{2} \\rho_{r} + 2 R T_{r} c \\delta + \\frac{R T_{r}}{\\rho_{r}}} + \\frac{\\tau}{- b \\delta^{2} \\rho_{r} \\tau + c \\delta^{2} \\rho_{r} \\tau + \\delta \\tau} - \\frac{1}{\\delta}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 217,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = Z.replace(v, 1/rho)\n",
    "display(Math('Z = '+latex(Z)))\n",
    "Z = Z.replace(T, T_r/tau)\n",
    "Z = Z.replace(rho, delta*rho_r)\n",
    "display(Math('Z = '+latex(Z)))\n",
    "dalphar_dDelta = (Z-1)/delta\n",
    "dalphar_dDelta = expand(simplify(dalphar_dDelta))\n",
    "Math('\\\\frac{d\\\\alpha^r}{d\\delta}|\\\\tau = ' + latex(dalphar_dDelta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\alpha^r = \\frac{1}{R T_{r} b} \\left(R T_{r} b \\log{\\left (- \\frac{1}{\\rho_{r} \\left(b - c\\right)} \\right )} - R T_{r} b \\log{\\left (\\frac{\\delta \\rho_{r} \\left(b - c\\right) - 1}{\\rho_{r} \\left(b - c\\right)} \\right )} + a{\\left (\\tau \\right )} \\log{\\left (\\frac{c^{\\tau} \\left(b + c \\delta \\rho_{r} \\left(b + c\\right) + c\\right)^{\\tau}}{\\left(c \\left(\\delta \\rho_{r} \\left(b + c\\right) + 1\\right)\\right)^{\\tau} \\left(b + c\\right)^{\\tau}} \\right )}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\lim_{c \\to 0}\\alpha^r = - \\log{\\left (b \\delta \\rho_{r} - 1 \\right )} + i \\pi - \\frac{\\tau a{\\left (\\tau \\right )}}{R T_{r} b} \\log{\\left (b \\delta \\rho_{r} + 1 \\right )}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "antideriv_pieces = 0\n",
    "for arg in expand(simplify(dalphar_dDelta))._args:\n",
    "    antideriv_pieces += integrate(arg,delta)\n",
    "\n",
    "antideriv_pieces = simplify(antideriv_pieces)\n",
    "alphar = antideriv_pieces - antideriv_pieces.subs(delta,0)\n",
    "\n",
    "display(Math('\\\\alpha^r = ' + latex(simplify(alphar))))\n",
    "display(Math('\\\\lim_{c \\\\to 0}\\\\alpha^r = ' + latex(simplify(limit(alphar,c,0)))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Maybe turn off volume translation\n",
    "if normal_cubic:\n",
    "    alphar = limit(alphar,c,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def format_deriv(arg, itau, idel, RHS):\n",
    "    \"\"\" \n",
    "    A function for giving a nice latex representation of \n",
    "    the partial derivative in question \n",
    "    \"\"\"\n",
    "    if itau+idel == 1:\n",
    "        numexp = ''\n",
    "    else:\n",
    "        numexp = '^{{{s:d}}}'.format(s=itau+idel)\n",
    "        \n",
    "    if itau == 0:\n",
    "        tau = ''\n",
    "    elif itau == 1:\n",
    "        tau = '\\\\partial \\\\tau'\n",
    "    else:\n",
    "        tau = '\\\\partial \\\\tau^{{{s:d}}}'.format(s=itau)\n",
    "        \n",
    "    if idel == 0:\n",
    "        delta = ''\n",
    "    elif idel == 1:\n",
    "        delta = '\\\\partial \\\\delta'\n",
    "    else:\n",
    "        delta = '\\\\partial \\\\delta^{{{s:d}}}'.format(s=idel)\n",
    "        \n",
    "    temp = '\\\\frac{{\\\\partial{{{numexp:s}}} {arg:s}}}{{{{{tau:s}}}{{{delta:s}}}}} = '\n",
    "    if itau == 0 and idel == 0:\n",
    "        as_latex = arg +' = ' + latex(RHS)\n",
    "    else:\n",
    "        as_latex = temp.format(**locals()) + latex(RHS)\n",
    "    return Math(as_latex), as_latex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{} \\alpha^r}{{}{\\partial \\delta}} = \\frac{1}{R T_{r} b} \\left(- \\frac{R T_{r} b}{\\delta - \\frac{1}{b \\rho_{r}}} - \\frac{b \\rho_{r} \\tau a{\\left (\\tau \\right )}}{b \\delta \\rho_{r} + 1}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{} \\alpha^r}{{\\partial \\tau}{}} = \\frac{1}{R T_{r} b} \\left(- \\tau \\log{\\left (b \\delta \\rho_{r} + 1 \\right )} \\frac{d}{d \\tau} a{\\left (\\tau \\right )} - a{\\left (\\tau \\right )} \\log{\\left (b \\delta \\rho_{r} + 1 \\right )}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{2}} \\alpha^r}{{}{\\partial \\delta^{2}}} = \\frac{1}{R T_{r}} \\left(\\frac{R T_{r}}{\\left(\\delta - \\frac{1}{b \\rho_{r}}\\right)^{2}} + \\frac{b \\rho_{r}^{2} \\tau a{\\left (\\tau \\right )}}{\\left(b \\delta \\rho_{r} + 1\\right)^{2}}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{2}} \\alpha^r}{{\\partial \\tau}{\\partial \\delta}} = \\frac{1}{R T_{r} b} \\left(- \\frac{b \\rho_{r} \\tau \\frac{d}{d \\tau} a{\\left (\\tau \\right )}}{b \\delta \\rho_{r} + 1} - \\frac{b \\rho_{r} a{\\left (\\tau \\right )}}{b \\delta \\rho_{r} + 1}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{2}} \\alpha^r}{{\\partial \\tau^{2}}{}} = - \\frac{1}{R T_{r} b} \\left(\\tau \\frac{d^{2}}{d \\tau^{2}}  a{\\left (\\tau \\right )} + 2 \\frac{d}{d \\tau} a{\\left (\\tau \\right )}\\right) \\log{\\left (b \\delta \\rho_{r} + 1 \\right )}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{3}} \\alpha^r}{{}{\\partial \\delta^{3}}} = - \\frac{1}{R T_{r}} \\left(\\frac{2 R T_{r}}{\\left(\\delta - \\frac{1}{b \\rho_{r}}\\right)^{3}} + \\frac{2 b^{2} \\rho_{r}^{3} \\tau a{\\left (\\tau \\right )}}{\\left(b \\delta \\rho_{r} + 1\\right)^{3}}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{3}} \\alpha^r}{{\\partial \\tau}{\\partial \\delta^{2}}} = \\frac{b \\rho_{r}^{2} \\left(\\tau \\frac{d}{d \\tau} a{\\left (\\tau \\right )} + a{\\left (\\tau \\right )}\\right)}{R T_{r} \\left(b \\delta \\rho_{r} + 1\\right)^{2}}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{3}} \\alpha^r}{{\\partial \\tau^{2}}{\\partial \\delta}} = - \\frac{\\rho_{r}}{R T_{r} \\left(b \\delta \\rho_{r} + 1\\right)} \\left(\\tau \\frac{d^{2}}{d \\tau^{2}}  a{\\left (\\tau \\right )} + 2 \\frac{d}{d \\tau} a{\\left (\\tau \\right )}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{3}} \\alpha^r}{{\\partial \\tau^{3}}{}} = - \\frac{1}{R T_{r} b} \\left(\\tau \\frac{d^{3}}{d \\tau^{3}}  a{\\left (\\tau \\right )} + 3 \\frac{d^{2}}{d \\tau^{2}}  a{\\left (\\tau \\right )}\\right) \\log{\\left (b \\delta \\rho_{r} + 1 \\right )}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{4}} \\alpha^r}{{}{\\partial \\delta^{4}}} = \\frac{1}{R T_{r}} \\left(\\frac{6 R T_{r}}{\\left(\\delta - \\frac{1}{b \\rho_{r}}\\right)^{4}} + \\frac{6 b^{3} \\rho_{r}^{4} \\tau a{\\left (\\tau \\right )}}{\\left(b \\delta \\rho_{r} + 1\\right)^{4}}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{4}} \\alpha^r}{{\\partial \\tau}{\\partial \\delta^{3}}} = - \\frac{2 b^{2} \\rho_{r}^{3}}{R T_{r} \\left(b \\delta \\rho_{r} + 1\\right)^{3}} \\left(\\tau \\frac{d}{d \\tau} a{\\left (\\tau \\right )} + a{\\left (\\tau \\right )}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{4}} \\alpha^r}{{\\partial \\tau^{2}}{\\partial \\delta^{2}}} = \\frac{b \\rho_{r}^{2}}{R T_{r} \\left(b \\delta \\rho_{r} + 1\\right)^{2}} \\left(\\tau \\frac{d^{2}}{d \\tau^{2}}  a{\\left (\\tau \\right )} + 2 \\frac{d}{d \\tau} a{\\left (\\tau \\right )}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{4}} \\alpha^r}{{\\partial \\tau^{3}}{\\partial \\delta}} = - \\frac{\\rho_{r}}{R T_{r} \\left(b \\delta \\rho_{r} + 1\\right)} \\left(\\tau \\frac{d^{3}}{d \\tau^{3}}  a{\\left (\\tau \\right )} + 3 \\frac{d^{2}}{d \\tau^{2}}  a{\\left (\\tau \\right )}\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{4}} \\alpha^r}{{\\partial \\tau^{4}}{}} = - \\frac{1}{R T_{r} b} \\left(\\tau \\frac{d^{4}}{d \\tau^{4}}  a{\\left (\\tau \\right )} + 4 \\frac{d^{3}}{d \\tau^{3}}  a{\\left (\\tau \\right )}\\right) \\log{\\left (b \\delta \\rho_{r} + 1 \\right )}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for deriv_count in range(1,5):\n",
    "    for dtau in range(deriv_count+1):\n",
    "        ddelta = deriv_count-dtau\n",
    "        #print dtau, ddelta\n",
    "        mth, tex = format_deriv('\\\\alpha^r', dtau, ddelta, \n",
    "                             diff(diff(alphar,tau,dtau),delta,ddelta))\n",
    "        #print tex\n",
    "        display(mth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Derivatives of $a_i(\\tau)$\n",
    "=============="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$a_i = a_{0} \\left(\\kappa \\left(1 - \\frac{1}{\\sqrt{\\tau}}\\right) + 1\\right)^{2}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{} a_i}{{\\partial \\tau}{}} = \\frac{a_{0} \\kappa}{\\tau^{\\frac{3}{2}}} \\left(\\kappa \\left(1 - \\frac{1}{\\sqrt{\\tau}}\\right) + 1\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{2}} a_i}{{\\partial \\tau^{2}}{}} = \\frac{a_{0} \\kappa}{2} \\left(\\frac{\\kappa}{\\tau^{3}} - \\frac{1}{\\tau^{\\frac{5}{2}}} \\left(3 \\kappa \\left(1 - \\frac{1}{\\sqrt{\\tau}}\\right) + 3\\right)\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{3}} a_i}{{\\partial \\tau^{3}}{}} = \\frac{3 a_{0}}{4} \\kappa \\left(- \\frac{3 \\kappa}{\\tau^{4}} + \\frac{1}{\\tau^{\\frac{7}{2}}} \\left(5 \\kappa \\left(1 - \\frac{1}{\\sqrt{\\tau}}\\right) + 5\\right)\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial{^{4}} a_i}{{\\partial \\tau^{4}}{}} = \\frac{3 a_{0}}{8} \\kappa \\left(\\frac{29 \\kappa}{\\tau^{5}} - \\frac{1}{\\tau^{\\frac{9}{2}}} \\left(35 \\kappa \\left(1 - \\frac{1}{\\sqrt{\\tau}}\\right) + 35\\right)\\right)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a0,omega,kappa = symbols('a0,omega,kappa')\n",
    "ai=a0*(1+kappa*(1-sqrt(1/tau)))**2\n",
    "                \n",
    "for diff_count in range(0,5):\n",
    "    mth, tex = format_deriv('a_i',diff_count, 0, diff(ai, tau, diff_count))\n",
    "    display(mth)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
