【第13話】箱型障壁に平面波を照射したときのアニメーション！【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

################
#　物理系の設定
################
#電子のエネルギー
E = 1.0 * eV
#波数
k1 = np.sqrt(2.0 * me * E / hbar**2 )
#角振動数
omega = E/hbar
#１周期
T = 2.0 * np.pi / omega
#時間間隔
#dt = 0.2 * 10**-16
dt = T / 100
#時間刻み数
NT = 100

#空間刻み間隔
nm = 1E-9
#壁の厚さ
d = 2.0 * nm
#壁の高さ
V1 = 0 * eV
V2 = 0.8 * eV
V3 = 0 * eV

#描画範囲
x_min = -4.0 * nm
x_max =  6.0 * nm
#描画区間数
NX = 500
#座標点配列の生成
x1 = np.linspace(x_min, 0, NX)
x2 = np.linspace(0, d,     NX)
x3 = np.linspace(d, x_max, NX)

################
#　反射係数と透過係数を計算する関数定義
################
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]

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 ))
#反射係数と透過係数の計算
# 転送行列の定義
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);

###############################
## 透過率と反射率の入射波のエネルギー依存性のグラフ描画
###############################
#領域２の係数
A2_p = M21[0][0] * 1.0 + M21[0][1] * rc
A2_m = M21[1][0] * 1.0 + M21[1][1] * rc

#アニメーション作成用
ims=[]
    
#各時刻における計算
for tn in range(NT):
    t = dt * tn

    #波動関数の計算
    #領域１
    phi_I =      np.exp(   I * k1 * x1 - I * omega * t )
    phi_R = rc * np.exp( - I * k1 * x1 - I * omega * t )
    phi1 = phi_I + phi_R
    #領域２
    phi2_p = A2_p * np.exp(   I * k2 * x2 - I * omega * t )
    phi2_m = A2_m * np.exp( - I * k2 * x2 - I * omega * t )
    phi2 = phi2_p + phi2_m 
    #領域３
    phi_T = tc * np.exp(   I * k3 * (x3 - d) - I * omega * t )

    #各コマを描画
    img  = plt.plot(x1/nm, phi1.real, colors[0], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x1/nm, phi1.imag, colors[1], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x1/nm, abs(phi1), colors[2], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x2/nm, phi2.real, colors[0], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x2/nm, phi2.imag, colors[1], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x2/nm, abs(phi2), colors[2], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x3/nm, phi_T.real, colors[0], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x3/nm, phi_T.imag, colors[1], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x3/nm, abs(phi_T), colors[2], linestyle='solid', linewidth = 5.0)
    ims.append( img )


#グラフの描画（波動関数）
plt.title( u"箱型障壁へ照射した平面波の波動関数（" + r"$ E = 1.0 [{\rm eV}], d = 2.0[{\rm nm}], V = 0.8[{\rm eV}] $" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=30)
plt.ylabel(r"$ \psi(x, t) $", fontsize=30)

#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.12, top = 0.95)

#罫線の描画
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.vlines([0], -0.05, 2.05, "black", linestyles='solid', linewidth = 10 )
plt.vlines([2], -0.05, 2.05, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0], -5,  0, "black", linestyles='solid', linewidth = 10 )
plt.hlines([2], 0,  2, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0], 2,  6, "black", linestyles='solid', linewidth = 10 )


#描画範囲を設定
plt.xlim([x_min/nm, x_max/nm])
plt.ylim([-2.0, 2.5])

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)

#アニメーションの保存
ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)

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