{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy import *\n",
      "init_session(quiet=True)\n",
      "init_printing()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "IPython console for SymPy 0.7.3 (Python 2.7.2-32-bit) (ground types: python)\n"
       ]
      }
     ],
     "prompt_number": 30
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Associating Fluids"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Yields a Helmholtz energy term in the form\n",
      "\n",
      "$$\\alpha^{ass} = A\\left(\\ln X-\\frac{X}{2}+\\frac{1}{2} \\right) $$\n",
      "where\n",
      "\n",
      "$$X = \\frac{2}{\\sqrt{1+4\\bar\\Delta\\delta}+1}$$\n",
      "\n",
      "$$\\bar\\Delta(\\delta,\\tau) = g(\\eta(\\delta))[\\exp(\\bar\\varepsilon\\tau)-1]\\bar\\kappa$$\n",
      "\n",
      "$$\\eta = \\bar v_n\\delta $$\n",
      "\n",
      "$$ g(\\eta) = \\frac{1}{2}\\frac{2-\\eta}{(1-\\eta)^3} $$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Derivatives of $g$ and $\\bar\\Delta$ with respect to $\\tau$ and $\\delta$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ \\frac{dg}{d\\eta} = 0.5\\frac{5-2\\eta}{(1-\\eta)^4} $$\n",
      "\n",
      "$$ \\frac{d^2g}{d\\eta^2} = 3\\frac{3-\\eta}{(1-\\eta)^5} $$\n",
      "\n",
      "$$ \\frac{d^3g}{d\\eta^3} = 6\\frac{7-2\\eta}{(1-\\eta)^6} $$\n",
      "\n",
      "$$\\left(\\frac{\\partial \\bar \\Delta}{\\partial \\delta} \\right)_{\\tau} = \\frac{dg(\\eta)}{d\\eta}[\\exp(\\bar\\varepsilon\\tau)-1]\\bar\\kappa\\bar v_n \\qquad \\mbox{(Eq. A35)}  $$ \n",
      "\n",
      "$$\\left(\\frac{\\partial \\bar \\Delta}{\\partial \\tau} \\right)_{\\delta} = g(\\eta)\\bar\\kappa\\exp(\\bar\\varepsilon\\tau)\\bar \\varepsilon \\qquad \\mbox{(Eq. A55)}$$\n",
      "\n",
      "$$\\left(\\frac{\\partial^2 \\bar \\Delta}{\\partial \\delta^2} \\right)_{\\tau} = \\frac{d^2g(\\eta)}{d\\eta^2}[\\exp(\\bar\\varepsilon\\tau)-1]\\bar\\kappa\\bar v_n^2 \\qquad \\mbox{(Eq. A76)}  $$ \n",
      "\n",
      "$$\\left(\\frac{\\partial^2 \\bar \\Delta}{\\partial \\tau^2} \\right)_{\\delta} = g(\\eta)\\bar\\kappa\\exp(\\bar\\varepsilon\\tau)\\bar \\varepsilon^2 \\qquad \\mbox{(Eq. A83)}$$\n",
      "\n",
      "$$\\frac{\\partial^2 \\bar \\Delta}{\\partial \\tau\\partial\\delta} = \\frac{dg(\\eta)}{d\\eta}\\exp(\\bar\\varepsilon\\tau)\\bar\\varepsilon\\bar\\kappa\\bar v_n \\qquad \\mbox{(Eq. A89)}$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Newly derived terms"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Following the patterns of the derivatives\n",
      "\n",
      "$$\\left(\\frac{\\partial^3 \\bar \\Delta}{\\partial \\tau^3} \\right)_{\\delta} = g(\\eta)\\bar\\kappa\\exp(\\bar\\varepsilon\\tau)\\bar \\varepsilon^3 $$\n",
      "\n",
      "$$\\frac{\\partial^3 \\bar \\Delta}{\\partial\\delta\\partial \\tau^2} = \\frac{dg(\\eta)}{d\\eta}\\bar\\kappa\\exp(\\bar\\varepsilon\\tau)\\bar \\varepsilon^2\\bar v_n$$\n",
      "\n",
      "$$\\frac{\\partial^3 \\bar \\Delta}{\\partial \\tau\\partial\\delta^2} = \\frac{d^2g(\\eta)}{d\\eta^2}\\exp(\\bar\\varepsilon\\tau)\\bar\\varepsilon\\bar\\kappa\\bar v_n^2$$\n",
      "\n",
      "$$\\left(\\frac{\\partial^3 \\bar \\Delta}{\\partial \\delta^3} \\right)_{\\tau} = \\frac{d^3g(\\eta)}{d\\eta^3}[\\exp(\\bar\\varepsilon\\tau)-1]\\bar\\kappa\\bar v_n^3 $$ "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Derivatives of $X$ with respect to $\\tau$ and $\\delta$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ \\left(\\frac{\\partial X}{\\partial \\bar\\Delta}\\right)_{\\delta} = -\\frac{\\delta X^2}{2\\bar\\Delta\\delta X+1} $$\n",
      "\n",
      "$$ \\left(\\frac{\\partial X}{\\partial \\delta}\\right)_{\\bar\\Delta} = -\\frac{\\bar\\Delta X^2}{2\\bar\\Delta\\delta X+1} $$\n",
      "\n",
      "$$X_{\\delta} = \\left(\\frac{\\partial X}{\\partial \\delta}\\right)_{\\bar\\Delta}+\\left(\\frac{\\partial X}{\\partial \\bar\\Delta}\\right)_{\\delta}\\left(\\frac{\\partial \\bar \\Delta}{\\partial\\delta}\\right)_{\\tau} \\qquad\\mbox{(Eq. A32)}$$\n",
      "\n",
      "$$X_{\\delta} = -\\frac{X^2}{2\\bar\\Delta\\delta X+1}\\left[\\bar\\Delta+\\delta\\left(\\frac{\\partial \\bar\\Delta}{\\partial\\delta}\\right)_{\\tau}\\right] \\qquad\\mbox{(Eq. A36)}$$\n",
      "\n",
      "$$X_{\\tau} = \\left(\\frac{\\partial X}{\\partial \\bar\\Delta}\\right)_{\\tau,\\delta}\\left(\\frac{\\partial \\bar \\Delta}{\\partial\\tau}\\right)_{\\delta} \\qquad\\mbox{(Eq. A54)}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$X_{\\delta\\delta} = \\left(\\frac{\\partial X_\\delta}{\\partial \\delta}\\right)_{X,\\bar\\Delta,\\bar\\Delta_\\delta}+\\frac{\\partial X_\\delta}{\\partial X}\\left(\\frac{\\partial X}{\\partial \\delta}\\right)_{\\bar\\Delta}\n",
      "+\\frac{\\partial X_\\delta}{\\partial X}\\frac{\\partial X}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\delta}+\\frac{\\partial X_\\delta}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\delta}+\\frac{\\partial X_\\delta}{\\partial \\bar\\Delta_\\delta}\\frac{\\partial \\bar\\Delta_\\delta}{\\partial \\delta}\\qquad \\mbox{(Eq. A74)}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$X_{\\tau\\tau} = \\frac{\\partial X_\\tau}{\\partial X}\\frac{\\partial X}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\tau}\n",
      "+\\frac{\\partial X_\\tau}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\tau}+\\frac{\\partial X_\\tau}{\\partial \\bar\\Delta_\\tau}\\frac{\\partial \\bar\\Delta_\\tau}{\\partial \\tau}\\qquad \\mbox{(Eq. A81)}$$\n",
      "\n",
      "$$X_{\\tau\\tau} = \\frac{2\\bar\\Delta_{\\tau}^2\\delta^2X^3(\\bar\\Delta\\delta X+1)}{(2\\bar\\Delta\\delta X+1)^3}+\\frac{2\\delta^2 X^3\\bar\\Delta_{\\tau}^2}{(2\\bar\\Delta\\delta X+1)^2}-\\frac{\\delta X^2\\bar\\Delta_{\\tau\\tau}}{1+2\\bar\\Delta\\delta X}$$\n",
      "\n",
      "$$X_{\\tau\\tau} = X_{\\tau\\tau}(X(\\bar\\Delta(\\tau)),\\bar\\Delta,\\bar\\Delta_\\tau,\\bar\\Delta_{\\tau\\tau})$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Third order partials"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Based on derivatives of $X_{\\tau\\tau}$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$X_{\\tau\\tau\\tau} = \\frac{\\partial X_{\\tau\\tau}}{\\partial X}\\frac{\\partial X}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\tau}\n",
      "+\\frac{\\partial X_{\\tau\\tau}}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\tau}+\\frac{\\partial X_{\\tau\\tau}}{\\partial \\bar\\Delta_\\tau}\\frac{\\partial \\bar\\Delta_\\tau}{\\partial \\tau}+\\frac{\\partial X_{\\tau\\tau}}{\\partial \\bar\\Delta_{\\tau\\tau}}\\frac{\\partial \\bar\\Delta_{\\tau\\tau}}{\\partial \\tau}$$\n",
      "\n",
      "$$X_{\\delta\\tau\\tau} = \\frac{\\partial X_{\\tau\\tau}}{\\partial\\delta} + \\frac{\\partial X_{\\tau\\tau}}{\\partial X}\\frac{\\partial X}{\\partial \\delta}+ \\frac{\\partial X_{\\tau\\tau}}{\\partial X}\\frac{\\partial X}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\delta}\n",
      "+\\frac{\\partial X_{\\tau\\tau}}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\delta}+\\frac{\\partial X_{\\tau\\tau}}{\\partial \\bar\\Delta_\\tau}\\frac{\\partial \\bar\\Delta_\\tau}{\\partial \\delta}+\\frac{\\partial X_{\\tau\\tau}}{\\partial \\bar\\Delta_{\\tau\\tau}}\\frac{\\partial \\bar\\Delta_{\\tau\\tau}}{\\partial \\delta}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X,delta,Delta,Delta_t,Delta_tt = symbols('X,delta,Delta,Delta_t,Delta_tt')\n",
      "X_tt = (2*Delta_t**2*delta**2*X**3*(Delta*delta*X+1)/(2*Delta*delta*X+1)**3\n",
      "        +(2*delta**2*X**3*Delta_t**2)/(2*Delta*delta*X+1)**2\n",
      "        -(delta*X**2*Delta_tt)/(2*Delta*delta*X+1)\n",
      "        )\n",
      "print 'dX = '+ccode(simplify(diff(X_tt,X)))\n",
      "print 'ddelta = '+ccode(simplify(diff(X_tt,delta)))\n",
      "print 'dDelta = '+ccode(simplify(diff(X_tt,Delta)))\n",
      "print 'dDelta_t = '+ccode(simplify(diff(X_tt,Delta_t)))\n",
      "print 'dDelta_tt = '+ccode(simplify(diff(X_tt,Delta_tt)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "dX = 2*X*delta*(-6*Delta*pow(Delta_t, 2)*pow(X, 2)*pow(delta, 2)*(Delta*X*delta + 1) + 3*pow(Delta_t, 2)*X*delta*(2*Delta*X*delta + 1) - Delta_tt*pow(2*Delta*X*delta + 1, 3) + X*delta*(Delta*Delta_tt + 3*pow(Delta_t, 2))*pow(2*Delta*X*delta + 1, 2))/pow(2*Delta*X*delta + 1, 4)\n",
        "ddelta = pow(X, 2)*(-12*Delta*pow(Delta_t, 2)*pow(X, 2)*pow(delta, 2)*(Delta*X*delta + 1) + 2*pow(Delta_t, 2)*X*delta*(-Delta*X*delta + 2)*(2*Delta*X*delta + 1) - Delta_tt*pow(2*Delta*X*delta + 1, 3) + 2*X*delta*(Delta*Delta_tt + 2*pow(Delta_t, 2))*pow(2*Delta*X*delta + 1, 2))/pow(2*Delta*X*delta + 1, 4)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "dDelta = 2*pow(X, 3)*pow(delta, 2)*(-6*pow(Delta_t, 2)*X*delta*(Delta*X*delta + 1) - 3*pow(Delta_t, 2)*X*delta*(2*Delta*X*delta + 1) + Delta_tt*pow(2*Delta*X*delta + 1, 2))/pow(2*Delta*X*delta + 1, 4)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "dDelta_t = 4*Delta_t*pow(X, 3)*pow(delta, 2)*(3*Delta*X*delta + 2)/pow(2*Delta*X*delta + 1, 3)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "dDelta_tt = -pow(X, 2)*delta/(2*Delta*X*delta + 1)\n"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Based on derivatives of $X_{\\delta\\delta}$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$X_{\\delta\\delta\\delta} = \\frac{X_{\\delta\\delta}}{\\partial \\delta} + \\frac{\\partial X_{\\delta\\delta}}{\\partial X}\\frac{\\partial X}{\\partial \\delta}+\\frac{\\partial X_{\\delta\\delta}}{\\partial X}\\frac{\\partial X}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\delta}\n",
      "+\\frac{\\partial X_{\\delta\\delta}}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\delta}+\\frac{\\partial X_{\\delta\\delta}}{\\partial \\bar\\Delta_\\delta}\\frac{\\partial \\bar\\Delta_\\delta}{\\partial \\delta}+\\frac{\\partial X_{\\delta\\delta}}{\\partial \\bar\\Delta_{\\delta\\delta}}\\frac{\\partial \\bar\\Delta_{\\delta\\delta}}{\\partial \\delta}$$\n",
      "\n",
      "$$X_{\\delta\\delta\\tau} = \\frac{\\partial X_{\\delta\\delta}}{\\partial X}\\frac{\\partial X}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\tau}\n",
      "+\\frac{\\partial X_{\\delta\\delta}}{\\partial \\bar\\Delta}\\frac{\\partial \\bar\\Delta}{\\partial \\tau}+\\frac{\\partial X_{\\delta\\delta}}{\\partial \\bar\\Delta_\\delta}\\frac{\\partial \\bar\\Delta_\\delta}{\\partial \\tau}+\\frac{\\partial X_{\\delta\\delta}}{\\partial \\bar\\Delta_{\\delta\\delta}}\\frac{\\partial \\bar\\Delta_{\\delta\\delta}}{\\partial \\tau}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X,delta,Delta,Delta_d,Delta_dd = symbols('X,delta,Delta,Delta_d,Delta_dd')\n",
      "            \n",
      "X_dd = (X**2*(2*Delta**2*X-Delta_d)/(2*Delta*delta*X+1)**2\n",
      "        -(Delta+delta*Delta_d)*2*(Delta*delta*X**2+X)/(2*Delta*delta*X+1)**2*(-Delta*X**2/(2*Delta*delta*X+1))\n",
      "        -(Delta+delta*Delta_d)*2*(Delta*delta*X**2+X)/(2*Delta*delta*X+1)**2*(-delta*X**2/(2*Delta*delta*X+1))*Delta_d\n",
      "        +X**2*(2*delta**2*X*Delta_d-1)/(2*Delta*delta*X+1)**2*Delta_d\n",
      "        -delta*X**2/(2*Delta*delta*X+1)*Delta_dd\n",
      "        )\n",
      "\n",
      "print 'dX = '+ccode(simplify(diff(X_dd,X)))\n",
      "print 'ddelta = '+ccode(simplify(diff(X_dd,delta)))\n",
      "print 'dDelta = '+ccode(simplify(diff(X_dd,Delta)))\n",
      "print 'dDelta_d = '+ccode(simplify(diff(X_dd,Delta_d)))\n",
      "print 'dDelta_dd = '+ccode(simplify(diff(X_dd,Delta_dd)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "dX = 2*X*(-6*Delta*pow(X, 2)*delta*pow(Delta + Delta_d*delta, 2)*(Delta*X*delta + 1) - Delta_dd*delta*pow(2*Delta*X*delta + 1, 3) + 2*X*(2*Delta*X*delta + 1)*(-Delta*Delta_d*delta*(2*Delta_d*X*pow(delta, 2) - 1) - Delta*delta*(2*pow(Delta, 2)*X - Delta_d) + Delta*(Delta + Delta_d*delta)*(Delta*X*delta + 1) + Delta_d*delta*(Delta + Delta_d*delta)*(Delta*X*delta + 1)) + pow(2*Delta*X*delta + 1, 2)*(3*pow(Delta, 2)*X + Delta*Delta_dd*X*pow(delta, 2) + Delta*X*(Delta + Delta_d*delta) + pow(Delta_d, 2)*X*pow(delta, 2) + Delta_d*X*delta*(Delta + Delta_d*delta) + Delta_d*(2*Delta_d*X*pow(delta, 2) - 1) - Delta_d))/pow(2*Delta*X*delta + 1, 4)\n",
        "ddelta = pow(X, 2)*(-24*pow(Delta, 4)*pow(X, 3)*delta - 8*pow(Delta, 3)*Delta_d*pow(X, 3)*pow(delta, 2) - 18*pow(Delta, 3)*pow(X, 2) + 8*pow(Delta, 2)*Delta_d*pow(X, 2)*delta - 4*pow(Delta, 2)*Delta_dd*pow(X, 2)*pow(delta, 2) + 10*Delta*pow(Delta_d, 2)*pow(X, 2)*pow(delta, 2) + 12*Delta*Delta_d*X - 4*Delta*Delta_dd*X*delta + 8*pow(Delta_d, 2)*X*delta - Delta_dd)/(16*pow(Delta, 4)*pow(X, 4)*pow(delta, 4) + 32*pow(Delta, 3)*pow(X, 3)*pow(delta, 3) + 24*pow(Delta, 2)*pow(X, 2)*pow(delta, 2) + 8*Delta*X*delta + 1)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "dDelta = pow(X, 3)*(-8*pow(Delta, 2)*Delta_d*pow(X, 2)*pow(delta, 3) + 8*pow(Delta, 2)*Delta_dd*pow(X, 2)*pow(delta, 4) + 10*pow(Delta, 2)*X*delta - 24*Delta*pow(Delta_d, 2)*pow(X, 2)*pow(delta, 4) + 8*Delta*Delta_d*X*pow(delta, 2) + 8*Delta*Delta_dd*X*pow(delta, 3) + 8*Delta - 18*pow(Delta_d, 2)*X*pow(delta, 3) + 12*Delta_d*delta + 2*Delta_dd*pow(delta, 2))/(16*pow(Delta, 4)*pow(X, 4)*pow(delta, 4) + 32*pow(Delta, 3)*pow(X, 3)*pow(delta, 3) + 24*pow(Delta, 2)*pow(X, 2)*pow(delta, 2) + 8*Delta*X*delta + 1)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "dDelta_d = 2*pow(X, 2)*(2*X*delta*(Delta + Delta_d*delta)*(Delta*X*delta + 1) + (2*Delta*X*delta + 1)*(2*Delta_d*X*pow(delta, 2) - 1))/pow(2*Delta*X*delta + 1, 3)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "dDelta_dd = -pow(X, 2)*delta/(2*Delta*X*delta + 1)\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Derivatives with respect to $\\tau$ and $\\delta$ of each term"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ \\frac{\\partial \\alpha^{ass}}{\\partial \\tau} = A\\left[\\frac{1}{X}-\\frac{1}{2}\\right]X_{\\tau} \\quad \\quad \\mathrm{(Eq. A52/A53)}$$ \n",
      "\n",
      "$$ \\frac{\\partial \\alpha^{ass}}{\\partial \\delta} = A\\left[\\frac{1}{X}-\\frac{1}{2}\\right]X_{\\delta}  \\quad \\quad \\mathrm{(Eq. A29/A30)}$$\n",
      "\n",
      "$$ \\frac{\\partial^2 \\alpha^{ass}}{\\partial \\delta^2} = A\\left[\\frac{1}{X}-\\frac{1}{2}\\right]X_{\\delta\\delta} -A\\frac{X_{\\delta}^2}{X^2} \\quad \\quad \\mathrm{(Eq. A72/A73)}$$\n",
      "\n",
      "$$ \\frac{\\partial^2 \\alpha^{ass}}{\\partial \\tau^2} = A\\left[\\frac{1}{X}-\\frac{1}{2}\\right]X_{\\tau\\tau} -A\\frac{X_{\\tau}^2}{X^2} \\quad \\quad \\mathrm{(Eq. A79/A80)}$$\n",
      "\n",
      "$$ \\frac{\\partial^2 \\alpha^{ass}}{\\partial \\tau\\partial\\delta} = A\\left[-\\frac{X_{\\tau}}{X^2}\\right]X_{\\delta}+AX_{\\delta\\tau}\\left[\\frac{1}{X}-\\frac{1}{2}\\right] \\quad \\quad \\mathrm{(Eq. A86/A87)}$$\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Newly Derived Terms"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ \\frac{\\partial^3 \\alpha^{ass}}{\\partial \\tau^3} =  A\\left[\\frac{1}{X}-\\frac{1}{2}\\right]X_{\\tau\\tau\\tau} + A\\left[-\\frac{X_{\\tau}}{X^2}\\right]X_{\\tau\\tau} -A\\frac{X^2(2X_{\\tau}X_{\\tau\\tau})-X_{\\tau}^2(2XX_{\\tau})}{X^4} $$\n",
      "\n",
      "$$ \\frac{\\partial^3 \\alpha^{ass}}{\\partial \\delta\\partial \\tau^2} =  A\\left[\\frac{1}{X}-\\frac{1}{2}\\right]X_{\\delta\\tau\\tau} + A\\left[-\\frac{X_{\\delta}}{X^2}\\right]X_{\\tau\\tau} -A\\frac{X^2(2X_{\\tau}X_{\\tau\\delta})-X_{\\tau}^2(2XX_{\\delta})}{X^4} $$\n",
      "\n",
      "$$ \\frac{\\partial^3 \\alpha^{ass}}{\\partial \\delta^2\\partial \\tau} =  A\\left[\\frac{1}{X}-\\frac{1}{2}\\right]X_{\\delta\\delta\\tau} + A\\left[-\\frac{X_{\\tau}}{X^2}\\right]X_{\\delta\\delta} -A\\frac{X^2(2X_{\\delta}X_{\\tau\\delta})-X_{\\delta}^2(2XX_{\\tau})}{X^4} $$\n",
      "\n",
      "$$ \\frac{\\partial^3 \\alpha^{ass}}{\\partial \\delta^3} =  A\\left[\\frac{1}{X}-\\frac{1}{2}\\right]X_{\\delta\\delta\\delta} + A\\left[-\\frac{X_{\\delta}}{X^2}\\right]X_{\\delta\\delta} -A\\frac{X^2(2X_{\\delta}X_{\\delta\\delta})-X_{\\delta}^2(2XX_{\\delta})}{X^4} $$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "from math import exp, log\n",
      "\n",
      "A = 2\n",
      "\n",
      "class Methanol(object):\n",
      "    \n",
      "    m = 0.977118832\n",
      "    vbarn = 0.204481952\n",
      "    kappabar = 0.148854832e-2\n",
      "    epsilonbar = 5.46341463\n",
      "    \n",
      "    def X(self, delta, Deltabar):\n",
      "        return 2/((1+4*Deltabar*delta)**0.5+1)\n",
      "    \n",
      "    def dX_dDeltabar__constdelta(self, delta, Deltabar):\n",
      "        X = self.X(delta,Deltabar)\n",
      "        return -delta*X**2/(2*Deltabar*delta*X+1)\n",
      "    \n",
      "    def dX_ddelta__constDeltabar(self, delta, Deltabar):\n",
      "        X = self.X(delta,Deltabar)\n",
      "        return -Deltabar*X**2/(2*Deltabar*delta*X+1)\n",
      "    \n",
      "    def dX_dtau(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        return self.dX_dDeltabar__constdelta(delta, Deltabar)*self.dDeltabar_dtau__constdelta(tau, delta)\n",
      "    \n",
      "    def dX_ddelta(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        return (self.dX_ddelta__constDeltabar(delta, Deltabar)\n",
      "               + self.dX_dDeltabar__constdelta(delta, Deltabar)*self.dDeltabar_ddelta__consttau(tau, delta))\n",
      "    \n",
      "    def d2X_dtau2(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        X = self.X(delta, Deltabar)\n",
      "        beta = self.dDeltabar_dtau__constdelta(tau, delta)\n",
      "        d_dXdtau_dbeta = -delta*X**2/(2*Deltabar*delta*X+1)\n",
      "        d_dXdtau_dDeltabar = 2*delta**2*X**3/(2*Deltabar*delta*X+1)**2*beta\n",
      "        d_dXdtau_dX = -2*beta*delta*X*(Deltabar*delta*X+1)/(2*Deltabar*delta*X+1)**2\n",
      "        dbeta_dtau = self.d2Deltabar_dtau2__constdelta(tau, delta)\n",
      "        dX_dDeltabar = self.dX_dDeltabar__constdelta(delta, Deltabar)\n",
      "        val = d_dXdtau_dX*dX_dDeltabar*beta+d_dXdtau_dDeltabar*beta+d_dXdtau_dbeta*dbeta_dtau\n",
      "        #print val, self.X_tt(X,delta,Deltabar,beta,dbeta_dtau)\n",
      "        return val\n",
      "    \n",
      "    def X_tt(self, X, delta, Delta, Delta_t, Delta_tt):\n",
      "        return (2*Delta_t**2*delta**2*X**3*(Delta*delta*X+1)/(2*Delta*delta*X+1)**3\n",
      "        +(2*delta**2*X**3*Delta_t**2)/(2*Delta*delta*X+1)**2\n",
      "        -(delta*X**2*Delta_tt)/(2*Delta*delta*X+1)\n",
      "        )\n",
      "                  \n",
      "    def d3X_dtau3(self, tau, delta):\n",
      "        Delta = self.Deltabar(tau, delta)\n",
      "        Delta_t = self.dDeltabar_dtau__constdelta(tau, delta)\n",
      "        Delta_tt = self.d2Deltabar_dtau2__constdelta(tau, delta)\n",
      "        dX_dDeltabar = self.dX_dDeltabar__constdelta(tau, delta)\n",
      "        dDeltabar_dtau = self.dDeltabar_dtau__constdelta(tau, delta)\n",
      "        dDeltabart_dtau = self.d2Deltabar_dtau2__constdelta(tau, delta)\n",
      "        dDeltabartt_dtau = self.d3Deltabar_dtau3__constdelta(tau, delta)\n",
      "        X = self.X(delta, Delta)\n",
      "        d = delta\n",
      "        t = tau\n",
      "        \n",
      "        dXtt_dDeltabartt = (-d*X**2)/(1+2*Delta*d*X)\n",
      "        dXtt_dDeltabart = ((4*Delta_t*d**2*X**3*(Delta*d*X+1))/(2*Delta*d*X+1)**3\n",
      "                           +(4*d**2*X**3*Delta_t)/(2*Delta*d*X+1)**2\n",
      "                           )\n",
      "        dXtt_dDeltabar = ((-12*Delta_t**2*d**3*X**4*(Delta*d*X+1))/(2*Delta*d*X+1)**4\n",
      "                        -(6*d**3*X**4*Delta_t**2)/(2*Delta*d*X+1)**3\n",
      "                        +(2*d**2*X**3*Delta_tt)/(2*Delta*d*X+1)**2\n",
      "                        )\n",
      "        dXtt_dX = (-(12*Delta_t**2*Delta*d**3*X**3*(Delta*d*X+1))/(2*Delta*d*X+1)**4\n",
      "                   +2*d**2*X**2*(3*Delta_t**2+Delta_tt*Delta)/(2*Delta*d*X+1)**2\n",
      "                   -(2*d*X*Delta_tt)/(2*Delta*d*X+1)\n",
      "                   )\n",
      "        \n",
      "        print '1', dXtt_dDeltabartt, (self.X_tt(X,delta,Delta,Delta_t,Delta_tt+1e-8)-self.X_tt(X,delta,Delta,Delta_t,Delta_tt-1e-8))/(2e-8)\n",
      "        print '2', dXtt_dDeltabart, (self.X_tt(X,delta,Delta,Delta_t+1e-8,Delta_tt)-self.X_tt(X,delta,Delta,Delta_t-1e-8,Delta_tt))/(2e-8)\n",
      "        print '3', dXtt_dDeltabar, (self.X_tt(X,delta,Delta+1e-8,Delta_t,Delta_tt)-self.X_tt(X,delta,Delta-1e-8,Delta_t,Delta_tt))/(2e-8)\n",
      "        print '4', dXtt_dX, (self.X_tt(X+1e-8,delta,Delta,Delta_t,Delta_tt)-self.X_tt(X-1e-8,delta,Delta,Delta_t,Delta_tt))/(2e-8)\n",
      "        \n",
      "        return (dXtt_dX*dX_dDeltabar*dDeltabar_dtau\n",
      "                +dXtt_dDeltabar*dDeltabar_dtau\n",
      "                +dXtt_dDeltabart*dDeltabart_dtau\n",
      "                +dXtt_dDeltabartt*dDeltabartt_dtau)\n",
      "    \n",
      "    def d2X_ddeltadtau(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        X = self.X(delta, Deltabar)\n",
      "        alpha = self.dDeltabar_ddelta__consttau(tau, delta)\n",
      "        beta = self.dDeltabar_dtau__constdelta(tau, delta)\n",
      "        dalpha_dtau = self.d2Deltabar_ddelta_dtau(tau, delta)\n",
      "        d_dXddelta_dDeltabar = X**2*(2*delta**2*X*alpha-1)/(2*Deltabar*delta*X+1)**2\n",
      "        d_dXddelta_dalpha = -delta*X**2/(2*Deltabar*delta*X+1)\n",
      "        d_dXddelta_dX = -(Deltabar+delta*alpha)*2*(Deltabar*delta*X**2+X)/(2*Deltabar*delta*X+1)**2\n",
      "        dX_dDeltabar = self.dX_dDeltabar__constdelta(delta, Deltabar)\n",
      "        return d_dXddelta_dX*dX_dDeltabar*beta+d_dXddelta_dDeltabar*beta+d_dXddelta_dalpha*dalpha_dtau\n",
      "    \n",
      "    def d2X_ddelta2(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        X = self.X(delta, Deltabar)\n",
      "        alpha = self.dDeltabar_ddelta__consttau(tau, delta)\n",
      "        dalpha_ddelta = self.d2Deltabar_ddelta2__consttau(tau, delta)\n",
      "        d_dXddelta_dDeltabar = X**2*(2*delta**2*X*alpha-1)/(2*Deltabar*delta*X+1)**2\n",
      "        d_dXddelta_dalpha = -delta*X**2/(2*Deltabar*delta*X+1)\n",
      "        d_dXddelta_dX = -(Deltabar+delta*alpha)*2*(Deltabar*delta*X**2+X)/(2*Deltabar*delta*X+1)**2\n",
      "        dX_dDeltabar = self.dX_dDeltabar__constdelta(delta, Deltabar)\n",
      "        dX_ddelta_constall = X**2*(2*Deltabar**2*X-alpha)/(2*Deltabar*delta*X+1)**2\n",
      "        return (dX_ddelta_constall\n",
      "                + d_dXddelta_dX*self.dX_ddelta__constDeltabar(delta, Deltabar)\n",
      "                + d_dXddelta_dX*dX_dDeltabar*alpha\n",
      "                + d_dXddelta_dDeltabar*alpha\n",
      "                + d_dXddelta_dalpha*dalpha_ddelta)\n",
      "        \n",
      "    def g(self, eta):\n",
      "        return 0.5*(2-eta)/(1-eta)**3\n",
      "    \n",
      "    def dg_deta(self, eta):\n",
      "        return 0.5*(5-2*eta)/(1-eta)**4\n",
      "    \n",
      "    def d2g_deta2(self, eta):\n",
      "        return 3*(3-eta)/(1-eta)**5\n",
      "    \n",
      "    def d3g_deta3(self, eta):\n",
      "        return 6*(7-2*eta)/(1-eta)**6\n",
      "    \n",
      "    def eta(self, delta):\n",
      "        return self.vbarn*delta\n",
      "    \n",
      "    def Deltabar(self, tau, delta):\n",
      "        return self.g(self.eta(delta))*(exp(self.epsilonbar*tau)-1)*self.kappabar\n",
      "    \n",
      "    def dDeltabar_ddelta__consttau(self, tau, delta):\n",
      "        return self.dg_deta(self.eta(delta))*(exp(self.epsilonbar*tau)-1)*self.kappabar*self.vbarn\n",
      "    \n",
      "    def d2Deltabar_ddelta2__consttau(self, tau, delta):\n",
      "        return self.d2g_deta2(self.eta(delta))*(exp(self.epsilonbar*tau)-1)*self.kappabar*self.vbarn**2\n",
      "    \n",
      "    def dDeltabar_dtau__constdelta(self, tau, delta):\n",
      "        return self.g(self.eta(delta))*self.kappabar*exp(self.epsilonbar*tau)*self.epsilonbar\n",
      "    \n",
      "    def d2Deltabar_dtau2__constdelta(self, tau, delta):\n",
      "        return self.g(self.eta(delta))*self.kappabar*exp(self.epsilonbar*tau)*self.epsilonbar**2\n",
      "    \n",
      "    def d3Deltabar_dtau3__constdelta(self, tau, delta):\n",
      "        return self.g(self.eta(delta))*self.kappabar*exp(self.epsilonbar*tau)*self.epsilonbar**3\n",
      "    \n",
      "    def d2Deltabar_ddelta_dtau(self, tau, delta):\n",
      "        return self.dg_deta(self.eta(delta))*exp(self.epsilonbar*tau)*self.epsilonbar*self.kappabar*self.vbarn\n",
      "    \n",
      "    def base(self, tau, delta):\n",
      "        X = self.X(delta, self.Deltabar(tau, delta))\n",
      "        return 2*(log(X)-X/2.0+0.5)\n",
      "    \n",
      "    def dDelta(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        X = self.X(delta, Deltabar)\n",
      "        return 2*(1/X-0.5)*self.dX_ddelta(tau, delta)\n",
      "    \n",
      "    def dTau(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        X = self.X(delta, Deltabar)\n",
      "        return 2*(1/X-0.5)*self.dX_dtau(tau, delta)\n",
      "    \n",
      "    def dTau2(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        X = self.X(delta, Deltabar)\n",
      "        X_tau = self.dX_dtau(tau, delta)\n",
      "        X_tautau = self.d2X_dtau2(tau, delta)\n",
      "        return 2*(1/X-0.5)*X_tautau-2*(X_tau/X)**2\n",
      "    \n",
      "    def dDeltadTau(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        X = self.X(delta, Deltabar)\n",
      "        X_delta = self.dX_ddelta(tau, delta)\n",
      "        X_tau = self.dX_dtau(tau, delta)\n",
      "        X_deltatau = self.d2X_ddeltadtau(tau, delta)\n",
      "        return 2*(-X_tau/X**2)*X_delta+2*X_deltatau*(1/X-0.5)\n",
      "    \n",
      "    def dDelta2(self, tau, delta):\n",
      "        Deltabar = self.Deltabar(tau, delta)\n",
      "        X = self.X(delta, Deltabar)\n",
      "        X_delta = self.dX_ddelta(tau, delta)\n",
      "        X_deltadelta = self.d2X_ddelta2(tau, delta)\n",
      "        return 2*(1/X-0.5)*X_deltadelta-2*(X_delta/X)**2\n",
      "    \n",
      "    def dTau3(self, tau, delta):\n",
      "        X = self.X(delta, Deltabar)\n",
      "        X_t = self.dX_dtau(tau, delta)\n",
      "        X_tt = self.d2X_dtau2(tau, delta)\n",
      "        X_ttt = self.d3X_dtau3(tau, delta)\n",
      "        return (A*(1/X-1/2)*X_ttt\n",
      "                +A*(-X_t/X**2)*X_tt\n",
      "                -2*A*(X**2*(X_t*X_tt)-X_t**2*(X*X_t))/X**4\n",
      "                )\n",
      "    \n",
      "M = Methanol()\n",
      "\n",
      "#print M.base(0.5,0.5)\n",
      "Deltabar = M.Deltabar(0.5,0.5)\n",
      "#print M.dX_dDeltabar__constdelta(0.5,Deltabar),(M.X(0.5,Deltabar+1e-8)-M.X(0.5,Deltabar))/1e-8\n",
      "\n",
      "print 'ttt', M.d3X_dtau3(0.5,0.5),(M.d2X_dtau2(0.5+1e-8,0.5)-M.d2X_dtau2(0.5-1e-8,0.5))/2e-8\n",
      "print 'd',M.dDelta(0.5,0.5), (M.base(0.5,0.5+1e-8)-M.base(0.5,0.5))/1e-8\n",
      "print 't',M.dTau(0.5,0.5), (M.base(0.5+1e-8,0.5)-M.base(0.5,0.5))/1e-8\n",
      "print 'tt',M.dTau2(0.5,0.5), (M.dTau(0.5+1e-8,0.5)-M.dTau(0.5,0.5))/1e-8\n",
      "print 'ttt',M.dTau3(0.5,0.5), (M.dTau2(0.5+1e-8,0.5)-M.dTau2(0.5,0.5))/1e-8\n",
      "print 'dt',M.dDeltadTau(0.5,0.5), (M.dTau(0.5,0.5+1e-8)-M.dTau(0.5,0.5))/1e-8\n",
      "print 'dd',M.dDelta2(0.5,0.5), (M.dDelta(0.5,0.5+1e-8)-M.dDelta(0.5,0.5))/1e-8"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "ttt 1 -0.473370012679 -0.473370015408\n",
        "2 0.295697260234 0.295697258346\n",
        "3 0.354430230606 0.354430232385\n",
        "4 -0.811327611566 -0.77524633102\n",
        "-1.95954772419 -1.93168185336\n",
        "d -0.0351186626009 -0.0351186635328\n",
        "t -0.0796836538477 -0.0796836374661\n",
        "tt -0.422816696531 -0.422816709422\n",
        "ttt 1 -0.473370012679 -0.473370015408\n",
        "2 0.295697260234 0.295697258346\n",
        "3 0.354430230606 0.354430232385\n",
        "4 -0.811327611566 -0.77524633102\n",
        "-2.20666413982 -2.17802790803\n",
        "dt -0.199708453268 -0.199708456716\n",
        "dd -0.0354391474199 -0.0354391481439\n"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}