In [1]:

    

from sympy import *

# symbolic variables for Helmholtz energy
fi, fr, delta, tau = symbols("fi fr delta tau")

# symbolic variables for parameters in bank of terms
li1, li2, ei1, ei2, ci1, ci2, si1, si2 = symbols("li1 li2 ei1 ei2 ci1 ci2 si1 si2")
pi1, pi2, pi3 = symbols("pi1 pi2 pi3")
bi1, bi2, bi3, bi4 = symbols("bi1 bi2 bi3 bi4")
gi1, gi2, gi3, gi4, gi5, gi6, gi7, gi8, gi9 = symbols("gi1 gi2 gi3 gi4 gi5 gi6 gi7 gi8 gi9")
ni1, ni2, ni3, ni4, ni5, ni6, ni7, ni8, ni9, ni10, ni11, ni12 = symbols("ni1 ni2 ni3 ni4 ni5 ni6 ni7 ni8 ni9 ni10 ni11 ni12")
Distance, Phi = symbols("Distance Phi")

# terms for ideal Helmholtz energy
logTerms = + li1*log(tau**li2);
idealPolyTerms = + pi1*tau**pi2;
EinsteinTerms = + ei1*log(1 - exp(ei2*tau));
coshTerms = - ci1*log((cosh(ci2*tau)));
sinhTerms = + si1*log((sinh(si2*tau)));

# terms for residual Helmholtz energy
residualPolyTerms = pi1*tau**pi2*delta**pi3;
BWRTerms = bi1*tau**bi2*delta**bi3*exp(-delta**bi4);
GaussTerms = gi1*tau**gi2*delta**gi3*exp(gi6*(delta-gi9)**2 + gi7*(tau-gi8)**2);
Distance = ((1-tau)+ni8*((delta-1)**2)**(1/(2*ni7)))**2 + ni11*((delta-1)**2)**ni12
Phi = exp(-ni9*(delta-1)**2-ni10*(tau-1)**2)
nonAnalyticalTerms = ni1*delta*Distance**ni6*Phi;


    

In [2]:

    

# Set the Helmholtz energy we want to look at
# fi = +log(delta) +logTerms +idealPolyTerms +EinsteinTerms +coshTerms +sinhTerms
fi = logTerms + coshTerms + sinhTerms
# fr = +residualPolyTerms +BWRTerms +GaussTerms +nonAnalyticalTerms
fr = nonAnalyticalTerms

# and print terms, just to make sure
print(fi)
print(fr)


    

    




-ci1*log(cosh(ci2*tau)) + li1*log(tau**li2) + si1*log(sinh(si2*tau))
delta*ni1*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)**ni6*exp(-ni10*(tau - 1)**2 - ni9*(delta - 1)**2)

In [3]:

    

## First derivs
# Calculate
fit = simplify(diff(fi,tau))
frt = simplify(diff(fr,tau))
frd = simplify(diff(fr,delta))

## Print results
print(fit)
print(frt)
print(frd)


    

    




-ci1*ci2*tanh(ci2*tau) + li1*li2/tau + si1*si2/tanh(si2*tau)
-2*delta*ni1*(ni10*(tau - 1)*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)**(ni6 + 1) + ni6*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)**ni6*(ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1))*exp(-ni10*(tau - 1)**2 - ni9*(delta - 1)**2)/(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)
ni1*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)**(ni6 - 1)*(2*delta*ni6*(ni11*ni12*ni7*((delta - 1)**2)**ni12 + ni8*(ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)*((delta - 1)**2)**(1/(2*ni7))) + ni7*(delta - 1)*(ni11*((delta - 1)**2)**ni12 + (ni8*((delta - 1)**2)**(1/(2*ni7)) - tau + 1)**2)*(-2*delta*ni9*(delta - 1) + 1))*exp(-ni10*(tau - 1)**2 - ni9*(delta - 1)**2)/(ni7*(delta - 1))

In [ ]:

    

# Second derivs
# Calculate
fitt = simplify(diff(fit,tau))
frtt = simplify(diff(frt,tau))
frtd = simplify(diff(frt,delta))
frdd = simplify(diff(frd,delta))

# Print results
print(fitt)
print(frtt)
print(frtd)
print(frdd)


    

In [ ]:

    

# Third derivs
#Calculate
fittt = simplify(diff(fitt,tau))
frttt = simplify(diff(frtt,tau))
frttd = simplify(diff(frtd,tau))
frtdd = simplify(diff(frtd,delta))
frddd = simplify(diff(frdd,delta))

# Print results
print(fittt)
print(frttt)
print(frttd)
print(frtdd
print(frddd)