RuntaScience diary

気象系データを扱う学生 旅が好きです

【Science】黒体放射の式Planckの法則(Planck's law)をPythonで計算&描写

この記事をシェアする

f:id:RuntaScience:20210501134934p:plain 黒体放射の式Planckの法則(Planck's law)の計算と描写 今回は波長と周波数の関数を使用。

Planck's law

プランクの法則とは、

物質と熱平衡状態の放射が出すエネルギー密度(放射輝度)分布を表す法則

である(引用:天文学辞典)。

式は以下のように表される(放射輝度:radiance[W m-3 sr-1]or [W m-2sr-1 Hz-1 ])。

f:id:RuntaScience:20210501134550p:plain f:id:RuntaScience:20210501134605p:plain

B: Spectral radiance of Black body (Energy)
λ: wavelength[m], f: frequency[Hz]
h: Plank's constant($=6.62606896×10^{-34}$)
k: Boltzmann constant($=1.3806504×10^{-23}$)
c: Speed of light($=2.99792458×108$)

コード

定数の定義

c = 2.99792458e8 #c=speed of light[m/s]
h = 6.62606896e-34 #h=Planck constant[Js]
k = 1.3806504e-23 #k=Boltzmann constant[J/K] 
T_arr = np.array([7000, 6000, 5000, 4000, 3000, 255]) #temprature[K]
lam_arr = np.linspace(1e-9,1e-2,1000000) #wavelength[m]
f_arr = np.linspace(1e10, 1e16,1000000) #frequency[Hz]

c,h,kは以下のサイトの値を使用 https://japanknowledge.com/contents/common/butsuriteisu.html

波長の関数

関数の定義
def planck_lam(lam, T): 
    """Definicion of Planck's law
        return B(lambda, temprature)[radiance:W/m^3/sr]
        calculation=>a * b
    """
    a = 2.0 * h * (c**2) / (lam**5)
    b = 1 / (np.exp(h*c/(k*lam*T)) - 1.0)
    radiance = a * b
    return radiance
描写
fig = plt.figure(figsize=(10,5))
plt.rcParams["font.size"] = 16
ax = plt.subplot(111)
ax.set_title("Blackbody spectrum (Radiance)")

for T in T_arr:
    B = planck_lam(lam_arr, T)
    # x-axis:lambda[m]→[nm](*1e9) y-axis: [W/m^3/sr]→[kW/m^2/sr/µm](*1e-3)
    ax.plot(lam_arr*1e9, B*1e-3*1e-6, label="T={}[K]".format(T))

#軸設定
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(1e2, 1e4)
ax.set_ylim(1e1, 1e5)
ax.set_ylabel("Spectral radiance [kW $\mathrm{m^{-2} sr^{-1} \mu m^{-1}}$]")
ax.set_xlabel("Wavelength [$\mathrm{n m}$]")

ax.grid(which="both",color="gray",linewidth=0.2)
ax.legend()

f:id:RuntaScience:20210501132931p:plain

周波数の関数

関数の定義
def planck_f(f, T): 
    """Definicion of Planck's law
        return B(frequency, temprature)[radiance:W/m^2/sr/Hz]
        calculation=>a * b
    """
    a = 2.0 * h * (f**3) / (c**2)
    b = 1 / (np.exp(h*f/(k*T)) - 1.0)
    radiance = a * b
    return radiance
描写
plt.rcParams["font.size"] = 16
ax = plt.subplot(111)
ax.set_title("Blackbody spectrum (Radiance)")

for T in T_arr:
    B = planck_f(f_arr, T)
    # x-axis:frequency[Hz], y-axis: [W/m^2/sr/Hz]→[kW/m^2/sr/Hz](*1e-3)
    ax.plot(f_arr, B*1e-3, label="T={}[K]".format(T))

#軸設定
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(1e12, 4e15)
ax.set_ylim(1e-15, 1e-9)
ax.set_ylabel("Spectral radiance [kW $\mathrm{m^{-2} sr^{-1} Hz^{-1}}$]")
ax.set_xlabel("Frequency [$\mathrm{Hz}$]")

ax.grid(which="both",color="gray",linewidth=0.2)
ax.legend()

f:id:RuntaScience:20210501133011p:plain

参考文献

それでは🌏

プライバシーポリシー
お問い合わせ