Wien's displacement law: where a black body's spectrum peaks¶

The wavelength at which a black body's spectral radiance per unit wavelength (lusine.radiation.planck) peaks, Wien's displacement law, lambda_max = b / T.

The library holds no b: it is h c / (k x), x = 4.965114231744276... the root of x = 5 (1 - e^-x) (the zero of Planck's law's derivative in lambda), and h, c and k are lusine.constants', exact, so b is derived here from them, 2.897771955e-3 m K (CODATA's "Wien wavelength displacement law constant", exact). x is a number of mathematics, x = 5 + W0(-5 e^-5), W0 the principal branch of Lambert's W. Source: Wien, Sitzungsberichte der Berliner Akademie (1893); Rybicki & Lightman, Radiative Processes in Astrophysics, 1.5 (eq. 1.56); CODATA 2018 (Tiesinga, Mohr, Newell & Taylor, Rev. Mod. Phys. 93, 025010 (2021)): b = (h c / k) / 4.965 114 231...; x as Lambert's W gives it, Corless, Gonnet, Hare, Jeffrey & Knuth, "On the Lambert W function", Adv. Comput. Math. 5, 329 (1996), which takes Wien's law as an application. Checked against a reference in numpy, which has its own tests.

Source

Wien (1893); Rybicki & Lightman, Radiative Processes in Astrophysics, 1.5 (eq. 1.56); CODATA 2018 (Tiesinga et al., Rev. Mod. Phys. 93, 025010): b = (h c / k) / 4.965 114 231...; Corless et al., On the Lambert W function, Adv. Comput. Math. 5, 329 (1996), for x.

Proof

reference, invariant.

wien_peak.trait

namespace lusine.radiation;

float<m> wien_peak(float<K> temperature)
{
  let x = 4.965114231744276;
  return lusine.constants.h * lusine.constants.c / (lusine.constants.k * x * temperature);
}

The wavelength at which Planck's radiance per unit wavelength peaks, λ_max = b/T. The library holds no Wien constant b ("a derivation and not a datum"), so the kernel derives it, b = hc/(kx), from lusine.constants' h, c and k and x = 4.965114231744276…, the root of x = 5(1 − e⁻ˣ), a number of mathematics written in the kernel with its source. Reading the constants makes it a binding's kernel, not a standalone function: it runs in a program, as this page runs it.

In [1]:

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

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.

In [2]:

print(trait_run.source('radiation/wien_peak.trait'))

print(trait_run.source('radiation/wien_peak.trait'))

// lusine.radiation.wien_peak: the wavelength at which a black body's spectral radiance per unit wavelength
// (lusine.radiation.planck) peaks, Wien's displacement law, lambda_max = b / T. The library holds no b: it is
// h c / (k x), x = 4.965114231744276... the root of x = 5 (1 - e^-x) (the zero of Planck's law's derivative in
// lambda), and h, c and k are lusine.constants', exact, so b is derived here from them, 2.897771955e-3 m K (CODATA's
// "Wien wavelength displacement law constant", exact). x is a number of mathematics, x = 5 + W0(-5 e^-5), W0 the
// principal branch of Lambert's W.
// Source: Wien, Sitzungsberichte der Berliner Akademie (1893); Rybicki & Lightman, Radiative Processes in
// Astrophysics, 1.5 (eq. 1.56); CODATA 2018 (Tiesinga, Mohr, Newell & Taylor, Rev. Mod. Phys. 93, 025010 (2021)):
// b = (h c / k) / 4.965 114 231...; x as Lambert's W gives it, Corless, Gonnet, Hare, Jeffrey & Knuth, "On the Lambert
// W function", Adv. Comput. Math. 5, 329 (1996), which takes Wien's law as an application.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.radiation;

float<m> wien_peak(float<K> temperature)
{
  let x = 4.965114231744276;
  return lusine.constants.h * lusine.constants.c / (lusine.constants.k * x * temperature);
}

The program¶

Each section runs the kernel over a lattice of items, one a case.

In [3]:

print(trait_run.law_program('radiation/wien_peak.trait', 'wien_peak'))

print(trait_run.law_program('radiation/wien_peak.trait', 'wien_peak'))

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host float<K> temperature;
  float<m> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.radiation.wien_peak(items.temperature); }

The root¶

x is where the derivative of Planck's law in λ is zero, (x − 5)eˣ + 5 = 0, whose positive root is 5 + W₀(−5e⁻⁵) (Lambert's W, scipy's): the kernel's digits are that root's, and it solves x = 5(1 − e⁻ˣ) to float64's rounding.

In [4]:

import re
from scipy.special import lambertw
x = float(re.search(r'let x = ([0-9.]+);', trait_run.source('radiation/wien_peak.trait')).group(1))
root = 5.0 + float(lambertw(-5.0 * np.exp(-5.0)).real)
assert abs(x - root) < 1e-15 and abs(x - 5 * (1 - np.exp(-x))) < 1e-15, (x, root)
print(x, root)

import re from scipy.special import lambertw x = float(re.search(r'let x = ([0-9.]+);', trait_run.source('radiation/wien_peak.trait')).group(1)) root = 5.0 + float(lambertw(-5.0 * np.exp(-5.0)).real) assert abs(x - root) < 1e-15 and abs(x - 5 * (1 - np.exp(-x))) < 1e-15, (x, root) print(x, root)

4.965114231744276 4.965114231744276

Reference: CODATA's b, and the Sun¶

At T = 1 K the kernel's λ_max is b itself: CODATA's Wien wavelength displacement law constant, 2.897 771 955 × 10⁻³ m K (exact since 2019, scipy's CODATA table), to binary32's rounding. The Sun's photosphere, 5778 K, peaks at 501.5 nm.

In [5]:

import scipy.constants as sc
codata = sc.physical_constants['Wien wavelength displacement law constant'][0]
b, sun = trait_run.law('radiation/wien_peak.trait', 'wien_peak', temperature=[1.0, 5778.0])
assert abs(b / codata - 1) <= 4 * EPS and abs(sun - 501.5e-9) < 0.05e-9, (b, sun)
print(f'b = {b:.9e} m K (CODATA {codata}); the Sun peaks at {sun * 1e9:.2f} nm')

import scipy.constants as sc codata = sc.physical_constants['Wien wavelength displacement law constant'][0] b, sun = trait_run.law('radiation/wien_peak.trait', 'wien_peak', temperature=[1.0, 5778.0]) assert abs(b / codata - 1) <= 4 * EPS and abs(sun - 501.5e-9) < 0.05e-9, (b, sun) print(f'b = {b:.9e} m K (CODATA {codata}); the Sun peaks at {sun * 1e9:.2f} nm')

b = 2.897771774e-03 m K; the Sun peaks at 501.52 nm

Invariant: where Planck's law peaks, at every temperature¶

For 300 K to 30 000 K, the emitted lusine.radiation.planck at 4001 wavelengths over ±1% of the kernel's λ_max, a cubic fitted through them (binary32's rounding makes the top too flat for the largest sample to say), peaks within 10⁻⁵ of it; λ_max T is the one constant b at every temperature, to a rounding.

In [6]:

ts = np.array([300.0, 1000.0, 3000.0, 5778.0, 10000.0, 30000.0])
peaks = trait_run.law('radiation/wien_peak.trait', 'wien_peak', temperature=ts)
for t, p in zip(ts, peaks):
    d = np.linspace(-0.01, 0.01, 4001)
    b = trait_run.law('radiation/planck.trait', 'planck', wavelength=p * (1 + d), temperature=t)
    c3, c2, c1, _ = np.polyfit(d, b / b.max(), 3)
    roots = np.roots([3 * c3, 2 * c2, c1]).real
    off = float(roots[np.argmin(np.abs(roots))])
    assert abs(off) < 1e-5, (t, off)
    print(f'{t:7g} K: lambda_max {p:.5e} m, Planck peaks at {off:+.1e} of it')
assert np.all(np.abs(peaks * ts / codata - 1) <= 4 * EPS)
fig, ax = plt.subplots(figsize=(7, 3.2))
t = np.geomspace(1000.0, 30000.0, 100)
ax.loglog(t, trait_run.law('radiation/wien_peak.trait', 'wien_peak', temperature=t) * 1e9)
ax.axhspan(380, 780, color='gold', alpha=0.3, label='visible')
ax.set_xlabel('T (K)')
ax.set_ylabel('lambda_max (nm)')
ax.legend()
plt.show()

ts = np.array([300.0, 1000.0, 3000.0, 5778.0, 10000.0, 30000.0]) peaks = trait_run.law('radiation/wien_peak.trait', 'wien_peak', temperature=ts) for t, p in zip(ts, peaks): d = np.linspace(-0.01, 0.01, 4001) b = trait_run.law('radiation/planck.trait', 'planck', wavelength=p * (1 + d), temperature=t) c3, c2, c1, _ = np.polyfit(d, b / b.max(), 3) roots = np.roots([3 * c3, 2 * c2, c1]).real off = float(roots[np.argmin(np.abs(roots))]) assert abs(off) < 1e-5, (t, off) print(f'{t:7g} K: lambda_max {p:.5e} m, Planck peaks at {off:+.1e} of it') assert np.all(np.abs(peaks * ts / codata - 1) <= 4 * EPS) fig, ax = plt.subplots(figsize=(7, 3.2)) t = np.geomspace(1000.0, 30000.0, 100) ax.loglog(t, trait_run.law('radiation/wien_peak.trait', 'wien_peak', temperature=t) * 1e9) ax.axhspan(380, 780, color='gold', alpha=0.3, label='visible') ax.set_xlabel('T (K)') ax.set_ylabel('lambda_max (nm)') ax.legend() plt.show()

    300 K: lambda_max 9.65924e-06 m, Planck peaks at -6.8e-08 of it
   1000 K: lambda_max 2.89777e-06 m, Planck peaks at -1.9e-07 of it
   3000 K: lambda_max 9.65924e-07 m, Planck peaks at +1.1e-07 of it
   5778 K: lambda_max 5.01518e-07 m, Planck peaks at -2.8e-07 of it
  10000 K: lambda_max 2.89777e-07 m, Planck peaks at +4.2e-08 of it
  30000 K: lambda_max 9.65924e-08 m, Planck peaks at -1.1e-08 of it