Radiation thermometer:
Calculate temperature from maximum intensity wavelength lm

Download script from .\radiation_thermometer.py

import sys
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 #9200 # nm
lmax = 10000 #9600
lstep = 100
nl = int((lmax - lmin) / lstep) + 1

# figure configration
figsize = (8, 12)


# 解を求める関数: 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(figsize = figsize)

# 上段のグラフは、近似式で得た値と厳密解をプロット
    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()