#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()