The precision at long wavelengths L is limited by the calculation of exp(x)-1, where x = hk/TL. 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 getBBfluxPerMetre(angstromValues, kelvin, powerValues)
real(double), intent(in) :: angstromValues(:), kelvin
real(double), intent(out) :: powerValues(size(angstromValues))
end subroutine getBBfluxPerMetre
function bbFluxPerMetre(angstroms, kelvin) result(power)
real(double), intent(in) :: angstroms, kelvin
real(double) :: power
end function bbFluxPerMetre