import sys
import random
from math import exp
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from matplotlib import pyplot as plt
"""
概要:
プランクの法則に基づき、最大強度波長 (λm) から物体の絶対温度 (T) を計算し、
その結果を近似値と比較してグラフ表示するスクリプト。
詳細説明:
このスクリプトは、プランクの放射法則から導かれるウィーンの変位則を数値的に解き、
特定の最大強度波長における物体の温度を求めます。
また、近似式による温度と比較し、ニュートン法による厳密解との差を可視化します。
物理定数には、プランク定数、光速、ボルツマン定数が使用されます。
関連リンク:
:doc:`radiation_thermometer_usage`
"""
# 定数
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
kB = 1.380658e-23 # JK<sup>-1</sup>";
# Wavelength range
lmin = 100 # nm
lmax = 10000
lstep = 100
nl = int((lmax - lmin) / lstep) + 1
# 解を求める関数: func(lm, T) = 0
[ドキュメント]
def func(lm, T):
"""
放射温度計の温度計算に使用する非線形方程式を定義します。
詳細説明:
この関数は、プランクの放射法則から導出されるウィーンの変位則
`h*c/(kB*T*lm) = 5.0 * (1 - exp(-h*c/(kB*T*lm)))`
を変形した `func(lm, T) = 0` の形式で定義されます。
ニュートン法 (`scipy.optimize.newton`) を用いて、
与えられた最大強度波長 `lm` に対する温度 `T` の厳密解を求めるために使用されます。
:param lm: float: 最大強度波長 (nm)。関数内部でメートルに変換されます。
:param T: float: 温度 (K)。
:returns: float: 方程式 `func(lm, T) = 0` の評価結果。
"""
beta = 1.0 / kB / T
bhcl = beta * h * c / (lm * 1.0e-9)
return -5.0 + bhcl / (1.0 - exp(-bhcl))
[ドキュメント]
def main():
"""
放射温度計の計算と結果のプロットを行うメイン処理です。
詳細説明:
この関数は以下の処理を実行します:
1. 指定された波長範囲 (`lmin` から `lmax`、ステップ `lstep`) で、
最大強度波長 `lm` を変化させながらループします。
2. 各 `lm` に対して、近似式 `Tapprox = h * c / 4.97 / kB / (lm * 1.0e-9)`
を用いて近似温度 `Tapprox` を計算します。
3. `scipy.optimize.newton` 関数と `func` 関数を使用して、
厳密な温度 `Texact` をニュートン法で計算します。
この際、`Tapprox` を初期値として使用します。
4. 近似温度と厳密温度の差 `Tdif` を計算します。
5. 各計算結果 (`lm`, `Tapprox`, `Texact`, `Tdif`) をコンソールに出力します。
6. 計算された温度（近似値と厳密解）および誤差を2つのサブプロットとしてグラフ表示します。
7. グラフ表示後、ユーザーからの入力待ち状態となり、Enterキーが押されるとプログラムを終了します。
:returns: None
"""
lx = []
Tapp = []
Texa = []
Tdif = []
print("lm (nm)\tTapproximation (K)\tTexact (K)\tdifference (K)")
for i in range(nl):
# ピーク波長
lm = lmin + i * lstep
# 近似式による温度
Tapprox = h * c / 4.97 / kB / (lm * 1.0e-9)
# newton関数を呼び出す。変数として、解を求める変数しか与えられないので、
# 変数をTとするlambda関数を使って、ピーク波長 lm をfuncに渡す
# 初期値として、近似式の値を渡す
Texact = optimize.newton(lambda T: func(lm, T), Tapprox)
lx.append(lm)
Tapp.append(Tapprox)
Texa.append(Texact)
# 誤差を計算
Tdif.append(Texact - Tapprox)
print("%6.1f\t%8.3f\t%8.3f\t%8.3g" % (lm, Tapprox, Texact, Texact - Tapprox))
#=============================
# グラフの表示
#=============================
fig = plt.figure()
# 上段のグラフは、近似式で得た値と厳密解をプロット
ax1 = fig.add_subplot(2, 1, 1)
# 下段のグラフは、誤差をプロット
ax2 = fig.add_subplot(2, 1, 2)
plt.subplots_adjust(wspace = 0.2, hspace = 0.3)
ax1.plot(lx, Tapp, label = 'approximation')
ax1.plot(lx, Texa, label = 'exact', color = 'red')
ax1.set_title('Temperature from maximum intensity wavelength')
ax1.set_xlabel("lambda_max (nm)")
ax1.set_ylabel("T (K)")
ax1.set_xlim([lmin, lmax])
ax1.set_ylim([0.0, max(Tapp)])
ax1.legend()
ax2.plot(lx, Tdif, label = 'difference')
ax2.set_xlabel("lambda_max (nm)")
ax2.set_ylabel("dT (K)")
ax2.set_xlim([lmin, lmax])
ax2.set_ylim([0.0, max(Tdif)])
ax2.legend()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
if __name__ == '__main__':
main()