bec.py

bec.py

from __future__ import print_function

import numpy as np
from scipy import integrate
from scipy import special
from scipy import optimize

# 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

# scipy.special.gamma
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gamma.html

# scipy.optimize.root_scalar
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root_scalar.html#scipy.optimize.root_scalar


def integrand(x, s, alpha):
    # x+alpha > 0
    e = np.exp(-x-alpha)
    return np.power(x, s-1) * e / (1.0 - e)


def bose_einstein_integral(s, alpha, verbose=False):
    y, err = integrate.quad(integrand, 0, np.inf, args=(s, alpha))
    if verbose:
        print(f"integral : {y}")
        print(f"error    : {err}")
    return y / special.gamma(s)


def func_alpha(alpha, t):
    return bose_einstein_integral(1.5, alpha) * np.power(t, 1.5) - special.zeta(1.5)


def calc_alpha(t, verbose=False):
    if t<=1:
        return 0

    a_min = 0
    a_max = 100
    sol = optimize.root_scalar(func_alpha, args=(t,), method='brentq', bracket=(a_min, a_max))
    if verbose:
        print(sol)
    return sol.root


def calc_specific_heat(alpha, t):
    if t<=1:
        z52 = special.zeta(2.5)
        z32 = special.zeta(1.5)
        return 15. * z52 / (4. * z32) * np.power(t, 1.5)
    else:
        f12 = bose_einstein_integral(0.5, alpha)
        f32 = bose_einstein_integral(1.5, alpha)
        f52 = bose_einstein_integral(2.5, alpha)
        return 15. * f52 / (4. * f32) - 9. * f32 / (4. * f12)


def solve_bec(t, verbose=False):
    alpha = calc_alpha(t, verbose=verbose)
    c = calc_specific_heat(alpha, t)
    return alpha, c


def main():
    t = 1.5
    alpha, c = solve_bec(t, verbose=True)
    print(f"alpha = {alpha}")
    print(f"c = {c}")


if __name__ == '__main__':
    main()