【第12話】箱型障壁に照射したときの反射率と透過率を計算！【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 透過率と反射率の入射波のエネルギー依存性のグラフ描画（E = 1.0[eV]）
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#虚数単位
I = 0.0 + 1.0j

######################################
#　物理定数
######################################
#プランク定数
h = 6.6260896 * 10**-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19

######################################
#　物理系の設定
######################################
#空間刻み間隔
nm = 1E-9
#壁の厚さ
d = 2.0 * nm
#壁の高さ
V1 = 0 * eV
V2 = 1.0 * eV
V3 = 0 * eV

######################################
#　反射係数と透過係数を計算する関数定義
######################################
def ReflectionCoefficient( MM ):
 return - MM[1][0] / MM[1][1] 
def TransmissionCoefficient(MM):
 return (MM[0][0] * MM[1][1] - MM[0][1] * MM[1][0] ) / MM[1][1]
def calculateCoefficient( k1, k2, k3, d):
 # 転送行列の定義
 M21 = np.zeros((2,2), dtype=np.complex)
 T2  = np.zeros((2,2), dtype=np.complex)
 M32 = np.zeros((2,2), dtype=np.complex)
 # 境界
 M21[0][0] = (1.0 + 0.0j + k1 / k2) / 2.0
 M21[0][1] = (1.0 + 0.0j - k1 / k2) / 2.0
 M21[1][0] = (1.0 + 0.0j - k1 / k2) / 2.0
 M21[1][1] = (1.0 + 0.0j + k1 / k2) / 2.0
 M32[0][0] = (1.0 + 0.0j + k2 / k3) / 2.0
 M32[0][1] = (1.0 + 0.0j - k2 / k3) / 2.0
 M32[1][0] = (1.0 + 0.0j - k2 / k3) / 2.0
 M32[1][1] = (1.0 + 0.0j + k2 / k3) / 2.0
 T2[0][0] = np.exp( I * k2 * d )
 T2[0][1] = 0
 T2[1][0] = 0
 T2[1][1] = np.exp( - I * k2 * d )
 #転送行列の計算
 M31 = M32 @ T2 @ M21;
 #反射係数と透過係数の計算
 rc = ReflectionCoefficient(M31);
 tc = TransmissionCoefficient(M31);
 return rc, tc

#################################################################
## 透過率と反射率の入射波のエネルギー依存性のグラフ描画
#################################################################
#エネルギーの計算範囲
E_min = 0.001
E_max = 3.0
#計算点数
NK = 500

#エネルギーの行列
Es = np.linspace(E_min, E_max, NK)

#反射率と透過率の配列
Rs = []
Ts = []
for E in Es:
 E = E * eV
 #波数の計算
 k1 = np.sqrt(np.complex128( 2.0 * me * (E - V1) / hbar**2 ))
 k2 = np.sqrt(np.complex128( 2.0 * me * (E - V2) / hbar**2 ))
 k3 = np.sqrt(np.complex128( 2.0 * me * (E - V3) / hbar**2 ))
 rc, tc = calculateCoefficient( k1, k2, k3, d)
 Rs.append( abs(rc)**2 )
 Ts.append( abs(tc)**2 )

#グラフの描画（固有関数）
plt.title( u"透過率と反射率の入射波のエネルギー依存性のグラフ（壁の厚さ：" + r"$2.0[{\rm nm}]$" + ", 壁の高さ：" + r"$1.0[{\rm eV}]$" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(u"入射波のエネルギー" + r"$[{\rm eV}]$"  ,  fontname="Yu Gothic", fontsize=20, fontweight=1000)
plt.ylabel(u"反射率と透過率", fontname="Yu Gothic", fontsize=20, fontweight=1000)

#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([0, 3])
plt.ylim([0, 1])

#波数分布グラフの描画
plt.plot(Es, Rs, linestyle='solid', linewidth = 5)
plt.plot(Es, Ts, linestyle='solid', linewidth = 5)

#グラフの表示
plt.show()