プランクの放射式を python (matplotlib) で描く

python のライブラリ matplotlib で、関数の図を描く遊び。
ここでは試しに黒体放射からのスペクトルを例にしてみた。
欧米で主流の、python でほぼ全てやってしまう流へのステップのひとつ。

import pylab

def planck(x, kt):
    norm = 2.1e-4
    y = norm*x**3/(pylab.exp(x/kt)-1)
    return y

xmin = 1e-2
xmax = 1e+4
ymin = 1e-8
ymax = 1e+3

fig = pylab.plt.figure()
x = pylab.arange(xmin, xmax, 0.1)

y1 = planck(x, 0.58)
pylab.loglog(x, y1, 'r', lw=1.5)

y2 = planck(x, 5.8)
pylab.loglog(x, y2, 'g', lw=1.5)

y3 = planck(x, 58)
pylab.loglog(x, y3, 'b', lw=1.5)

pylab.xlim(xmin, xmax)
pylab.ylim(ymin, ymax)
pylab.xlabel("Energy $\epsilon$ (eV)", fontsize=12, fontname='serif')
pylab.ylabel("Specific Radiative Intensity \
(erg s$^{-1}$ cm$^{-2}$ Hz$^{-1}$ str$^{-1}$)", 
             fontsize=12, fontname='serif')
pylab.title("Planck's law", fontsize=12, fontname='serif')
pylab.legend(("0.58 eV", "5.8 eV", "58 eV"), "upper left")
fig.savefig("planck.png", format='png')

結果、生成される図は

黒体放射の Planck の法則のメモ*1。
せっかくだから、物理も復習しておこう。

温度 $T$ (K) の黒体での周波数 $\nu$ における、単位体積、単位周波数あたりのエネルギー密度 $u (\nu, T)$ は、

$\quad u(\nu, T)=\frac{8\pi \nu^2}{c^3}\frac{h\nu}{\exp(\frac{h\nu}{k_BT})-1} \quad \quad (\frac{\textrm{erg}}{\textrm{cm^3} \quad \textrm{Hz}})$

と表される。ここで、

$\quad \frac{8\pi \nu^2}{c^3}d\nu$

は、単位体積あたりで周波数 $\nu \sim \nu+d\nu$ にあるモード（独立な粒子/波）の数であり、残りの部分は $h\nu$ 単位で離散化した平均エネルギーに対応する。

全方位 $4\pi$ に等しく放射されているとき、単位時間あたり、単位面積から、単位立体角に放出される、単位周波数での放射輝度は

$\quad I(\nu, T)=\frac{2 \nu^2}{c^2}\frac{h\nu}{\exp(\frac{h\nu}{k_BT})-1} \quad \quad (\frac{\textrm{erg}}{\textrm{s} \quad \textrm{cm^2} \quad \textrm{Hz} \quad \textrm{str}})$

である。これを光子のエネルギー $\epsilon$ で書き表すと、

$\quad 2.1\times 10^{-4} \quad \left( \frac{\epsilon}{\textrm{1\quad eV}} \right)^3 \quad \frac{1}{\exp $\epsilon/k_BT$-1} \quad \quad \quad (\frac{\textrm{erg}}{\textrm{s} \quad \textrm{cm^2} \quad \textrm{Hz} \quad \textrm{str}})$

になる。

*1:はてなダイアリーで tex 記法を使ってみよう、という試し。はてなで tex を各のは現状ではそれほど綺麗ではない@2010/08/22。