Find the Planck Curve fraction from wavelengths 0 to λ

(c) 2017 by Barton Paul Levenson

Language:
Fortran-95

Dialect:
Salford Fortran Personal Edition

Discussion:
The Planck Curve relates the "irradiance" at a given temperature and wavelength to those factors:

I_λ = c₁ / [λ⁵ [1 - exp(c₂ / {λ T})]

where c₁ and c₂ are the 1st and 2nd radiation constants, respectively. To find the fraction of irradiance between two wavelengths, then, you need only integrate this equation between those wavelengths. Unfortunately, this equation falls into a class that cannot be integrated! But an iterative approximation was worked out some time in the years 1910-1930 (I'd be grateful if anyone can point me to a primary source). The algorithm is produced below in the Fortran-95 programming language. Please note that what this routine does is give you the fraction of the Planck curve betweeen 0 and a given wavelength, known as "D." To get the fraction between two wavelengths, you need D(λ₂) - D(λ₁), where λ₂ > λ₁.

This version is 31 times faster than the C# version I posted four years ago.

! DD finds the Planck radiation fraction from 0 to lambda.
real(d) function DD(lambda, T)
    implicit   none

    real(d), parameter  :: c2  = 14387.77d0
    real(d), parameter  :: pi  = 3.141592653589793d0
    real(d), parameter  :: k   = 15d0 / (pi * pi * pi * pi)
    real(d), intent(in) :: lambda, T

    integer  :: m
    real(d)  :: f, q, sum, term1, term2
    real(d)  :: lastSum, eps
    
    lastSum = -1d0
    m       = 0
    q       = c2 / (lambda * T)
    sum     = 1d-18

    do
        m     = m + 1
        if (m .gt. 200) exit

	f     = m * q
        term1 = exp(-f) / (m * m * m * m)
    	term2 = f * (f * (f + 3d0) + 6d0) + 6d0
   	sum   = sum + term1 * term2
        eps   = abs(sum - lastSum) / sum

	if (eps .lt. 1d-6) exit
   	    lastSum = sum
        end do

	DD = k * sum
end function DD