投稿日: 投稿者: simu-kun

【第26話】【結果は渋いです】２次元電子波束を斜めに薄膜へ（一応トンネル効果）

#################################################################
## 薄膜へ照射したガウス波束（トンネル効果）
#################################################################
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 = False
flag_animation = True

#図全体
aspect = 0.9
fig = plt.figure(figsize=(6, 6 * 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 = 6
#フォルダ指定
dir_path = "output30/"
#フォルダ生成
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 = 30.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 = 150
NY = 150
#壁の高さ
V1 = 0.0 * eV
V2 = 10.5 * eV
V3 = 0.0 * eV
#壁の幅
d = 1.0 * nm

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

#空間幅
x_min = -30.0 * nm
x_max =  30.0 * nm
y_min = x_min
y_max = x_max

######################################
#　反射係数と透過係数を計算する関数定義
######################################
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 = 0 #-20.0 * nm * np.cos(theta0)
y0 = 0 #-20.0 * nm * np.sin(theta0)

#ガウス波束の値を計算する関数
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 - I * k1y * y0) * 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] * 1000

	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()

コメントを残す

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