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