【第19話】【はやくち解説】エバネッセント波ってなに？【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 波動関数空間分布の角度依存性
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mat
import matplotlib.animation as animation

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

#################################################################
## 物理定数
#################################################################
#プランク定数
h = 6.62607015 * 1.0E-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.1093837015 * 1.0E-31
#電子ボルト
eV = 1.602176634 * 1.0E-19
#虚数単位
I = 0.0 + 1.0j

######################################
#　物理系の設定
######################################
#電子のエネルギー
E = 2.0 * eV
#波数
k = np.sqrt(2.0 * me * E / hbar**2 )
#空間刻み間隔
nm = 1E-9
#壁の高さ
V1 = 0.0 * eV
V2 = 1.0 * 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 ))

#角振動数
omega = E/hbar
#空間刻み数
NX = 500
#空間幅
x_min = -2.0 * nm
x_max = 2.0 * nm



######################################
#　反射係数と透過係数を計算する関数定義
######################################
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, cos1, cos2):
    # 転送行列の定義
    M21 = np.zeros((2,2), dtype=np.complex)

    # 境界
    M21[0][0] = (1.0 + 0.0j + ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0
    M21[0][1] = (1.0 + 0.0j - ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0
    M21[1][0] = (1.0 + 0.0j - ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0
    M21[1][1] = (1.0 + 0.0j + ( k1 * cos1 ) / ( k2 * cos2 ) ) / 2.0

    #反射係数と透過係数の計算
    rc = ReflectionCoefficient(M21)
    tc = TransmissionCoefficient(M21)
    return rc, tc, M21

######################################
# 空間分布の計算
######################################
#実時間
t = 0
thetaN = 90
#x軸・y軸のメッシュ生成
xs = np.linspace(x_min, x_max, NX)
ys = np.linspace(x_min, x_max, NX)
X, Y = np.meshgrid(xs, ys)

#アニメーションフレーム作成用
def update( theta_n ):
    #実時間
    theta = theta_n * np.pi / 2.0 / thetaN
    sin1 = np.sin(theta) * k / k1
    sin2 = np.sin(theta) * k / k2
    cos1 = np.sqrt(np.complex128( 1.0 - sin1**2 ))
    cos2 = np.sqrt(np.complex128( 1.0 - sin2**2 ))

    #波数の指定
    k1x = k1 * cos1
    k1y = k1 * sin1
    k2x = k2 * cos2
    k2y = k2 * sin2
    rc, tc, M21 = calculateCoefficient( k1, k2, cos1, cos2 )

    Z = np.sin( np.complex128(X/nm) )
    #グラフ消去
    fig.clear()
    for ix in range(NX):
        for iy in range(NX):
            x = xs[ix]
            y = ys[iy]
            if( x < 0 ):
                Z[iy][ix] = np.exp( I * ( k1x * x + k1y * y ) - I * omega * t )
                Z[iy][ix] += rc * np.exp( I * ( - k1x * x + k1y * y ) - I * omega * t )
            elif( 0 <= x ):
                Z[iy][ix] = tc * np.exp( I * ( k2x * x + k2y * y ) - I * omega * t )

    Z = Z[:-1, :-1]
    Z = Z.real

    #両軸位置の設定
    axes = fig.add_axes(    [0.07, 0.05/aspect, 0.78, 0.78/aspect])
    #カラーバー位置の設定
    bar_axes = fig.add_axes([0.87, 0.05/aspect, 0.05, 0.78/aspect])
    #グラフ描画
    img = axes.pcolormesh( X/nm, Y/nm, Z, vmin=-1.0, vmax=1.0, cmap=cmap)
    #カラーバーの追加
    fig.colorbar(img, cax=bar_axes)
    #見出し
    axes.set_title( u"角度:" + r"$" + str(round(theta_n/thetaN*90,2)) + r"{}^{\circ}$" , fontname="Meiryo" )


#カラーマップの配色配列
color_list = []
color_list.append( [ 0,   "blue"  ] )
color_list.append( [ 0.5, "black" ] )
color_list.append( [ 1,   "red"   ] )
#カラーマップ用オブジェクトの生成
cmap = mat.colors.LinearSegmentedColormap.from_list('cmap', color_list)

#アニメーションの設定
ani = animation.FuncAnimation( fig, update, range(thetaN), 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()