fix return type precision for n_thermal

Author: ytdHuang

src/utilities.jl

@@ -42,10 +42,10 @@ n(\omega, \omega_{\textrm{th}}) = \frac{1}{e^{\omega/\omega_{\textrm{th}}} - 1},
 ```
 where ``\hbar`` is the reduced Planck constant, and ``k_B`` is the Boltzmann constant.
 """
-function n_thermal(ω::Real, ω_th::Real)::Float64
+function n_thermal(ω::T1, ω_th::T2) where {T1<:Real,T2<:Real}
     x = exp(ω / ω_th)
-    n = ((x != 1) && (ω_th > 0)) ? (1.0 / (x - 1.0)) : 0.0
-    return n
+    n = ((x != 1) && (ω_th > 0)) ? 1 / (x - 1) : 0
+    return _FType(promote_type(T1, T2))(n)
 end
 
 _get_dense_similar(A::AbstractArray, args...) = similar(A, args...)