Source code for burnman.eos.einstein

# 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.


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(nopython=True):
        def decorator(fn):
            return fn

        return decorator


"""
Functions for the Einstein model of a solid.
"""

eps = np.finfo(float).eps




@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





@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





@jit(nopython=True)
def helmholtz_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





@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





@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





@jit(nopython=True)
def dentropy_dTheta(T, einstein_T, n):
    """
    Derivative of Einstein entropy with respect to Einstein temperature [J/K^2]
    """
    if T <= eps:
        return 0.0
    x = einstein_T / T
    f = np.exp(x) - 1.0
    return -3.0 * n * constants.gas_constant * (x * np.exp(x) / (f * f * T))





@jit(nopython=True)
def dhelmholtz_dTheta(T, einstein_T, n):
    """
    First derivative of Helmholtz free energy with respect to Einstein temperature [J/K]
    """
    if T <= eps:
        return 0.0
    x = einstein_T / T
    f = np.exp(x) - 1.0
    return 3.0 * n * constants.gas_constant / f





@jit(nopython=True)
def d2helmholtz_dTheta2(T, einstein_T, n):
    """
    Second derivative of Helmholtz free energy with respect to Einstein temperature [J/K^2]
    """
    if T <= eps:
        return 0.0
    x = einstein_T / T
    f = np.exp(x) - 1.0
    return -3.0 * n * constants.gas_constant * np.exp(x) / (f * f * T)