pub fn bose_einstein(x: f64, mu: f64) -> f64
Bose-Einstein distribution: n_BE(x, μ) = 1/(e^(x+μ) - 1)