import random
from math import exp
from scipy import optimize # newton関数はscipy.optimizeモジュールに入っている
from matplotlib import pyplot as plt
"""
Radiation thermometer:
Calculate temperature from maximum intensity wavelength lm
"""
# 定数
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
kB = 1.380658e-23 # JK-1";
# Wavelength range
lmin = 100 # nm
lmax = 10000
lstep = 100
nl = int((lmax - lmin) / lstep) + 1
# 解を求める関数: func(lm, T) = 0
def func(lm, T):
beta = 1.0 / kB / T
bhcl = beta * h * c / (lm * 1.0e-9)
return -5.0 + bhcl / (1.0 - exp(-bhcl))
def main():
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()
