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).
Solution P8.4.3
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.
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