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