

In [30]:

    
from sympy import *
init_session(quiet=True)
init_printing()









    



IPython console for SymPy 0.7.3 (Python 2.7.2-32-bit) (ground types: python)

Associating Fluids

Yields a Helmholtz energy term in the form

$$\alpha^{ass} = A\left(\ln X-\frac{X}{2}+\frac{1}{2} \right) $$

where

$$X = \frac{2}{\sqrt{1+4\bar\Delta\delta}+1}$$$$\bar\Delta(\delta,\tau) = g(\eta(\delta))[\exp(\bar\varepsilon\tau)-1]\bar\kappa$$$$\eta = \bar v_n\delta $$$$ g(\eta) = \frac{1}{2}\frac{2-\eta}{(1-\eta)^3} $$

Derivatives of $g$ and $\bar\Delta$ with respect to $\tau$ and $\delta$

$$ \frac{dg}{d\eta} = 0.5\frac{5-2\eta}{(1-\eta)^4} $$$$ \frac{d^2g}{d\eta^2} = 3\frac{3-\eta}{(1-\eta)^5} $$$$ \frac{d^3g}{d\eta^3} = 6\frac{7-2\eta}{(1-\eta)^6} $$$$\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)} $$

$$\left(\frac{\partial \bar \Delta}{\partial \tau} \right)_{\delta} = g(\eta)\bar\kappa\exp(\bar\varepsilon\tau)\bar \varepsilon \qquad \mbox{(Eq. A55)}$$$$\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)} $$

$$\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)}$$$$\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)}$$

Newly derived terms

Following the patterns of the derivatives

$$\left(\frac{\partial^3 \bar \Delta}{\partial \tau^3} \right)_{\delta} = g(\eta)\bar\kappa\exp(\bar\varepsilon\tau)\bar \varepsilon^3 $$$$\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$$$$\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$$$$\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 $$

Derivatives of $X$ with respect to $\tau$ and $\delta$

$$ \left(\frac{\partial X}{\partial \bar\Delta}\right)_{\delta} = -\frac{\delta X^2}{2\bar\Delta\delta X+1} $$$$ \left(\frac{\partial X}{\partial \delta}\right)_{\bar\Delta} = -\frac{\bar\Delta X^2}{2\bar\Delta\delta X+1} $$$$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)}$$$$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)}$$$$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)}$$

$$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} +\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)}$$

$$X_{\tau\tau} = \frac{\partial X_\tau}{\partial X}\frac{\partial X}{\partial \bar\Delta}\frac{\partial \bar\Delta}{\partial \tau} +\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)}$$$$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}$$$$X_{\tau\tau} = X_{\tau\tau}(X(\bar\Delta(\tau)),\bar\Delta,\bar\Delta_\tau,\bar\Delta_{\tau\tau})$$

Third order partials

Based on derivatives of $X_{\tau\tau}$

$$X_{\tau\tau\tau} = \frac{\partial X_{\tau\tau}}{\partial X}\frac{\partial X}{\partial \bar\Delta}\frac{\partial \bar\Delta}{\partial \tau} +\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}$$$$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} +\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}$$



In [31]:

    
X,delta,Delta,Delta_t,Delta_tt = symbols('X,delta,Delta,Delta_t,Delta_tt')
X_tt = (2*Delta_t**2*delta**2*X**3*(Delta*delta*X+1)/(2*Delta*delta*X+1)**3
        +(2*delta**2*X**3*Delta_t**2)/(2*Delta*delta*X+1)**2
        -(delta*X**2*Delta_tt)/(2*Delta*delta*X+1)
        )
print 'dX = '+ccode(simplify(diff(X_tt,X)))
print 'ddelta = '+ccode(simplify(diff(X_tt,delta)))
print 'dDelta = '+ccode(simplify(diff(X_tt,Delta)))
print 'dDelta_t = '+ccode(simplify(diff(X_tt,Delta_t)))
print 'dDelta_tt = '+ccode(simplify(diff(X_tt,Delta_tt)))









    



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)
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)
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)
dDelta_t = 4*Delta_t*pow(X, 3)*pow(delta, 2)*(3*Delta*X*delta + 2)/pow(2*Delta*X*delta + 1, 3)
dDelta_tt = -pow(X, 2)*delta/(2*Delta*X*delta + 1)

Based on derivatives of $X_{\delta\delta}$

$$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} +\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}$$$$X_{\delta\delta\tau} = \frac{\partial X_{\delta\delta}}{\partial X}\frac{\partial X}{\partial \bar\Delta}\frac{\partial \bar\Delta}{\partial \tau} +\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}$$



In [32]:

    
X,delta,Delta,Delta_d,Delta_dd = symbols('X,delta,Delta,Delta_d,Delta_dd')
            
X_dd = (X**2*(2*Delta**2*X-Delta_d)/(2*Delta*delta*X+1)**2
        -(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+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
        +X**2*(2*delta**2*X*Delta_d-1)/(2*Delta*delta*X+1)**2*Delta_d
        -delta*X**2/(2*Delta*delta*X+1)*Delta_dd
        )

print 'dX = '+ccode(simplify(diff(X_dd,X)))
print 'ddelta = '+ccode(simplify(diff(X_dd,delta)))
print 'dDelta = '+ccode(simplify(diff(X_dd,Delta)))
print 'dDelta_d = '+ccode(simplify(diff(X_dd,Delta_d)))
print 'dDelta_dd = '+ccode(simplify(diff(X_dd,Delta_dd)))









    



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)
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)
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)
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)
dDelta_dd = -pow(X, 2)*delta/(2*Delta*X*delta + 1)

Derivatives with respect to $\tau$ and $\delta$ of each term

$$ \frac{\partial \alpha^{ass}}{\partial \tau} = A\left[\frac{1}{X}-\frac{1}{2}\right]X_{\tau} \quad \quad \mathrm{(Eq. A52/A53)}$$

$$ \frac{\partial \alpha^{ass}}{\partial \delta} = A\left[\frac{1}{X}-\frac{1}{2}\right]X_{\delta} \quad \quad \mathrm{(Eq. A29/A30)}$$$$ \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)}$$$$ \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)}$$$$ \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)}$$

Newly Derived Terms

$$ \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} $$$$ \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} $$$$ \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} $$$$ \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} $$



In [28]:

    
from __future__ import division
from math import exp, log

A = 2

class Methanol(object):
    
    m = 0.977118832
    vbarn = 0.204481952
    kappabar = 0.148854832e-2
    epsilonbar = 5.46341463
    
    def X(self, delta, Deltabar):
        return 2/((1+4*Deltabar*delta)**0.5+1)
    
    def dX_dDeltabar__constdelta(self, delta, Deltabar):
        X = self.X(delta,Deltabar)
        return -delta*X**2/(2*Deltabar*delta*X+1)
    
    def dX_ddelta__constDeltabar(self, delta, Deltabar):
        X = self.X(delta,Deltabar)
        return -Deltabar*X**2/(2*Deltabar*delta*X+1)
    
    def dX_dtau(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        return self.dX_dDeltabar__constdelta(delta, Deltabar)*self.dDeltabar_dtau__constdelta(tau, delta)
    
    def dX_ddelta(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        return (self.dX_ddelta__constDeltabar(delta, Deltabar)
               + self.dX_dDeltabar__constdelta(delta, Deltabar)*self.dDeltabar_ddelta__consttau(tau, delta))
    
    def d2X_dtau2(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        X = self.X(delta, Deltabar)
        beta = self.dDeltabar_dtau__constdelta(tau, delta)
        d_dXdtau_dbeta = -delta*X**2/(2*Deltabar*delta*X+1)
        d_dXdtau_dDeltabar = 2*delta**2*X**3/(2*Deltabar*delta*X+1)**2*beta
        d_dXdtau_dX = -2*beta*delta*X*(Deltabar*delta*X+1)/(2*Deltabar*delta*X+1)**2
        dbeta_dtau = self.d2Deltabar_dtau2__constdelta(tau, delta)
        dX_dDeltabar = self.dX_dDeltabar__constdelta(delta, Deltabar)
        val = d_dXdtau_dX*dX_dDeltabar*beta+d_dXdtau_dDeltabar*beta+d_dXdtau_dbeta*dbeta_dtau
        #print val, self.X_tt(X,delta,Deltabar,beta,dbeta_dtau)
        return val
    
    def X_tt(self, X, delta, Delta, Delta_t, Delta_tt):
        return (2*Delta_t**2*delta**2*X**3*(Delta*delta*X+1)/(2*Delta*delta*X+1)**3
        +(2*delta**2*X**3*Delta_t**2)/(2*Delta*delta*X+1)**2
        -(delta*X**2*Delta_tt)/(2*Delta*delta*X+1)
        )
                  
    def d3X_dtau3(self, tau, delta):
        Delta = self.Deltabar(tau, delta)
        Delta_t = self.dDeltabar_dtau__constdelta(tau, delta)
        Delta_tt = self.d2Deltabar_dtau2__constdelta(tau, delta)
        dX_dDeltabar = self.dX_dDeltabar__constdelta(tau, delta)
        dDeltabar_dtau = self.dDeltabar_dtau__constdelta(tau, delta)
        dDeltabart_dtau = self.d2Deltabar_dtau2__constdelta(tau, delta)
        dDeltabartt_dtau = self.d3Deltabar_dtau3__constdelta(tau, delta)
        X = self.X(delta, Delta)
        d = delta
        t = tau
        
        dXtt_dDeltabartt = (-d*X**2)/(1+2*Delta*d*X)
        dXtt_dDeltabart = ((4*Delta_t*d**2*X**3*(Delta*d*X+1))/(2*Delta*d*X+1)**3
                           +(4*d**2*X**3*Delta_t)/(2*Delta*d*X+1)**2
                           )
        dXtt_dDeltabar = ((-12*Delta_t**2*d**3*X**4*(Delta*d*X+1))/(2*Delta*d*X+1)**4
                        -(6*d**3*X**4*Delta_t**2)/(2*Delta*d*X+1)**3
                        +(2*d**2*X**3*Delta_tt)/(2*Delta*d*X+1)**2
                        )
        dXtt_dX = (-(12*Delta_t**2*Delta*d**3*X**3*(Delta*d*X+1))/(2*Delta*d*X+1)**4
                   +2*d**2*X**2*(3*Delta_t**2+Delta_tt*Delta)/(2*Delta*d*X+1)**2
                   -(2*d*X*Delta_tt)/(2*Delta*d*X+1)
                   )
        
        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)
        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)
        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)
        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)
        
        return (dXtt_dX*dX_dDeltabar*dDeltabar_dtau
                +dXtt_dDeltabar*dDeltabar_dtau
                +dXtt_dDeltabart*dDeltabart_dtau
                +dXtt_dDeltabartt*dDeltabartt_dtau)
    
    def d2X_ddeltadtau(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        X = self.X(delta, Deltabar)
        alpha = self.dDeltabar_ddelta__consttau(tau, delta)
        beta = self.dDeltabar_dtau__constdelta(tau, delta)
        dalpha_dtau = self.d2Deltabar_ddelta_dtau(tau, delta)
        d_dXddelta_dDeltabar = X**2*(2*delta**2*X*alpha-1)/(2*Deltabar*delta*X+1)**2
        d_dXddelta_dalpha = -delta*X**2/(2*Deltabar*delta*X+1)
        d_dXddelta_dX = -(Deltabar+delta*alpha)*2*(Deltabar*delta*X**2+X)/(2*Deltabar*delta*X+1)**2
        dX_dDeltabar = self.dX_dDeltabar__constdelta(delta, Deltabar)
        return d_dXddelta_dX*dX_dDeltabar*beta+d_dXddelta_dDeltabar*beta+d_dXddelta_dalpha*dalpha_dtau
    
    def d2X_ddelta2(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        X = self.X(delta, Deltabar)
        alpha = self.dDeltabar_ddelta__consttau(tau, delta)
        dalpha_ddelta = self.d2Deltabar_ddelta2__consttau(tau, delta)
        d_dXddelta_dDeltabar = X**2*(2*delta**2*X*alpha-1)/(2*Deltabar*delta*X+1)**2
        d_dXddelta_dalpha = -delta*X**2/(2*Deltabar*delta*X+1)
        d_dXddelta_dX = -(Deltabar+delta*alpha)*2*(Deltabar*delta*X**2+X)/(2*Deltabar*delta*X+1)**2
        dX_dDeltabar = self.dX_dDeltabar__constdelta(delta, Deltabar)
        dX_ddelta_constall = X**2*(2*Deltabar**2*X-alpha)/(2*Deltabar*delta*X+1)**2
        return (dX_ddelta_constall
                + d_dXddelta_dX*self.dX_ddelta__constDeltabar(delta, Deltabar)
                + d_dXddelta_dX*dX_dDeltabar*alpha
                + d_dXddelta_dDeltabar*alpha
                + d_dXddelta_dalpha*dalpha_ddelta)
        
    def g(self, eta):
        return 0.5*(2-eta)/(1-eta)**3
    
    def dg_deta(self, eta):
        return 0.5*(5-2*eta)/(1-eta)**4
    
    def d2g_deta2(self, eta):
        return 3*(3-eta)/(1-eta)**5
    
    def d3g_deta3(self, eta):
        return 6*(7-2*eta)/(1-eta)**6
    
    def eta(self, delta):
        return self.vbarn*delta
    
    def Deltabar(self, tau, delta):
        return self.g(self.eta(delta))*(exp(self.epsilonbar*tau)-1)*self.kappabar
    
    def dDeltabar_ddelta__consttau(self, tau, delta):
        return self.dg_deta(self.eta(delta))*(exp(self.epsilonbar*tau)-1)*self.kappabar*self.vbarn
    
    def d2Deltabar_ddelta2__consttau(self, tau, delta):
        return self.d2g_deta2(self.eta(delta))*(exp(self.epsilonbar*tau)-1)*self.kappabar*self.vbarn**2
    
    def dDeltabar_dtau__constdelta(self, tau, delta):
        return self.g(self.eta(delta))*self.kappabar*exp(self.epsilonbar*tau)*self.epsilonbar
    
    def d2Deltabar_dtau2__constdelta(self, tau, delta):
        return self.g(self.eta(delta))*self.kappabar*exp(self.epsilonbar*tau)*self.epsilonbar**2
    
    def d3Deltabar_dtau3__constdelta(self, tau, delta):
        return self.g(self.eta(delta))*self.kappabar*exp(self.epsilonbar*tau)*self.epsilonbar**3
    
    def d2Deltabar_ddelta_dtau(self, tau, delta):
        return self.dg_deta(self.eta(delta))*exp(self.epsilonbar*tau)*self.epsilonbar*self.kappabar*self.vbarn
    
    def base(self, tau, delta):
        X = self.X(delta, self.Deltabar(tau, delta))
        return 2*(log(X)-X/2.0+0.5)
    
    def dDelta(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        X = self.X(delta, Deltabar)
        return 2*(1/X-0.5)*self.dX_ddelta(tau, delta)
    
    def dTau(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        X = self.X(delta, Deltabar)
        return 2*(1/X-0.5)*self.dX_dtau(tau, delta)
    
    def dTau2(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        X = self.X(delta, Deltabar)
        X_tau = self.dX_dtau(tau, delta)
        X_tautau = self.d2X_dtau2(tau, delta)
        return 2*(1/X-0.5)*X_tautau-2*(X_tau/X)**2
    
    def dDeltadTau(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        X = self.X(delta, Deltabar)
        X_delta = self.dX_ddelta(tau, delta)
        X_tau = self.dX_dtau(tau, delta)
        X_deltatau = self.d2X_ddeltadtau(tau, delta)
        return 2*(-X_tau/X**2)*X_delta+2*X_deltatau*(1/X-0.5)
    
    def dDelta2(self, tau, delta):
        Deltabar = self.Deltabar(tau, delta)
        X = self.X(delta, Deltabar)
        X_delta = self.dX_ddelta(tau, delta)
        X_deltadelta = self.d2X_ddelta2(tau, delta)
        return 2*(1/X-0.5)*X_deltadelta-2*(X_delta/X)**2
    
    def dTau3(self, tau, delta):
        X = self.X(delta, Deltabar)
        X_t = self.dX_dtau(tau, delta)
        X_tt = self.d2X_dtau2(tau, delta)
        X_ttt = self.d3X_dtau3(tau, delta)
        return (A*(1/X-1/2)*X_ttt
                +A*(-X_t/X**2)*X_tt
                -2*A*(X**2*(X_t*X_tt)-X_t**2*(X*X_t))/X**4
                )
    
M = Methanol()

#print M.base(0.5,0.5)
Deltabar = M.Deltabar(0.5,0.5)
#print M.dX_dDeltabar__constdelta(0.5,Deltabar),(M.X(0.5,Deltabar+1e-8)-M.X(0.5,Deltabar))/1e-8

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
print 'd',M.dDelta(0.5,0.5), (M.base(0.5,0.5+1e-8)-M.base(0.5,0.5))/1e-8
print 't',M.dTau(0.5,0.5), (M.base(0.5+1e-8,0.5)-M.base(0.5,0.5))/1e-8
print 'tt',M.dTau2(0.5,0.5), (M.dTau(0.5+1e-8,0.5)-M.dTau(0.5,0.5))/1e-8
print 'ttt',M.dTau3(0.5,0.5), (M.dTau2(0.5+1e-8,0.5)-M.dTau2(0.5,0.5))/1e-8
print 'dt',M.dDeltadTau(0.5,0.5), (M.dTau(0.5,0.5+1e-8)-M.dTau(0.5,0.5))/1e-8
print 'dd',M.dDelta2(0.5,0.5), (M.dDelta(0.5,0.5+1e-8)-M.dDelta(0.5,0.5))/1e-8









    



ttt 1 -0.473370012679 -0.473370015408
2 0.295697260234 0.295697258346
3 0.354430230606 0.354430232385
4 -0.811327611566 -0.77524633102
-1.95954772419 -1.93168185336
d -0.0351186626009 -0.0351186635328
t -0.0796836538477 -0.0796836374661
tt -0.422816696531 -0.422816709422
ttt 1 -0.473370012679 -0.473370015408
2 0.295697260234 0.295697258346
3 0.354430230606 0.354430232385
4 -0.811327611566 -0.77524633102
-2.20666413982 -2.17802790803
dt -0.199708453268 -0.199708456716
dd -0.0354391474199 -0.0354391481439



In [ ]: