【第18話】【はやくち解説】スネルの法則ってなに？【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
#時間間隔
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()