matplotlib でコンプトン散乱断面積の図を描いてみよう

こんにちは。 山師です。

今回は少し計算機寄りのことを書こうと思います。 具体的には、反応断面積を \(\theta\) 方向の極座標にプロットする方法についてです。

何をしたいか

目的

修士論文を書くとき、電子の反応断面積の図を作る必要がありました。 具体的には、コンプトン散乱の微分散乱断面積の図を作って、レビューパートに載せました。 どうせ載せるなら綺麗な図で、しかも先輩からのお下がりではない図を載せたいと思ったので、自分でコードを書くことにしました。

物理学的背景

場の量子論によれば、無偏光に対する静止した自由電子のコンプトン散乱断面積は、次に示される Klein-Nishina の公式で与えられます。 $$\frac{d\sigma}{d\Omega} = \frac{1}{2} {r_0}^2 \left( \frac{E}{E_0} \right)^2 \left( \frac{E}{E_0} + \frac{E_0}{E} – \sin^2 \theta \right)$$ ここで、\(r_0\) は電子の古典半径、 \(E_0\) および \(E\) はそれぞれ散乱前後の光子のエネルギーです。 \(\theta\) は光子の入射方向に対する散乱方向を表します。 \(E_0\) と \(E\) との間には次の関係があります。 $$E_0 – E = \frac{E_0}{m_e c^2} \left( 1 – \cos \theta \right)$$ ここで、\(m_e = 511 {\rm\ keV} / c^2\) は電子の静止質量です。

コード

以下に Python でのコードを示します。 numpy と matplotlib を使用しています。

import math
import numpy as np
import matplotlib.pyplot as plt

# energy of incoming photons
E_INIT = [10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0] # keV

def energy_ratio(e_init, theta) : 
	return 1.0 / (1.0 + (e_init / 511.0) * (1 - math.cos(theta)))

# r0 [m]
def KN_cross_diff(e_init, theta, r0 = 2.818e-15) :
	hoge = energy_ratio(e_init, theta)
	s = math.sin(theta)
	return 0.5 * r0 * r0 * hoge * hoge * (hoge + 1.0 / hoge - s * s)

def plot_KN_cross(e_init) : 
	for e in e_init :
		x = np.linspace(0, 360, 361)
		y = np.array([KN_cross_diff(e, math.pi * ang / 180.0) for ang in x])
		plt.polar(x * math.pi / 180.0, y, label='%.1f keV'%(e))
	plt.grid(b = True, color='silver')
	plt.xlabel('\u03B8 (deg)')
	plt.ylim([0.0, 9.0e-30])
	plt.legend(loc='lower right', borderaxespad=0, ncol=1, bbox_to_anchor=(1.33, 0.0))
	plt.show()

def main() : 
	plot_KN_cross(E_INIT)

# magic code
if __name__ == '__main__' : 
	main()

plt.polar を使うのがミソです。 蟹ミソ食べたい。

結果

python3 plot_KN.py

単位は \({\rm m^2 / Sr}\) で出力されます。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です