statistical mechanicsthermal-radiationthermodynamics
In the above image, the curves for different temperature dont intersect anywhere. Stefan-Boltzmann law and Wien displacement law dont preclude the intersection.
Is it because if, for example, they cross each other at some longer wavelength somewhere (after the peak), then that would imply that the body with the lower temperature emits more power than that with the higher temperature for the wavelengths after the crossing point? If yes, then why is that not possible?
The answer can be calculated by referring to band emission tables of any introductory heat transfer book. I use the heat transfer book by Incropera and DeWitt where are concept is explained and the tables given in chapter 12.
Here is how its done. The fraction of total energy emitted as wavelengths from 0 to a certain wavelength $\lambda$ can be obtained from the band emission table ($F_{(0-\lambda)}$). This way the emissions from 0 to 400 nm ($F_{(0-0.4\mu m)}$) and between 0 to 800 nm ($F_{(0-0.8\mu m)}$) can be obtained. The difference in these two quantities gives the emission from 400 nm to 800 nm and is the required answer.
The energy fractions $F_{(0-\lambda)}$ as a function of $\lambda$ and $T$ are given in the table
Here $T$ = 2,573 K, while $\lambda$ is the upper limit of wavelength.
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:
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...
Best Answer
Here are two reasons they can't cross, one from statistical mechanics and one from pure thermo.
A more direct proof is to go all the way to Planck's law, but that's not necessary: we can prove the curves can't cross with general principles, without knowing much about the details.