The Wien displacement law · scipython.com

P8.4.3The Wien displacement law

The Wien displacement law predicts that the wavelength of maximum emission from a black body described by Planck's law is proportional to $$1/T$$ :

\lambda_\mathrm{max} T = b,

where $$b$$ is a constant known as Wien's displacement constant. Given the Planck distribution of emitted energy density as a function of wavelength,

u(\lambda, T) = \frac{8\pi^2hc}{\lambda^5}\frac{1}{e^{hc/\lambda k_\mathrm{B}T} - 1},

determine the constant $$b$$ by using scipy.optimize.minimize_scalar to find the maximum in $u(\lambda, T)$ for temperatures in the range $500\;\mathrm{K} \le T \le 6000\;\mathrm{K}$ and fitting $\lambda_\mathrm{max}$ to a straight line against $$1/T$$ . Compare with the "exact" value of $$b$$ , which is available within scipy.constants (see Section 8.1.1 of the book).

Show solution

The following code fits a striaght line (first-order polynomial) to $\lambda_\mathrm{max} = b/T$ obtained by finding the maximum in $u(\lambda, T)$ (minimum in $-u(\lambda, T)$ ) for a range of temperatures.

wiendisplacementlaw

import numpy as np
import scipy.constants as pc
from scipy.optimize import minimize_scalar

h = pc.value("Planck constant")
c = pc.value("speed of light in vacuum")
k = pc.value("Boltzmann constant")
c2 = pc.value("second radiation constant")

# The Planck function as a function of wavelength, u(lambda, T)
planck = lambda wv, T: 8 * np.pi * h * c / wv**5 / (np.exp(c2 / wv / T) - 1)


def lambda_max(T):
    """
    Return the maximum wavelength of the Planck function. This function
    brackets the maximum between 100 nm and 10 um so is only valid
    for 300 K <= T <= 28000 K.

    """
    res = minimize_scalar(
        lambda wv, T: -planck(wv, T), args=(T,), bracket=(1.0e-7, 1.0e-5)
    )
    return res["x"]


Tmin, Tmax, n = 500, 6000, 12
Tgrid = np.linspace(Tmin, Tmax, n)
lam_max = np.zeros(n)
for i, T in enumerate(Tgrid):
    lam_max[i] = lambda_max(T)

p, stats = np.polynomial.Polynomial.fit(
    1 / Tgrid, lam_max, 1, window=(Tmin, Tmax), domain=(Tmin, Tmax), full=True
)
resid, _, _, _ = stats
print("Fitted straight line m/T + c with:")
m, c = p.coef[::-1]
print(f"m = {m:.10f} m.K, c = {c:.10e} m")
print("rms residual:", resid[0])
b = p.coef[1]
print(f"Estimate of b = {b:.10f} m.K")
b_ref = pc.value("Wien wavelength displacement law constant")
print(f"Reference value of Wien displacement law constant: {b_ref:.10f} m.K")

Note that we need to bracket the minimum in order for Brent's method to converge correctly. The output from this code is:

Fitted straight line m/T + c with:
m = 0.0028977723 m.K, c = -7.7785266041e-13 m
rms residual: 3.7496492463185676e-23
Estimate of b = 0.0028977723 m.K
Reference value of Wien displacement law constant: 0.0028977720 m.K