debye.py

debye.py

from __future__ import print_function

import numpy as np
from scipy import integrate

# scipy.integrate
# https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html

# scipy.integrate.quad
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad


def integrand(x):
    if x==0:
        return 0
    return x**4 / np.sinh(x/2.)**2


def specific_heat(t, verbose=False):
    if t==0:
        return 0
    y, err = integrate.quad(integrand, 0, 1./t)
    if verbose:
        print(f"integral : {y}")
        print(f"error    : {err}")
    return (3./4.) * t**3 * y


def main():
    t = 0.1
    c = specific_heat(t, verbose=True)
    print(f"T = {t}")
    print(f"c = {c}")


if __name__ == '__main__':
    main()