Debye's Temperature calculator basing on cv/3R value

debye_calculator.py

#!/usr/bin/env python3
import numpy as np
from scipy.integrate import quad
from scipy.optimize import root_scalar

# Funkcja Debye'a D_3(y) = \int_0^y (x^3 / (exp(x) - 1)) dx
def D3(y):
    integrand = lambda x: x**3 / (np.exp(x) - 1.0)
    val, _ = quad(integrand, 0, y, limit=200)
    return val

# Równanie implicitne f(T_D) = 0
def debye_equation(TD, cv_over_R, T1, T2):
    term = (
        T2**4 * D3(TD / T2)
        - T1**4 * D3(TD / T1)
    )
    rhs = (9.0 / (cv_over_R * (T2 - T1))) * term
    return TD**3 - rhs

# Rozwiązanie numeryczne
def solve_TD(cv_over_3R, T1, T2, TD_min=1.0, TD_max=2000.0):
    cv_over_R = cv_over_3R * 3.0
    sol = root_scalar(
        debye_equation,
        args=(cv_over_R, T1, T2),
        bracket=[TD_min, TD_max],
        method='brentq'
    )

    if not sol.converged:
        raise RuntimeError("Nie udało się znaleźć T_D")

    return sol.root

if __name__ == "__main__":
    # dane eksperymentalne
    #cv_over_3R = 3*0.7343    # przykładowa wartość
    d = [
["grafit (podpisany)",0.1939,0.001186],
["grafit",0.197,0.00164],
["krzem",0.5681,0.004529],
["SiC  (walec)",0.304,0.004742],
["SiC (łuska)",0.3054,0.003076],
["WC (kulka)",0.5851,0.007655],
["German",0.7791,0.004959],
["żelazo",0.6633,0.004107],
["Ni",0.7454,0.005065],
["Zr",0.854,0.006326],
["Nb",0.838,0.006082],
["Cu (taśma)",0.6007,0.006208],
["Cu (wstążka)",0.8804,0.009904],
["Cr",0.7229,0.005013],
["Al",0.7343,0.00477],
["Co",0.61,0.00487],
["Tungsten (W)",0.8356,0.008167]
    ]
    T1 = 77.0          # K
    T2 = 295.0          # K

    for x in d:
        TD = solve_TD(x[1], T1, T2)
        TD2 = solve_TD(x[1] + x[2], T1, T2)
        uTD = np.abs((TD2 - TD))
        print(f"{x[0]}, {TD:.2f}, {uTD: 2f}")