Find the Planck Curve fraction from wavelengths 0 to λ

(c) 2013 by Barton Paul Levenson

Language:
C#

Dialect:
Microsoft Visual C# Express 2010

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 C# 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 λ₂ > λ₁.

        // D returns the Planck fraction from 0 to a given wavelength.
        static double D(double lambda, double T)
        {
            const double c2 = 14387.77;
            const double pi = 3.141592653589793;
            const double k = 15.0 / (pi * pi * pi * pi);

            double q = c2 / (lambda * T);
            double sum = 0.0;

            for (int m = 1; m < 200; m++)
            {
                double f = m * q;
                double term1 = Math.Exp(-f) / (m * m * m * m);
                double term2 = f * (f * (f + 3.0) + 6.0) + 6.0;
                sum += (term1 * term2);
            }
            return k * sum;
        } // D