# Implementing a black-body using Planck's Law in a 3D ray tracer

1_ In one of my software projects, I need to implement Planck's Law in terms of wavelength. I can't be sure if it's right because I've seen different variations of it, and all the implementations that I've tried give different results (besides, it's been years since I've had to do physics). Could someone verify that it's correct? Have I got the units right?

It looks something like this, and it's the implementation of Planck's Law in terms of wavelength found on this wikipedia page:

h = 6.6260695729e-34 J*s
k = 1.3806488e-23 J*K^-1
c = 2.998e8 m/s
c2 = pow(c, 2) m^2/s^2
h_c = h * c
_2_h_c2 = 2.0 * h * c2

planck(wavelength, temp) {
return _2_h_c2 / (pow(wavelength, 5.0) * (exp(h_c / (wavelength * k * temp)) - 1.0))
}


temp is the temperature in Kelvins, and wavelength is in nanometers.

2_ If I wanted it normalized about 555nm, similar to this wikipedia page, what would it look like? This is what I have so far, but, again, I can't be sure if it's right:

planck(wavelength, temp) {
Real norm = pow(555.0, 5.0) * (exp(h_c / (555.0 * k * temp)) - 1.0);
return norm / (pow(wavelength, 5.0) * (exp(h_c / (wavelength * k * temp)) - 1.0));
}


3_ Say I have an array of wavelengths from 360nm to 830nm, in increments of 5nm, and I apply the planck function with temperature of 6504k to each wavelength, and store the results in a new array. If I use this new array (spectral power distribution) and calculate the XYZ values for it, it should match CIE Standard Illuminant D65, should it not?

-
Pro tip: Either use a unit support system like Geant implements or always, always document the units of figures in the comments. You'll thank me later. – dmckee Jul 31 '12 at 3:35
This question can't be answered, because the different forms are answers to slightly different questions. Are you asking how much light energy there is in a given wavelength going in a certain direction, or summed over all the directions? Averaging over directions gives you geometric power laws in front. Also you need to know if you are looking for the distribution of energy per increment of wavelength, or per increment of frequency, the two distributions differ by a power-factor in front. You have to be specific about the physical problem you want to solve. – Ron Maimon Jul 31 '12 at 7:17
@dmckee I documented the units of the ones that I know. – Arlen Jul 31 '12 at 15:05
@RonMaimon I'm trying to implement blackbody radiation, which is emitted according to Planck's Law as far as I can tell. And it needs to be in terms of wavelength, not frequency. – Arlen Jul 31 '12 at 15:08
The resolution of Ron's question is you care about the most physically and biologically correct model for the light transport and the response of the human eye and optical interpretation system. All of which is very complicated. – dmckee Jul 31 '12 at 15:09

From the Wikipedia link you posted, it looks like CIE Standard Illuminant D65 does not exactly model a black body spectrum. Look at the animation in the second figure. The red line is a perfect black body spectrum and the black line is the D65 one. I believe this is because the purpose of D65 is to represent the spectrum of daylight on an average day at Earth's surface. Thus it's approximately a black body spectrum (because the Sun is, to a good approximation, a black body), but the D65 spectrum also includes the effect of absorption by various things in the atmosphere.

For this reason, you shouldn't expect your XYZ values based on a perfect black body spectrum to perfectly match those obtained from the D65 standard. I haven't checked your code in detail, but it looks plausible, so if you're getting approximately correct answers then it's probably right.

-
My XYZ values in part 3 of my question are way off. I know my implementation of Spectrum to XYZ conversion is correct because I've tested it against other data sets. Link to the formula for those interested: brucelindbloom.com/Eqn_Spect_to_XYZ.html – Arlen Jul 31 '12 at 22:11