The precision at low frequencies f is limited by the calculation of exp(x)-1, where x = hfk/T. At low values of x, exp(x) is close to 1 and thus the difference between exp(x) and 1 is a number of low precision. Thus for x values less than 0.1, a series expansion of (exp(x)-1)/x is used instead. The minimum precision occurs at x = 0.1 and is approximately equal to precision(1d0)-1.
subroutine getBBfluxPerHertz(hertzValues, kelvin, powerValues)
real(double), intent(in) :: hertzValues(:), kelvin
real(double), intent(out) :: powerValues(size(hertzValues))
end subroutine getBBfluxPerHertz
function bbFluxPerHertz(hertz, kelvin) result(power)
real(double), intent(in) :: hertz, kelvin
real(double) :: power
end function bbFluxPerHertz