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.e-7, 1.e-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:')
print('m = {:.10f} m.K, c = {:.10e} m'.format(*p.coef[::-1]))
print('rms residual:', resid[0])
print('Estimate of b = {:.10f} m.K'.format(p.coef[1]))
print('Exact value of Wien displacement law constant: {:.10f} m.K'.format(
pc.value('Wien displacement law constant')))
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.0028977725 m.K, c = -7.7787808484e-13 m
rms residual: 3.74963676734e-23
Estimate of b = 0.0028977725 m.K
Exact value of Wien displacement law constant: 0.0028977685 m.K