I'm working on a small programming project involving Planck's Law, and I keep getting errors. I'm fairly certain this is due to a misunderstanding of physics on my behalf.

Basically, I am trying to do numerical integration over a range of wavelengths (1nm to 10 000nm) for Planck's Law. I have a function written for this, and I am comparing the results with the Stefan Boltzmann Law results for the same temperatures.

I am working under the assumption that the radiance over the range of wavelengths from 1nm to 10 000nm will more or less approximate the total spectrum, as spectral radiance approaches zero quite quickly after ~5000nm. Is this a valid assumption...? I've been comparing the error ratios (Stefan Boltzmann "true" value / integrated "approximate" value) between the two functions (over a wide range of temperatures) and here is what they look like:

  • 500K: 3.9795
  • 1000K: 2.87267
  • 1500K: 2.73022
  • 2000K: 2.69029
  • ...
  • 9500K: 2.65730
  • 10000K: 2.65723

As you can see, the Stefan Boltzmann Law result converges to be around 2.65 times higher than the integrated result as temperature increases. My theory on this is that I am severely under-estimating the amount of radiance accounted for in wavelengths between 10000nm and as it goes to infinity. As seen in this graph:

As the temperature gets lower, the radiance in the range I'm checking (1nm to 10 000nm) presumably becomes a lower overall part of the total spectrum radiance. I think this may explain my higher ratios (i.e lower denominator, which is my approximate integration function) at lower temperatures? I cannot check this however, as I am not prepared to numerically integrate over an infinite range of wavelengths. Am I correct on this? If so, how can I (preferably computationally) conduct a definite integration of Planck's Law over a range of wavelengths?



#include "Sunlight.h"
#include "Constants.h"
#include "si.h"

#include <iostream>

double Sunlight::spectrumPLNK(int temp_body, double wavelength){
    double spect_rad = ((2 * planck_c * pow(speed_light_c, 2)) / pow(wavelength, 5)) * 1 / (expm1((planck_c*speed_light_c) / (wavelength * boltzmann_c * temp_body)));
return spect_rad;

double Sunlight::integratePLNK(int temp_body, double n, double wave_len_low, double wave_len_high){
    double power = 0;

    double a = wave_len_low;
    double b = wave_len_high;

    double h = (b-a)/n;

    for(double i = 0; i<n; i++){
        power += spectrumPLNK(temp_body, a + i*h) + spectrumPLNK(temp_body, a + (i+1));
        std::cout << "New power: " << power << std::endl;

    power *= h;

    std::cout << "Total spectrum power: " << power << std::endl;
    return power;


double Sunlight::powerSBZL(int temp_body){
    double power_body = sigma_c* pow(temp_body, 4);

    std::cout << "Total SBL power: " << power_body << std::endl;
    return power_body;


double Sunlight::integratePLNK(int temp_body, double n, double wave_len_low, double wave_len_high){
    double power = 0;

    double a = wave_len_low;
    double b = wave_len_high;

    double h = (b-a)/n;

    for(double i = 0; i<n; i++){
        power += spectrumPLNK(temp_body, a + (i+1)*h);

    power *= pi_c*h;

    std::cout << "Total spectrum power: " << power << std::endl;
    return power;
  • $\begingroup$ Your problem with the proportionality constant might be in the way you are trying to integrate over the wavelength, rather than the frequency. The volume element that is proportional to the energy in a band is proportional to $df$, not $d\lambda$. The transformation from frequency to wavelength will also distort the integration at the short end of the spectrum, I believe, because your sampling intervals are getting very crude. $\endgroup$
    – CuriousOne
    Commented Jun 28, 2015 at 2:38
  • $\begingroup$ What is the integration scheme you are using? It seems you are evaluating at a+i*h and a+i+1 for i=0:n. That second part makes no sense at all. Did you mean a + (i+1)*h perhaps? Even that is not very good - but definitely better... I would highly recommend using something other than trapezoid - even Simpson's rule (1-4-2-4-2-...4-1) will deal better with the curvature than (1-2-2-2-2-...2-1). $\endgroup$
    – Floris
    Commented Jun 28, 2015 at 2:41
  • $\begingroup$ Sorry @Floris, I have changed my integration method multiple times (trying to see if using better methods would yield substantially more accurate results). This was supposed to be trapezoid I believe. I tried Simpson's Rule earlier and there was no real difference in accuracy that I could see. $\endgroup$
    – Brad Power
    Commented Jun 28, 2015 at 2:51
  • $\begingroup$ Do you see the error - you are using (i+1) instead of (i+1)*h... $\endgroup$
    – Floris
    Commented Jun 28, 2015 at 3:06
  • $\begingroup$ Yes, I changed that. However, the error between my functions persists. This is odd, I feel as if I am missing a factor or something somewhere as you have mentioned in the comments on your answer. $\endgroup$
    – Brad Power
    Commented Jun 28, 2015 at 3:13

1 Answer 1


Without seeing your code, it is hard to know where you went wrong. I did write a few lines of code myself to see what might be going on - and found that the values converge quite nicely as you approach higher temperatures (for the lowest temperatures of 500 K and 1000 K, significant amounts of energy do indeed extend beyond 10,000 nm).

updated There is a small problem with scaling: the SB law gives radiation from a black body in W/m$^2$ which is different than Planck's Law which gives radiance per unit area, per steradian.

Now there is a modification for the SB law: if you divide the usual equation by $\pi$ you get the result in "per steradian" (see the Wiki article).

This gives me the following code snippet (Python - not my strongest language but it's good for this kind of thing because of the scipy.integrate functions which do a lot of the hard work; the built in physical constants and plotting routines are nice too):

# integrate Planck's Law and compare results with Stefan-Boltzman

from math import pow,exp,pi
from scipy.constants import codata
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt

# get physical constants we need:
D = codata.physical_constants
h = D['Planck constant'][0]
k = D['Boltzmann constant'][0]
c = D['speed of light in vacuum'][0]
s = D['Stefan-Boltzmann constant'][0]

def planck(l, T):
    # reference https://en.wikipedia.org/wiki/Planck's_law
    if (l <= 0):
        return 0  # no valid values <=0
    p = c*h/(k*l*T)
    if (p > 700):
        return 1e-99 # preventing under/overflow in the exp function
        return (2*h*c*c)/(pow(l, 5.0) * (exp(c*h/(k*l*T))-1))  # <<< fixed

def SB(T, emissivity = 1.0):
    #  Stefan-Boltzmann law : returns radiance = total power per steradian
    # reference https://en.wikipedia.org/wiki/Stefan-Boltzmann_law
    return emissivity * s * pow(T, 4.0) / pi  # <<< expression "per steradian"

# T values to use:
Tvec = np.arange(500,10500, 500)

max_wavelength = 10000*1e-9 # 10,000 nm
for T in Tvec:
    y,abserr = integrate.quad(planck, 0, max_wavelength, args=(T))

plt.plot(Tvec, np.array(Planck_integrals)/np.array(SB_integrals))
plt.xlabel('Temperature (K)')
plt.ylabel('Ratio integral/SB value')
plt.title('Comparing direct integration and Stefan-Boltzmann - integrating 0-10 um')

And this plot:

enter image description here

I played around a bit with discrete integration to see if I can get the results to match (so you don't need to use existing libraries... although that's usually a really good idea). This was the result (using the same definition of functions planck and SB as before):

# stepwise integration
n = 1000
simpson = 0.
trapezoid = 0.
max_wavelength = 10000.0
min_wavelength = 0.
H = (max_wavelength - min_wavelength)/n
wavelength_scale = 1e-9
T = 2000.0

for i in range(n):
    if ((i==0) | (i == n-1)):
        w_trap = 1
        w_simpson = 1
        w_trap = 2
        if (i%2==0):
            w_simpson = 4
            w_simpson = 2
    wavelength = (min_wavelength + H * i) * wavelength_scale
    hp = H * planck(wavelength, T) * wavelength_scale
    trapezoid = trapezoid + hp * w_trap
    simpson = simpson + hp * w_simpson

print 'trapezoid: ', trapezoid / 2
print 'simpson: ', simpson / 3
print 'SB: ', SB(T)

This gave the following output:

trapezoid:  284606.470099
simpson:  284606.473505
SB:  288789.72548

Clearly the two integration methods give nearly the same result; and it's very comparable with the result from the Stefan-Boltzman formula. I have to conclude the error is in your integration formula (as I pointed out in my comment), and/or in extremely coarse sampling of the function. But even when I set n=100 in the above I get comparable results...

  • $\begingroup$ Thank you! I have updated my post with my C++ code, and I have a small question about yours. I noticed in your planck() function, you return: "return (hcc)/(pow(l, 5.0) * (exp(ch/(klT))-1))". Isn't the Planck Function supposed to be (2*hc*c), or am I missing something? $\endgroup$
    – Brad Power
    Commented Jun 28, 2015 at 2:44
  • $\begingroup$ You are right - I will look more carefully and update. I was a bit puzzled about the $2\pi$ factor. In fact there is a "SB per steradian" function that has a factor $\pi$. Combined with the factor 2 you spotted we can clear this up. I will also include some comments on your code in the update. $\endgroup$
    – Floris
    Commented Jun 28, 2015 at 2:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.