pub fn _4(x: f64) -> f64
This routine computes the fourth-order Debye function D_4(x) = (4/x^4) \int_0^x dt (t^4/(e^t - 1)).