Planck's law: a black body's spectral radiance¶
A black body's spectral radiance per unit wavelength, Planck's law, B(lambda, T) = (2 h c^2 / lambda^5) / (exp(h c / (lambda k T)) - 1), in W per square metre, per steradian, per metre of wavelength; h, c and k are lusine.constants', exact.
It is computed as 2 c E / lambda^4 / (exp(E / (k T)) - 1) with E = h c / lambda, the photon's energy, so no factor of it falls below binary32's smallest normal at the wavelengths light is measured at. Integrated over every wavelength it is sigma T^4 / pi (lusine.radiation.radiated_power over pi); its peak is at lusine.radiation.wien_peak. Far into the Rayleigh-Jeans tail (x = h c / (lambda k T) small) exp(x) - 1 keeps binary32's rounding over x, about 6e-6 at x = 0.02. Source: Planck, Annalen der Physik 4, 553 (1901); Rybicki & Lightman, Radiative Processes in Astrophysics, 1.5 (eq. 1.52). Checked against a reference in numpy, which has its own tests.
- Source
- Planck, Annalen der Physik 4, 553 (1901); Rybicki & Lightman, Radiative Processes in Astrophysics, 1.5 (eq. 1.52).
- Proof
- reference, invariant, limit.
planck.trait
namespace lusine.radiation;
float<W/(m3 sr)> planck(float<m> wavelength, float<K> temperature)
{
let photon = lusine.constants.h * lusine.constants.c / wavelength;
let x = photon / (lusine.constants.k * temperature);
return 2.0 * lusine.constants.c * photon / (wavelength * wavelength * wavelength * wavelength) / (exp(x) - 1.0);
}
A black body's spectral radiance per unit wavelength, B(λ, T) = (2hc²/λ⁵)/(exp(hc/(λkT)) − 1), in W/(m² sr m), with h, c and k read from lusine.constants. Reading the constants makes it a binding's kernel, not a standalone function: it runs in a program, as this page runs it.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.library_sources import parse_rows
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
import re
def constant(name):
"""A constant of lusine.constants as the library prints it."""
return float(re.search(r' ' + name + r' = ([-0-9.e]+) ', trait_run.source('constants.trait')).group(1))
The kernel¶
The law as the library states it.
print(trait_run.source('radiation/planck.trait'))
// lusine.radiation.planck: a black body's spectral radiance per unit wavelength, Planck's law,
// B(lambda, T) = (2 h c^2 / lambda^5) / (exp(h c / (lambda k T)) - 1),
// in W per square metre, per steradian, per metre of wavelength; h, c and k are lusine.constants', exact. It is
// computed as 2 c E / lambda^4 / (exp(E / (k T)) - 1) with E = h c / lambda, the photon's energy, so no factor of it
// falls below binary32's smallest normal at the wavelengths light is measured at. Integrated over every wavelength it
// is sigma T^4 / pi (lusine.radiation.radiated_power over pi); its peak is at lusine.radiation.wien_peak. Far into
// the Rayleigh-Jeans tail (x = h c / (lambda k T) small) exp(x) - 1 keeps binary32's rounding over x, about 6e-6 at
// x = 0.02.
// Source: Planck, Annalen der Physik 4, 553 (1901); Rybicki & Lightman, Radiative Processes in Astrophysics, 1.5
// (eq. 1.52).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.radiation;
float<W/(m3 sr)> planck(float<m> wavelength, float<K> temperature)
{
let photon = lusine.constants.h * lusine.constants.c / wavelength;
let x = photon / (lusine.constants.k * temperature);
return 2.0 * lusine.constants.c * photon / (wavelength * wavelength * wavelength * wavelength) / (exp(x) - 1.0);
}
The program¶
Each section runs the kernel over a lattice of items, one a case.
print(trait_run.law_program('radiation/planck.trait', 'planck'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float<m> wavelength;
host float<K> temperature;
float<W/(m3 sr)> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.radiation.planck(items.wavelength, items.temperature); }
Reference: the closed form, across the spectrum¶
At 300, 1000, 5778 and 10 000 K, at 4001 wavelengths over 50 nm to 1 cm, against Planck's law in float64 with lusine.constants' h, c and k: within 4ε(1 + x + 1/x), x = hc/(λkT): a few parts in 10⁶ about the peak, the kernel's exp(x) − 1 losing binary32's rounding over x far into the Rayleigh–Jeans tail, as its file says, and exp's argument's rounding times x far into the Wien tail. Where exp overflows binary32 (x above 88.7, far into the Wien tail, where the radiance is below 10⁻²⁰ of its peak) the kernel gives 0, never NaN.
h, c, k = constant('h'), constant('c'), constant('k')
lam = np.float32(np.geomspace(50e-9, 1e-2, 4001)).astype(float)
for t in (300.0, 1000.0, 5778.0, 10000.0):
b = trait_run.law('radiation/planck.trait', 'planck', wavelength=lam, temperature=t)
x = h * c / (lam * k * t)
with np.errstate(over='ignore'):
want = 2 * h * c * c / lam ** 5 / np.expm1(x)
assert np.all(np.isfinite(b)) and np.all(b[x > 88.8] == 0.0)
m = x < 88.0
rel = np.abs(b[m] / want[m] - 1)
assert np.all(rel <= 4 * EPS * (1 + x[m] + 1 / x[m])), (t, rel.max())
near = (x[m] > 1) & (x[m] < 10)
print(f'{t:7g} K: worst {rel[near].max():.1e} for 1 < x < 10, {rel.max():.1e} anywhere')
300 K: worst 8.5e-07 for 1 < x < 10, 1.2e-05 anywhere 1000 K: worst 1.7e-06 for 1 < x < 10, 3.8e-05 anywhere 5778 K: worst 8.5e-07 for 1 < x < 10, 2.3e-04 anywhere 10000 K: worst 1.6e-06 for 1 < x < 10, 3.5e-04 anywhere
Invariant: integrated over every wavelength it is σT⁴/π¶
The kernel's radiance at 4001 wavelengths from 50 nm to 1 cm, spaced evenly in log λ, integrated by the trapezoid rule over ln λ (B λ d ln λ), against lusine.constants.sigma T⁴/π: within 10⁻⁴, the integral's accuracy at this spacing and range.
sigma = constant('sigma')
lam = np.geomspace(50e-9, 1e-2, 4001)
u = np.log(lam)
fig, ax = plt.subplots(figsize=(7, 3.4))
for t in (1000.0, 3000.0, 5778.0):
b = trait_run.law('radiation/planck.trait', 'planck', wavelength=lam, temperature=t)
y = b * lam
total = float(np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(u)))
want = sigma * t ** 4 / np.pi
assert abs(total / want - 1) < 1e-4, (t, total / want - 1)
print(f'{t:6g} K: integral / (sigma T^4 / pi) - 1 = {total / want - 1:+.1e}')
ax.loglog(lam * 1e9, b, label=f'{t:g} K')
ax.set_xlim(100, 1e5)
ax.set_ylim(1e3, 1e14)
ax.set_xlabel('wavelength (nm)')
ax.set_ylabel('B (W / (m2 sr m))')
ax.legend()
plt.show()
1000 K: integral / (sigma T^4 / pi) - 1 = +3.4e-07 3000 K: integral / (sigma T^4 / pi) - 1 = +3.2e-07 5778 K: integral / (sigma T^4 / pi) - 1 = +4.5e-08
Limit: the peak and the Rayleigh–Jeans tail¶
At 5778 K, the Sun's photosphere, the kernel's radiance at 4001 wavelengths within 1% of 501.5 nm, a cubic fitted through them (binary32's rounding makes the top too flat for the largest sample to say), peaks at 501.52 nm, Wien's b/T (lusine.radiation.wien_peak). Far into the infrared, at 1 mm and 300 K (x = 0.048), it is Rayleigh and Jeans's 2ckT/λ⁴ times x/(eˣ − 1), to 10⁻⁵.
d = np.linspace(-0.01, 0.01, 4001)
b = trait_run.law('radiation/planck.trait', 'planck', wavelength=501.5e-9 * (1 + d), temperature=5778.0)
c3, c2, c1, _ = np.polyfit(d, b / b.max(), 3)
roots = np.roots([3 * c3, 2 * c2, c1]).real
peak = 501.5e-9 * (1 + float(roots[np.argmin(np.abs(roots))]))
assert abs(peak - 501.518e-9) < 0.005e-9, peak
far = trait_run.law('radiation/planck.trait', 'planck', wavelength=1e-3, temperature=300.0)[0]
x = h * c / (1e-3 * k * 300.0)
rj = 2 * c * k * 300.0 / 1e-3 ** 4
assert abs(float(far) / (rj * x / np.expm1(x)) - 1) < 1e-5
print(f'peak at {peak * 1e9:.2f} nm; at 1 mm, {float(far) / rj:.4f} of Rayleigh-Jeans')
peak at 501.52 nm; at 1 mm, 0.9762 of Rayleigh-Jeans