[Physics] How to find the radiance over a finite range of wavelengths using Planck’s Law

electromagnetic-radiationthermal-radiation

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?

Thanks!

EDIT 1:

#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;
}

EDIT 2:

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;
}

Best Answer

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
    else:
        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)

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

plt.figure()
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')
plt.show()

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
    else:
        w_trap = 2
        if (i%2==0):
            w_simpson = 4
        else:
            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...

Related Question