# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for the Earth and Planetary Sciences
# Copyright (C) 2012 - 2017 by the BurnMan team, released under the GNU
# GPL v2 or later.
from __future__ import absolute_import
import numpy as np
from .. import constants
# Try to import the jit from numba. If it is
# not available, just go with the standard
# python interpreter
try:
import os
if "NUMBA_DISABLE_JIT" in os.environ and int(os.environ["NUMBA_DISABLE_JIT"]) == 1:
raise ImportError("NOOOO!")
from numba import jit
except ImportError:
def jit(fn):
return fn
"""
Functions for the Einstein model of a solid.
"""
eps = np.finfo(float).eps
[docs]
@jit(nopython=True)
def thermal_energy(T, einstein_T, n):
"""
calculate the thermal energy of a substance. Takes the temperature,
the Einstein temperature, and n, the number of atoms per molecule.
Returns thermal energy in J/mol
"""
if T <= eps:
return 0.0
x = einstein_T / T
E_th = 3.0 * n * constants.gas_constant * einstein_T * (1.0 / (np.exp(x) - 1.0))
return E_th
[docs]
@jit(nopython=True)
def molar_heat_capacity_v(T, einstein_T, n):
"""
Heat capacity at constant volume. In J/K/mol
"""
if T <= eps:
return 0.0
x = einstein_T / T
C_v = (
3.0
* n
* constants.gas_constant
* (x * x * np.exp(x) / np.power(np.exp(x) - 1.0, 2.0))
)
return C_v
[docs]
@jit(nopython=True)
def helmholtz_free_energy(T, einstein_T, n):
"""
Helmholtz free energy of lattice vibrations in the Einstein model [J].
It is important to note that this does NOT include the zero
point energy for the lattice. As long as you are
calculating relative differences in F, this should cancel anyway.
"""
E = thermal_energy(T, einstein_T, n)
S = entropy(T, einstein_T, n)
return E - T * S
[docs]
@jit(nopython=True)
def entropy(T, einstein_T, n):
"""
Entropy due to lattice vibrations in the Einstein model [J/K]
"""
if T <= eps:
return 0.0
x = einstein_T / T
S = (
3.0
* n
* constants.gas_constant
* (-x * np.exp(-x) / (np.exp(-x) - 1.0) - np.log(1.0 - np.exp(-x)))
)
return S
[docs]
@jit(nopython=True)
def dmolar_heat_capacity_v_dT(T, einstein_T, n):
"""
First temperature derivative of the heat capacity at constant volume
according to the Einstein model [J/K^2/mol].
"""
if T <= eps:
return 0.0
x = einstein_T / T
dCvdT = (
3.0
* n
* constants.gas_constant
* x
* x
* np.exp(x)
* ((x - 2.0) * np.exp(x) + (x + 2.0))
/ (T * np.power(np.exp(x) - 1.0, 3.0))
)
return dCvdT