カテゴリー: 量子力学

#################################################################
## 有限の高さの障壁へ照射した平面波の時間発展（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 = 1000
#空間刻み間隔
dx = 1E-9
#描画範囲
x_min = -5.0 * dx
x_max =  5.0 * dx
#描画区間数
NX = 500
#座標点配列の生成
x1 = np.linspace(x_min, 0, NX)
x2 = np.linspace(0, x_max, NX)

#壁の高さ
V = 0.8 * eV
k2 = np.sqrt(np.complex128( 2.0 * me * (E - V) / hbar**2 ))

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

    #反射係数
    rc = ( k1 - k2 ) / ( k1 + k2 )
    #透過係数
    tc = 2.0 * k1 / ( k1 + k2 )

    #波動関数の計算
    phi_I = np.exp( I * k1 * x1 - I * omega * t )
    phi_R = rc * np.exp( - I * k1 * x1 - I * omega * t )
    phi_T = tc * np.exp( I * k2 * x2 - I * omega * t )
    phi = phi_I + phi_R

    #各コマを描画
    img  = plt.plot(x1/dx, phi_I.real, colors[0], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x1/dx, phi_R.real, colors[1], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x1/dx, phi.real, colors[2], linestyle='solid', linewidth = 5.0)
    img += plt.plot(x2/dx, phi_T.real, colors[3], linestyle='solid', linewidth = 5.0)
    ims.append( img )

#グラフの描画（波動関数）
plt.title( u"有限の高さの障壁へ照射した平面波の波動関数（" + r"$ E = 1.0 [{\rm eV}] $" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=30)
plt.ylabel(r"$ {\rm Re}\{\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, 1.05, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0], -5,  0, "black", linestyles='solid', linewidth = 10 )
plt.hlines([1], 0,  5, "black", linestyles='solid', linewidth = 10 )


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

#アニメーションの生成
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()

#################################################################
## 無限に高い障壁に向けてガウス波束を照射したときの時間発展（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

#################################################################
## 物理系の設定
#################################################################
#波束の中心エネルギー
E0 = 10.0 * eV
#重ね合わせの数
NK = 200
#ガウス分布の幅
sigma = 2 / ( 1.25 * 10**-9 )
#波数の間隔
dk = 10.0 * sigma / NK
#ガウス分布の中心波数
k0 = np.sqrt( 2.0 * me * E0 / hbar**2)
#空間分割サイズ
dx = 1.0 * 10**-9
#空間分割数
NX = 500
#描画範囲
x_min = -10.0 * dx
x_max =  0.0 * dx

#ガウス波束の初期位置
x0 = -7.5 * dx

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

#ガウス波束の値を計算する関数
def Psi( x, t ):
    #波動関数値の初期化
    psi = x * (0.0 + 0.0j)
    #各波数ごとの寄与を足し合わせる
    for kn in range(NK):
        #各波数を取得
        k = k0 + dk * (kn - NK/2)
        #波数から各振動数を取得
        omega = hbar / (2.0 * me) * k**2
        #入射波と反射波
        phi_I = np.exp( I * k * x - I * omega * t )
        phi_R = - np.exp( - I * k * x - I * omega * t )
        #平面波を足し合わせる
        psi += (phi_I + phi_R) * np.exp( - I * k * x0 ) * np.exp( -( (k - k0) / (2.0 * sigma) )**2 )

    return psi

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

#座標点配列の生成
x = np.linspace(x_min, x_max, NX)

#アニメーション作成用
ims=[]

### 波束の空間分布の計算
for tn in range(ts, te + 1):
    #実時間の
    t = dt * tn
    #時刻t=0の波動関数を取得
    psi = Psi( x, t )
    #波動関数の規格化
    psi /= psi_abs
    #各コマを描画
    img  = plt.plot(x/dx, psi.real, colors[0], linestyle='solid', linewidth = 3.0)
    img += plt.plot(x/dx, psi.imag, colors[1], linestyle='solid', linewidth = 3.0)
    img += plt.plot(x/dx, abs(psi),  colors[2], linestyle='solid', linewidth = 3.0)
    time = plt.text(8, 1.15, "t = " +  str(round(t/10**-15, 1)) + r"$[{\rm fs}]$" , ha='left', va='center', fontsize=18)
    #テキストをグラフに追加
    img.append(time)
    #アニメーションに追加
    ims.append(img)

#グラフの描画（波動関数）
plt.title( u"無限に高い障壁に向けて照射したガウス波束の波動関数（" + r"$ E = 10.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.15, 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], -5, 5, "black", linestyles='solid', linewidth = 10 )

#描画範囲を設定
plt.xlim([x_min/dx, x_max/dx + 1])
plt.ylim([-1.1, 1.5])

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=50)
#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
ani.save("output.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

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

#グラフの描画（固有関数）
plt.title( u"ガウス分布（正規分布）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$ (k_n - \bar{k})/\sigma $", 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([(k_min- k0)/sigma, (k_max- k0)/sigma])
plt.ylim([0, 1.05])

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

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

※計算時間がかなり掛かります。

#################################################################
## 薄膜へ照射したガウス波束（トンネル効果）
#################################################################
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mat
import matplotlib.animation as animation
from multiprocessing import Pool

flag_calculate = True
flag_animation = True

#図全体
aspect = 0.9
fig = plt.figure(figsize=(8, 4 * 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']

#並列数
ParallelNumber = 12
#フォルダ指定
dir_path = "output/"
#フォルダ生成
os.makedirs(dir_path, exist_ok = True)
#################################################################
## 物理定数
#################################################################
#プランク定数
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

######################################
#　物理系の設定
######################################
#波束の中心エネルギー
E0 = 10.0 * eV
#入射角度
theta0 = 0 / 180 * np.pi
#重ね合わせの数
NK = 100
#ガウス分布の幅
sigma = 0.1 / ( 1.25 * 10**-9 )
#波数の間隔
dk = 10.0 * sigma / NK
#ガウス分布の中心波数
k0 = np.sqrt( 2.0 * me * E0 / hbar**2)
kx0 = k0 * np.cos(theta0)
ky0 = k0 * np.sin(theta0)

#空間刻み間隔
nm = 1E-9
#空間刻み数
NX = 200
NY = 100
#壁の高さ
V1 = 0.0 * eV
V2 = 10.5 * eV
V3 = 0.0 * eV
#壁の幅
d = 1.0 * nm

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

#空間幅
x_min = -50.0 * nm
x_max =  50.0 * nm
y_min = x_min/2
y_max = x_max/2

######################################
#　反射係数と透過係数を計算する関数定義
######################################
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, cos1, cos2, cos3, 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 * 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
    M32[0][0] = (1.0 + 0.0j + ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
    M32[0][1] = (1.0 + 0.0j - ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
    M32[1][0] = (1.0 + 0.0j - ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
    M32[1][1] = (1.0 + 0.0j + ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
    T2[0][0] = np.exp( I * k2 * cos2 * d )
    T2[0][1] = 0
    T2[1][0] = 0
    T2[1][1] = np.exp( - I * k2 * cos2 * d )
    #転送行列の計算
    M31 = M32 @ T2 @ M21
    #反射係数と透過係数の計算
    rc = ReflectionCoefficient(M31)
    tc = TransmissionCoefficient(M31)
    return rc, tc, M21

#透過係数と反射係数と転送行列の配列
rcs = np.zeros((NK, NK))
tcs = np.zeros((NK, NK))
rcs = np.array(rcs, dtype=complex)
tcs = np.array(tcs, dtype=complex)
M21s = [ 0 ] * NK
for n1 in range(NK):
    M21s[ n1 ] = [ np.zeros((2,2), dtype=np.complex) ] * NK
for kxn in range(NK):
    for kyn in range(NK):
        kx = kx0 + dk * (kxn - NK/2)
        ky = ky0 + dk * (kyn - NK/2)
        k = np.sqrt( kx**2 + ky**2 )
        k1 = np.sqrt(np.complex128( k**2 - 2.0 * me * V1 / hbar**2 ))
        k2 = np.sqrt(np.complex128( k**2 - 2.0 * me * V2 / hbar**2 ))
        k3 = np.sqrt(np.complex128( k**2 - 2.0 * me * V3 / hbar**2 ))
        sin = ky / k
        cos = kx / k
        sin1 = sin * k / k1
        sin2 = sin * k / k2
        sin3 = sin * k / k3
        cos1 = np.sqrt(np.complex128( 1.0 - sin1**2 ))
        cos2 = np.sqrt(np.complex128( 1.0 - sin2**2 ))
        cos3 = np.sqrt(np.complex128( 1.0 - sin3**2 ))
        rcs[kxn][kyn], tcs[kxn][kyn], M21s[kxn][kyn] = calculateCoefficient( k1, k2, k3, cos1, cos2, cos3, d)

######################################
# 波動関数の計算
######################################
#ガウス波束の初期位置
x0 = -35.0 * nm

#ガウス波束の値を計算する関数
def Psi( x, y, t ):
    #波動関数値の初期化
    psi = 0 + 0j
    #各波数ごとの寄与を足し合わせる
    for kxn in range(NK):
        for kyn in range(NK):
            #各波数を取得
            kx = kx0 + dk * (kxn - NK/2)
            ky = ky0 + dk * (kyn - NK/2)
            k = np.sqrt( kx**2 + ky**2 )
            k1 = np.sqrt(np.complex128( k**2 - 2.0 * me * V1 / hbar**2 ))
            k2 = np.sqrt(np.complex128( k**2 - 2.0 * me * V2 / hbar**2 ))
            k3 = np.sqrt(np.complex128( k**2 - 2.0 * me * V3 / hbar**2 ))
            sin = ky / k
            cos = kx / k
            sin1 = sin * k / k1
            sin2 = sin * k / k2
            sin3 = sin * k / k3
            cos1 = np.sqrt(np.complex128( 1.0 - sin1**2 ))
            cos2 = np.sqrt(np.complex128( 1.0 - sin2**2 ))
            cos3 = np.sqrt(np.complex128( 1.0 - sin3**2 ))
            k1x = k1 * cos1
            k1y = k1 * sin1
            k2x = k2 * cos2
            k2y = k2 * sin2
            k3x = k3 * cos3
            k3y = k3 * sin3

            #波数から各振動数を取得
            omega = hbar / (2.0 * me) * ( kx**2  + ky**2 )

            W = np.exp( - I * k1x * x0 ) * np.exp( -( (kx - kx0) / (2.0 * sigma) )**2  -( (ky - ky0) / (2.0 * sigma) )**2 )
            if( x < 0 ):
                psi += np.exp( I * ( k1x * x + k1y * y ) - I * omega * t ) * W
                psi += rcs[kxn][kyn] * np.exp( I * ( - k1x * x + k1y * y ) - I * omega * t ) * W
            elif( x < d ):
                A2_p = M21s[kxn][kyn][0][0] + M21s[kxn][kyn][0][1] * rcs[kxn][kyn]
                A2_m = M21s[kxn][kyn][1][0] + M21s[kxn][kyn][1][1] * rcs[kxn][kyn]
                psi += A2_p * np.exp( I * (   k2x * x + k2y * y ) - I * omega * t ) * W
                psi += A2_m * np.exp( I * ( - k2x * x + k2y * y ) - I * omega * t ) * W
            elif( d <= x ):
                psi += tcs[kxn][kyn] * np.exp( I * ( k3x * (x - d) + k3y * y ) - I * omega * t ) * W

    return psi

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

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

#並列用ラッパー関数
def wrap_timeEvolution( data ):
    return timeEvolution( *data ) 

#時間発展の計算
def timeEvolution( filepath , tn ):
    #実時間
    t = ( ts + tn ) * dt
    Z = np.sin( np.complex128(X/nm) )
    for ix in range(NX):
        for iy in range(NY):
            x = xs[ix]
            y = ys[iy]
            Z[iy][ix] = Psi( x, y, t ) / psi_abs

    #ファイルへ保存
    np.savetxt( filepath, abs(Z))

#アニメーションフレーム作成用
def update( tn ):
    #グラフ消去
    fig.clear()
    #両軸位置の設定
    axes = fig.add_axes(    [0.07, 0.1/aspect, 0.78, 0.78/aspect])
    #カラーバー位置の設定
    bar_axes = fig.add_axes([0.87, 0.1/aspect, 0.05, 0.78/aspect])
    #データ読み込み    
    filepath = dir_path + str(tn + ts) + ".txt"
    Z = np.loadtxt( filepath )
    #透過領域の振幅を定数倍
    for ny in range(len(X)):
        for nx in range(len(X[ny])):
            if(X[ny][nx] > 1.0 * nm ):
                Z[ny][nx] = Z[ny][nx] * 5

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

    #２次元カラーグラフ描画
    img = axes.pcolormesh( X/nm, Y/nm, Z, vmin=0.0, vmax=1.0, cmap=cmap)
    #カラーバーの追加
    fig.colorbar(img, cax=bar_axes)



####### メイン関数として実行 ########
if __name__ == "__main__":

    #計算開始
    if(flag_calculate):
        list = []
        for tn in range( te - ts + 1 ):
            filepath = dir_path + str(tn + ts) + ".txt"
            list.append( (filepath, tn) )

        with Pool(processes=ParallelNumber) as p:
            p.map(wrap_timeEvolution, list)

    #アニメーション作成
    if( flag_animation ):
        #カラーマップの配色配列
        color_list = []
        color_list.append( [ 0,   "black" ] )
        color_list.append( [ 1,   "white" ] )

        #カラーマップ用オブジェクトの生成
        cmap = mat.colors.LinearSegmentedColormap.from_list('cmap', color_list)

        #アニメーションの設定
        ani = animation.FuncAnimation( fig, update, range(te - ts + 1), interval=10)
        #アニメーションの保存
        ani.save("output.html", writer=animation.HTMLWriter())
        #ani.save('anim.gif', writer='pillow')
        ani.save("output.mp4", writer="ffmpeg", dpi=300)

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

#################################################################
## 無限に深い量子井戸の波動関数
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(8, 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']

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

#################################################################
## 物理系の設定
#################################################################
#量子井戸の幅（m）
L = 1.0 * nm

#固有関数
def verphi(n, x):
    kn = np.pi * (n + 1) / L
    return np.sqrt(2.0 / L) * np.sin( kn * ( x + L / 2.0 ) )
#固有エネルギー（eV）
def Energy(n):
    kn = np.pi * (n + 1) / L
    return hbar**2 * kn**2 / (2.0 * me)
#波動関数
def phi(n, x, t):
    kn = np.pi * (n + 1) / L
    omega = hbar * kn**2 / (2.0 * me)
    return verphi(n,x) * np.exp(- I * omega * t)

#状態数
NS = 4
#固有エネルギーの表示（eV）
for n in range( NS + 1 ):
    if( n == 0 ): text = "基底状態"
    else: text = "第" + str(n) + "励起状態"
    print( text + "：" + str(Energy( n )/eV) + " [eV]" )

#計算時間の幅
ts = 0
te = 1

#基底状態の周期
T0 = 2.0 * np.pi * hbar / Energy(0)
print( "基底状態の周期：" + str(T0) + " [s]" )

#########################################################################
# 波動関数 アニメーション
#########################################################################
#時間間隔
dt = T0 / (te - ts + 1)

#描画範囲
x_min = -L/2
x_max = L/2
#描画区間数
NX = 500

#座標点配列の生成
x = np.linspace(x_min, x_max, NX)

plt.title( u"無限に深い量子井戸の波動関数", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=30)
plt.ylabel(r"$ {\rm Enegry}\,  [{\rm eV}]$", fontsize=30)

#アニメーション作成用
ims = []

#各時刻における波動関数の計算
for tn in range(ts, te):
    #実時間の取得
    t = dt * tn

    #各コマを描画
    for n in range( NS + 1 ):
        #波動関数の計算
        y = 0.5 * phi(n, x, t).real / (np.sqrt(2.0 / L)) + Energy(n) / eV
        #波動関数の表示
        if( n == 0 ): img  = plt.plot( x/nm, y, colors[n], linestyle='solid', linewidth = 5 )
        else:         img += plt.plot( x/nm, y, colors[n], linestyle='solid', linewidth = 5 )

    #アニメーションに追加
    ims.append( img )

#罫線の描画
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.6, 0.6])
plt.ylim([-0, 11])

#井戸の概形
plt.vlines([-0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.vlines([0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0],  -0.5, 0.5, "black", linestyles='solid', linewidth = 10 )

#E_0～E_4の表示
for n in range( NS + 1 ):
    #式の表示
    plt.text(0.62, Energy(n)/eV - 0.20, r"$E_" + str(n) + "$", fontsize = 30)
    #エネルギー準位の表示
    plt.hlines([Energy(n)/eV], -1, 2, colors[n], linestyles='dashed', linewidth = 2 )

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

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#アニメーションの保存
ani.save("output.html", writer=animation.HTMLWriter())
#ni.save('anim.gif', writer='pillow')
#ani.save("output.mp4", writer="ffmpeg", dpi=300)
#グラフの表示
#plt.show()

#################################################################
## 透過率と反射率の入射角依存性のグラフ描画（E = 2.0[eV]）
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mat
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 = 2.0 * eV
#波数
k = np.sqrt(2.0 * me * E / hbar**2 )
#空間刻み間隔
nm = 1E-9
#壁の高さ
V1 = 0.0 * eV
V2 = 1.0 * eV
V3 = 0.0 * eV
#壁の幅
d = 2.0 * nm
#波数の計算
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 ))

######################################
#　反射係数と透過係数を計算する関数定義
######################################
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, cos1, cos2, cos3, 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 * 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
    M32[0][0] = (1.0 + 0.0j + ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
    M32[0][1] = (1.0 + 0.0j - ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
    M32[1][0] = (1.0 + 0.0j - ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
    M32[1][1] = (1.0 + 0.0j + ( k2 * cos2 ) / ( k3 * cos3 ) ) / 2.0
    T2[0][0] = np.exp( I * k2 * cos2 * d )
    T2[0][1] = 0
    T2[1][0] = 0
    T2[1][1] = np.exp( - I * k2 * cos2 * d )
    #転送行列の計算
    M31 = M32 @ T2 @ M21
    #反射係数と透過係数の計算
    rc = ReflectionCoefficient(M31)
    tc = TransmissionCoefficient(M31)
    return rc, tc

#################################################################
## 透過率と反射率の壁の厚さ依存性のグラフ描画
#################################################################

#角度の計算範囲
theta_min = 0
theta_max = 89
#計算点数
NK = 891

#壁の厚さの行列
thetas = np.linspace(theta_min, theta_max, NK)
#
Rs = []
Ts = []
for theta in thetas:
    theta = theta / 180 * np.pi

    sin1 = np.sin(theta) * k / k1
    sin2 = np.sin(theta) * k / k2
    sin3 = np.sin(theta) * k / k3
    cos1 = np.sqrt(np.complex128( 1.0 - sin1**2 ))
    cos2 = np.sqrt(np.complex128( 1.0 - sin2**2 ))
    cos3 = np.sqrt(np.complex128( 1.0 - sin3**2 ))

    rc, tc = calculateCoefficient( k1, k2, k3, cos1, cos2, cos3, d )
    print( round( theta / np.pi * 180, 2 ) , abs(tc)**2 )
    Rs.append( abs(rc)**2 )
    Ts.append( abs(tc)**2 )

#グラフの描画（固有関数）
plt.title( u"透過率と反射率の入射角依存性（" + r"$E=2.0[{\rm eV}]$" + ", " + r"$V=1.0[{\rm eV}]$" + ", " + r"$d=2.0[{\rm nm}]$" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(u"入射角" + r"$[{}^{\circ}]$"  ,  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, 90])
plt.ylim([0, 1])

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

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

#################################################################
## 波動関数空間分布の角度依存性
#################################################################
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()

#################################################################
## スネルの法則
#################################################################
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
#時間間隔
dt = 10**-16
#空間刻み数
NX = 500
#１周期
T = 2.0 * np.pi / omega
#時間間隔
#dt = 0.2 * 10**-16
dt = T / 100
#時間刻み数
NT = 100
#空間幅
x_min = -2.0 * nm
x_max = 2.0 * nm
#角度
#theta = 0
#theta = 30.0 / 180 * np.pi
#theta = 45.0 / 180 * np.pi + 0.00001
theta = 50.0 / 180 * np.pi

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

######################################
#　反射係数と透過係数を計算する関数定義
######################################
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
#x軸・y軸のメッシュ生成
xs = np.linspace(x_min, x_max, NX)
ys = np.linspace(x_min, x_max, NX)
X, Y = np.meshgrid(xs, ys)
rc, tc, M21 = calculateCoefficient( k1, k2, cos1, cos2 )


#アニメーションフレーム作成用
def update( tn ):
	#実時間
	t = tn * dt
	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(tn) + r"[{\rm fs}]$" , 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(NT), 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()

#################################################################
## ２次元自由粒子の運動（E = 1.0[eV]）
#################################################################
import numpy as np
import matplotlib as mat
import matplotlib.pyplot as plt
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 = 1.0 * eV
#波数
k = 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
#空間刻み数
NX = 500
#空間刻み間隔
nm = 10**-9
#空間幅
x_min = -2.0
x_max = 2.0
#入射角度
theta = np.pi/6

#アニメーションフレーム作成用
def update( tn ):
    #グラフ消去
    fig.clear()
    #両軸位置の設定
    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])

    #実時間
    t = tn * dt
    #描画間隔
    delta_x = (x_max - x_min) / NX
  #x軸・y軸のメッシュ生成
    x = np.arange(x_min, x_max + delta_x, delta_x) * nm
    y = np.arange(x_min, x_max + delta_x, delta_x) * nm
    X, Y = np.meshgrid(x, y)
    #波数の指定
    kx = k * np.cos(theta)
    ky = k * np.sin(theta)
    #平面波の実部
    Z = np.exp( I * ( kx * X + ky * Y ) - I * omega * t ).real
    Z = Z[:-1, :-1]

    #グラフ描画
    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(tn/10) + r"[{\rm fs}]$" , fontname="Meiryo" )
    #軸ラベル
    #axes.set_xlabel("x")
    #axes.set_ylabel("y")

#カラーマップの配色配列
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(NT), 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()

#################################################################
## クローニッヒ・ペニーモデルのエネルギーバンド
#################################################################
import numpy as np
import numpy.linalg as LA
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
#井戸の形状
d1 = 1.5 * nm
d2 = 1.0 * nm
#周期
l = d1 + d2
#井戸の深さ
V1 = 0.0 * eV
V2 = -10.0 * eV

######################################
#　転送行列の計算
######################################
def TransferMatrix( k ):
    #各領域の波数
    k1 = np.sqrt(np.complex128( k**2 - 2.0 * me * V1 / hbar**2 ))
    k2 = np.sqrt(np.complex128( k**2 - 2.0 * me * V2 / hbar**2 ))
    # 転送行列の定義
    T1  = np.zeros((2,2), dtype=np.complex)
    T2  = np.zeros((2,2), dtype=np.complex)
    M21 = np.zeros((2,2), dtype=np.complex)
    M12 = 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
    M12[0][0] = (1.0 + 0.0j + k2 / k1) / 2.0
    M12[0][1] = (1.0 + 0.0j - k2 / k1) / 2.0
    M12[1][0] = (1.0 + 0.0j - k2 / k1) / 2.0
    M12[1][1] = (1.0 + 0.0j + k2 / k1) / 2.0
    T1[0][0] = np.exp( I * k1 * d1 )
    T1[0][1] = 0
    T1[1][0] = 0
    T1[1][1] = np.exp( - I * k1 * d1 )
    T2[0][0] = np.exp( I * k2 * d2 )
    T2[0][1] = 0
    T2[1][0] = 0
    T2[1][1] = np.exp( - I * k2 * d2 )
    #転送行列の計算
    M21T1 = M21 @ T1 
    M12T2 = M12 @ T2 
    M = M12T2 @ M21T1 

    #転送行列
    return M , M21T1, M12T2


#################################################################
## 転送行列の固有値
#################################################################
#バンド端エネルギー
BandEdgeEs = [0.158706, 0.266136, 0.615037, 0.849982]
#バンド中心エネルギー
BandCenterEs = [0.205879, 0.7251]
#エネルギー
E = ( BandCenterEs[1] )* eV 
#E = ( BandEdgeEs[3] )* eV 

k = np.sqrt( 2.0 * me * E / hbar ** 2)
M, M21T1, M12T2 = TransferMatrix( k )
print("------------------------")
eig, A = LA.eig(M)
print( eig )
print( A )
print("------------------------")
print( abs( M @ A -  A * eig ) )

LN = 5
LNs = np.linspace(0, LN, LN + 1)

Ap = [ 0 + 0j ] * ( 2 * LN + 1 )
Am = [ 0 + 0j ] * ( 2 * LN + 1 )



#################################################################
## 波動関数の描画
#################################################################

#描画範囲
x_min = 0.0 * nm
x_max = l * LN
#角振動数
omega = E/hbar
#１周期
T = 2.0 * np.pi / omega
#時間間隔
#dt = 0.2 * 10**-16
dt = T / 100
#時間刻み数
NT = 100

#アニメーション作成用
ims=[]

#描画区間数
NX = 100

xs = [0] * LN * 2
xK =np.linspace( x_min,  x_max, 500 )

#座標点配列の生成
for n in range(LN):
    xs[ 2 * n ] = np.linspace( l * n,  l * n + d1, NX )
    xs[ 2 * n + 1 ] = np.linspace( l * n + d1,  l * n + l, NX )

Ap[0] = A[0][0]
Am[0] = A[1][0]
for n in range(LN):
    Ap[2*n+1] = M21T1[0][0] * Ap[2*n] + M21T1[0][1] * Am[2*n]
    Am[2*n+1] = M21T1[1][0] * Ap[2*n] + M21T1[1][1] * Am[2*n]
    Ap[2*n+2] = M12T2[0][0] * Ap[2*n+1] + M12T2[0][1] * Am[2*n+1]
    Am[2*n+2] = M12T2[1][0] * Ap[2*n+1] + M12T2[1][1] * Am[2*n+1]

k1 = np.sqrt(np.complex128( k**2 - 2.0 * me * V1 / hbar**2 ))
k2 = np.sqrt(np.complex128( k**2 - 2.0 * me * V2 / hbar**2 ))

phi_p = [ 0 + 0j ] * ( 2 * LN )
phi_m = [ 0 + 0j ] * ( 2 * LN )
phi = [ 0 + 0j ] * ( 2 * LN )

#各時刻における計算
for tn in range(NT):
    t = dt * tn

    #波動関数の計
    for n in range(LN):
        #領域１
        phi_p[2*n] = Ap[2*n] * np.exp(   I * k1 * ( xs[2*n] - l * n ) - I * omega * t )
        phi_m[2*n] = Am[2*n] * np.exp( - I * k1 * ( xs[2*n] - l * n ) - I * omega * t )
        phi[2*n] = phi_p[2*n] + phi_m[2*n]
        #領域２
        phi_p[2*n+1] = Ap[2*n+1] * np.exp(   I * k2 * ( xs[2*n+1] - ( l * n + d1 ) ) - I * omega * t )
        phi_m[2*n+1] = Am[2*n+1] * np.exp( - I * k2 * ( xs[2*n+1] - ( l * n + d1 ) ) - I * omega * t )
        phi[2*n+1] = phi_p[2*n+1] + phi_m[2*n+1]

        if( n == 0 ):
            img  = plt.plot(xs[2*n]/nm, phi[2*n].real, colors[0], linestyle='solid', linewidth = 5.0)
        else:
            img  += plt.plot(xs[2*n]/nm, phi[2*n].real, colors[0], linestyle='solid', linewidth = 5.0)

        img += plt.plot(xs[2*n]/nm, phi[2*n].imag, colors[1], linestyle='solid', linewidth = 5.0)
        img += plt.plot(xs[2*n+1]/nm, phi[2*n+1].real, colors[0], linestyle='solid', linewidth = 5.0)
        img += plt.plot(xs[2*n+1]/nm, phi[2*n+1].imag, colors[1], linestyle='solid', linewidth = 5.0)


    ims.append( img )


#グラフの描画（波動関数）
plt.title( u"クローニッヒ・ペニーモデルの波動関数（" + r"$l=2.5[{\rm nm}]$" + " , " + r"$d=1.0[{\rm nm}]$"+ " , " + r"$V=-10.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)


#境界の描画
for n in range(LN):
    plt.hlines([0], n*l/nm,  (n*l+d1)/nm, "black", linestyles='solid', linewidth = 5 )
    plt.vlines([(n*l+d1)/nm], -1, 0, "black", linestyles='solid', linewidth = 5 )
    plt.hlines([-1], (n*l+d1)/nm,  (n*l + l)/nm, "black", linestyles='solid', linewidth = 5 )
    plt.vlines([(n*l + l)/nm], -1, 0, "black", linestyles='solid', linewidth = 5 )



#描画範囲を設定
plt.xlim([x_min/nm, x_max/nm])
plt.ylim([-1.5, 1.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()