黒体輻射のスペクトル
黒体輻射のスペクトル（分光放射輝度）をプロットするPythonコードです．計算方法：黒体輻射のスペクトル - Notes_JP
import numpy as np import matplotlib.pyplot as plt import matplotlib.ticker as ptick
def Planck(lam, T): """ 分光放射輝度を返す (W/sr/m3) lam: 波長 (m) T: 温度 (K) """ h = 6.62607015e-34 # Js c = 299792458 # m/s kB = 1.380649e-23 # J/K spectrum = 2*h*(c**2)/lam**5 inv_exp = np.exp(-h*c/(kB*T*lam)) # オーバーフロー対策 spectrum *= inv_exp / (1 - inv_exp) return spectrum def Planck_peak_lam(T): """ ピーク波長を返す (m) T: 温度 (K) """ h = 6.62607015e-34 # Js c = 299792458 # m/s kB = 1.380649e-23 # J/K peak_lam = h*c/(4.96511*kB*T) return peak_lam
いくつかプロットしてみます．
英語版Wikiの図（Planck's law - Wikipedia）：
# wiki lam = np.linspace(0.1, 3000, 5000, endpoint=True) * 1e-9 T_all = [3000, 4000, 5000] peak_lam = [Planck_peak_lam(T) for T in T_all] peak_spec = [Planck(l, T)/1e9/1e3 for (l, T) in zip(peak_lam, T_all)] num = 1 fig = plt.figure(figsize = (6.4*num, 4.8), tight_layout = True) axes = [fig.add_subplot(1, num, i + 1) for i in range(num)] for T in T_all: spec = Planck(lam=lam, T=T)/1e9/1e3 axes[0].plot(lam*1e9, spec, label=f'{T} K') axes[0].set_xlabel('Wavelength (nm)') axes[0].set_ylabel('Spectral radiance (kW/sr m$^2$ nm)') axes[0].legend() axes[0].scatter(np.array(peak_lam)*1e9, peak_spec) for ax in axes: ax.set_xlim(0, 3000) ax.set_ylim(0, 14) ax.legend() ax.grid() plt.show()
次のページの図：プランクの法則 | 天文学辞典
lam = np.linspace(0.1, 100, 5000, endpoint=True) * 1e-6 T_all = [250, 500, 1000, 2000, 4000, 8000] peak_lam = [Planck_peak_lam(T) for T in T_all] peak_spec = [Planck(l, T)/1e6 for (l, T) in zip(peak_lam, T_all)] num = 1 fig = plt.figure(figsize = (6.4*num, 4.8), tight_layout = True) axes = [fig.add_subplot(1, num, i + 1) for i in range(num)] for T in T_all: spec = Planck(lam=lam, T=T)/1e6 axes[0].plot(lam*1e6, spec, label=f'{T} K') axes[0].set_xlabel('Wavelength ($\mu$m)') axes[0].set_ylabel('Spectral radiance (W/sr m$^2$ $\mu$m)') axes[0].legend() axes[0].scatter(np.array(peak_lam)*1e6, peak_spec) for ax in axes: ax.set_xlim(0.1, 100) ax.set_ylim(1e-1, 1e9) ax.set_xscale('log') ax.set_yscale('log') ax.legend() ax.grid() plt.show()
ある波長範囲のエネルギーの割合
全波長に占めるエネルギーに対する，ある波長範囲のエネルギーの割合を計算する．Stefan-Boltzmannの法則（黒体輻射のスペクトル - Notes_JP）の計算から，$x_{i}=h\nu_{i}/k_{\mathrm{B}}T$として
\begin{aligned}
\frac{15}{\pi^{4}}
\int_{x_{1}}^{x_{2}}
\frac{x^{3}}{e^{x} - 1}
\,\mathrm{d}x
\end{aligned}
が，「周波数$\nu_{1}$から$\nu_{2}$の間にあるエネルギー/全周波数のエネルギー」である．\frac{15}{\pi^{4}}
\int_{x_{1}}^{x_{2}}
\frac{x^{3}}{e^{x} - 1}
\,\mathrm{d}x
\end{aligned}
あとは$\lambda_{2} = c/\nu_{1}$，$\lambda_{1} = c/\nu_{2}$と波長に変換すれば計算できる．
import scipy.integrate as integrate def Planck_pct(lam1, lam2, T): """ [lam1, lam2]の波長範囲で全体の何%を占めるか lam: m T: K """ h = 6.62607015e-34 # Js c = 299792458 # m/s kB = 1.380649e-23 # J/K # 波長と周波数は大小関係が逆転する x1 = h*(c/lam2)/(kB*T) if lam2 != np.inf else 0 x2 = h*(c/lam1)/(kB*T) if lam1 !=0 else np.inf res = integrate.quad(lambda x:(x**3)*np.exp(-x)/(1-np.exp(-x)), x1, x2) return res[0]*15/(np.pi**4), res[1] # 計算例 ## (1.0000000000000002, 2.6284701417309788e-09) print(Planck_pct(lam1=0, lam2=np.inf, T=5000)) ## (0.6337258719159103, 5.794351198717696e-11) print(Planck_pct(lam1=0, lam2=1000e-9, T=5000))
