#encoding: utf8
from __future__ import division
from numpy import array, atleast_1d, pi, exp, vectorize, sqrt, log
from scipy.optimize import leastsq
from scipy import interpolate

class Adeb3(object):
    left = -1
    right = 1
    coefficients = [
        2.707737068327440945,
        0.340068135211091751,
        -0.12945150184440869e-01,
        0.7963755380173816e-03,
        -0.546360009590824e-04,
        0.39243019598805e-05,
        -0.2894032823539e-06,
        0.217317613962e-07,
        -0.16542099950e-08,
        0.1272796189e-09,
        -0.987963460e-11,
        0.7725074e-12,
        -0.607797e-13,
        0.48076e-14,
        -0.3820e-15,
        0.305e-16,
        -0.24e-17]

def cheb_eval(cs, x):
    d = 0
    dd = 0
    y = (2 * x - cs.right - cs.left) / (cs.right - cs.left)
    y2 = 2 * y

    for c in cs.coefficients[:0:-1]:
        d, dd = y2 * d - dd + c, d
    return y * d - dd + 0.5 * cs.coefficients[0]

@vectorize
def debye3(x):
    val_infinity = 19.4818182068004875
    xcut = 7.0839641853226408e+02

    if x < 0:
        raise ValueError("x must be > 0!")
    elif x < 2.0 * sqrt(2) * 1.5e-8:
        return 1 - (3 / 8) * x + x ** 2 / 20
    elif x <= 4:
        t = x ** 2 / 8 - 1
        c = cheb_eval(Adeb3, t)
        return c - 0.375 * x
    elif x < 36 - log(2):
        nexp = int(xcut / x)
        ex = exp(-x)
        sum = 0

        for rk in range(nexp, 0, -1):
            xk_inv = 1 / (rk * x)
            sum *= ex
            sum += (((6 * xk_inv + 6) * xk_inv + 3) * xk_inv + 1) / rk
        return val_infinity / x ** 3 - 3 * sum * ex
    elif x < xcut:
        sum = ((x + 3) * x + 6) * x + 6
        return (val_infinity - 3 * sum * exp(-x)) / x ** 3
    else:
        return val_infinity / x ** 3

def capacity(x):
    return 4 * debye3(x) - 3 * x / (exp(x) - 1)

#wikipedia
m = 157.25 # g / mol
density = 7.886e6 # g / m^3

# following Tsang et al., PRB 31, 235 (1985)
#td = 163.4 # K
#gamma = 6.38e-3 # J / mol K^2

# Lounasmaa et al., PR 150, 399 (1966) (= Tb)
#gamma = 10.5e-3 # J / mol K^2

# Hill et al., J Phys F 17, 1867 (1987)
gamma = 4.48e-3 # J / mol K^2
td = 169 # K

# CRC Handbook, 87 ed
# http://ingemeca.org/docs/Genie_chimique/Handbook%20of%20chemistry%20and%20physics/Section%2004/04_03_86.pdf
gamma = 4.48e-3 # J / mol K^2
td = 169 # K
lamda = 0.30 

# from wikipedia
# very large error bars since temperature dependent
kappa = 11 # W / m K

# Palik (ed.), Handbook of Optical Constants of Solids, Academic Press, 1998
# very large error bars
n = 2.5+2.5j

# from wikipedia, Dulong-Petit
gasR = 8.314472 # J / mol K

# from Lewis, PRB 1, 4368 (1970)
tc = 291 # K
alpha = -0.09
alphap = -0.32

cal = 4.184 # J

# from Griffel et al., PR 93, 657 (1954)
# data deleted for copyright reasons
#data = [ ]
#data = array(data)
#t = data[:, 0]
#
#rdata = data[:, 1] * cal - (gamma * t + capacity(td / t) * 3 * gasR)
#
#spline1 = interpolate.splrep(t[t < tc], rdata[t < tc], s=0.01)
#spline2 = interpolate.splrep(t[t >= tc], rdata[t >= tc], s=0.01)

# the resulte of the above spline interpolation:
spline1 = (array([15., 15., 15., 15., 35., 50., 60., 65., 70.,
         80.,   85.,  105.,  120.,  155.,  175.,  190.,  210.,  225.,
        260.,  275.,  285.,  290.,  290.,  290.,  290.]),
      array([0.49389679, 1.27420175, 2.77687614, 4.02375527,
         4.76956233, 5.22356252, 5.59837477, 5.92009627,
         6.49799508, 6.85174211, 7.75186406, 8.48986024,
         9.70070657, 10.58583945, 11.7457047, 13.56042647,
        16.1996443, 20.93242938, 24.73094546, 28.80272272,
        30.5090826, 0., 0., 0., 0.]), 3)
spline2 = (array([ 295.,  295.,  295.,  295.,  305.,  310.,  325., 340., 355.,
        355., 355., 355.]),
  array([ 13.45959205,  10.78000421,   9.00637154,   7.23686324,
         6.11750324,   5.3514927 ,   4.38785726,   4.33129557,
         0.        ,   0.        ,   0.        ,   0.        ]), 3)

def critical(t, alpha, a, b):
    return b + a * abs((t - tc) / tc) ** -alpha

left = [-1, 1]
right = [1, 1]

def spincapacity(t): # in J / mol K
    t = array(t)
    return (t > tc).choose(interpolate.splev(t, spline1),
        interpolate.splev(t, spline2))

def electroncapacity(t): # in J / mol K
    return gamma * t

def phononcapacity(t): # in J / mol K
    return capacity(td / t) * 3 * gasR

if __name__ == "__main__":
    from pylab import plot, show, xlabel, ylabel, savefig
    from numpy import linspace
    tt = linspace(15, 355, 300)
    xlabel("Temperatur (K)")
    ylabel(u"Wärmekapazität (J / mol)")
    plot(tt, spincapacity(tt))
    plot(tt, spincapacity(tt) + phononcapacity(tt) + electroncapacity(tt))
    tt = linspace(0, 355, 300)
    plot(tt, phononcapacity(tt))
    plot(tt, electroncapacity(tt))

    energy = 4 * pi * n.imag * 1.2 / 780e-7 # mJ / cm^3
    energy *= 1e3 # J / m^3
    energy *= m / density # J / mol
    t = 100 # K

    savefig("gadolinium.svg")
    show()