投稿日: 投稿者: simu-kun

【第14話】箱型の障壁へ照射した電子パルスの運動【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 箱型障壁へ照射したガウス波束の時間発展
#################################################################
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

######################################
#　物理系の設定
######################################
#波束の中心エネルギー
E0 = 20.0 * eV
#重ね合わせの数
NK = 200
#ガウス分布の幅
sigma = 2 / ( 1.25 * 10**-9 )
#波数の間隔
dk = 10.0 * sigma / NK
#ガウス分布の中心波数
k0 = np.sqrt( 2.0 * me * E0 / hbar**2)

#計算時間の幅
ts = 0
te = 1000
#時間間隔
dt = 0.05 * 10**-16

#空間刻み間隔
nm = 1E-9
#壁の厚さ
d = 0.2 * nm
#壁の高さ
V1 = 0.0 * eV
V2 = 40.0 * eV
V3 = 0.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 ):
    
    #反射係数と透過係数の計算
    # 転送行列の定義
    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    

    return rc, tc, A2_p, A2_m

######################################
#　波動関数の計算
######################################
#ガウス波束の初期位置
x0 = -3.0 * nm
#描画範囲
x_min = -5.0 * nm
x_max =  5.0 * nm
#描画区間数
NX = 500
#座標点配列の生成
x1 = np.linspace(x_min, 0, NX)
x2 = np.linspace(0, d,     5)
x3 = np.linspace(d, x_max, NX)

#ガウス波束の値を計算する関数
def Psi( x, t, region ):
    #波動関数値の初期化
    psi = x * (0.0 + 0.0j)
    #各波数ごとの寄与を足し合わせる
    for kn in range(NK):
        #各波数を取得
        k = k0 + dk * (kn - NK/2)
        E = (hbar * k)**2 / (2.0 * me)
        #波数から各振動数を取得
        omega = hbar / (2.0 * me) * k**2
        #各領域の波数
        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, A2_p, A2_m = calculateCoefficient( k1, k2, k3 )
        #ガウス分布
        Ck = np.exp( - I * k1 * x0 ) * np.exp( -( (k - k0) / (2.0 * sigma) )**2 )
        #平面波を足し合わせる
        if(region == 1):
            #入射波と反射波と透過波
            psi_I = np.exp( I * k1 * x - I * omega * t )
            psi_R = rc * np.exp( - I * k1 * x - I * omega * t )
            psi += (psi_I + psi_R) * Ck
        elif(region == 2):
            psi2_p = A2_p * np.exp(   I * k2 * x - I * omega * t )
            psi2_m = A2_m * np.exp( - I * k2 * x - I * omega * t )
            psi += (psi2_p + psi2_m) * Ck
        elif(region == 3):
            psi_T = tc * np.exp( I * k3 * ( x - d ) - I * omega * t )
            psi += psi_T * Ck

    return psi

#基準となる振幅を取得
psi_abs = abs(Psi( x0, 0, 1 ))

#################################################################
## ガウス波束の時間発展アニメーションの作成
#################################################################
#アニメーション作成用
ims=[]
    
#各時刻における計算
for tn in range(ts, te + 1):
    t = dt * tn

    #波動関数の計算
    psi1 = Psi( x1, t, 1 ) / psi_abs
    psi2 = Psi( x2, t, 2 ) / psi_abs
    psi3 = Psi( x3, t, 3 ) / psi_abs

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


#グラフの描画（波動関数）
plt.title( u"箱型障壁へ照射したガウス波束の波動関数（" + r"$ E_0 = 20.0 [{\rm eV}], d = 0.2[{\rm nm}], V = 40.0[{\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.02, 2.02, "black", linestyles='solid', linewidth = 10 )
plt.vlines([0.2], -0.02, 2.02, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0], -5,  0, "black", linestyles='solid', linewidth = 10 )
plt.hlines([2], 0,  0.2, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0], 0.2,  5, "black", linestyles='solid', linewidth = 10 )


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

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

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

#################################################################
## ガウス分布の描画（横軸：エネルギー）
#################################################################
fig_gaussian = plt.figure(figsize=(15, 8))

#波数の計算範囲
k_min = k0 - dk * NK / 2.0
k_max = k0 + dk * NK / 2.0

E_min = ( hbar * k_min )**2 / (2.0 * me) / eV
E_max = ( hbar * k_max )**2 / (2.0 * me) / eV

#波数座標点配列の生成
k = np.linspace(k_min, k_max, NK)
E = ( hbar * k )**2 / (2.0 * me) / eV
#正規分布
c_n = np.exp( -((k - k0) / (2.0 * sigma))**2 )

#グラフの描画（固有関数）
plt.title( u"ガウス分布（正規分布）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$ E[{\rm eV}] $", fontsize=30)
plt.ylabel(r"$ c_n/c_0 $", fontsize=30)

#罫線の描画
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([E_min, E_max])
plt.ylim([0, 1.05])

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

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

コメントを残す

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