【第32話】コヒーレント状態ってどんな状態？【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## [032]コヒーレント状態の時間発展
#################################################################
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']

#################################################################
## 物理定数
#################################################################
#プランク定数
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

#################################################################
## 規格化エルミート多項式
#################################################################
def NormalizedHermitePolynomial( n, x ):
	H0 = (1.0 / np.pi)**(1.0/4.0) * np.exp( - x**2 / 2.0 ) 
	H1 = (4.0 / np.pi)**(1.0/4.0) * np.exp( - x**2 / 2.0 )  * x
	if(n==0): return H0
	if(n==1): return H1
	for m in range(2, n+1):
		H2 = np.sqrt( 2.0 / m ) * x * H1 -  np.sqrt( (m - 1) / m )  * H0
		H0 = H1
		H1 = H2
	return H2

#################################################################
## 物理系の設定
#################################################################
#調和ポテンシャルの強さ
omega = 1.0 * 1.0E+15
A = np.sqrt( me * omega / hbar)
#固有関数
def verphi(n, x):
	barX = A * x
	return np.sqrt(A) * NormalizedHermitePolynomial( n, barX )
#固有エネルギー（eV）
def Energy(n):
	return (n + 1.0/2.0) * hbar * omega

#################################################################
## コヒーレント状態
#################################################################
NC = 400

def CoherentState(v, x, t):
	_psi = 0.0 + 0.0j
	for n in range(NC):
		_psi += ff(n, v) * verphi(n, x) * np.exp( - I * omega * (n + 1.0/2.0) * t )
	return _psi	* np.exp( -abs(v)**2 / 2.0)

def ff(n, v, ln = 1):
	if( n==0 ): return 1
	if( n==1 ): return v * ln
	return ff( n-1, v,  ln * v / np.sqrt(n) )


#########################################################################
# 波動関数 アニメーション
#########################################################################
#計算時間の幅
ts = 0
te = 1000

#調和振動子の周期
T0 = 2 * 2.0 * np.pi / omega

#時間間隔
dt = T0 / (te - ts + 1)

plt.title( u"コヒーレント状態の時間発展", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=20)
plt.ylabel(r"$ \psi_v(x,t) $", fontsize=30)

L = 9 * nm
#アニメーション作成用
ims = []

#描画範囲
x_min = -L/2
x_max = L/2
#描画区間数
NX = 500
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)

#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.2, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.2, linestyle = "-", linewidth = 1)

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

#調和ポテンシャルの概形
yE = 1.0 / 2.0 * me * omega**2 * x**2 /eV /30.0

ims = []
#各時刻における波動関数の計算
for tn in range(ts, te+1):
	#実時間の取得
	t = dt * tn
	#波動関数の計算
	v = 4.0
	vx1 = CoherentState(1, x, t)
	vx4 = CoherentState(4, x, t)
	vx7 = CoherentState(7, x, t)
	vx1_real = vx1.real / np.sqrt(A)
	vx1_imag = vx1.imag / np.sqrt(A)
	vx4_real = vx4.real / np.sqrt(A)
	vx4_imag = vx4.imag / np.sqrt(A)
	vx7_real = vx7.real / np.sqrt(A)
	vx7_imag = vx7.imag / np.sqrt(A)

	#波動関数の表示
	img  = plt.plot( x/nm, yE, linestyle='dotted', color="black", linewidth = 1 )
	img += plt.plot( x/nm, vx1_real, colors[0], linestyle='solid', linewidth = 3 )
	img += plt.plot( x/nm, vx1_imag, colors[1], linestyle='solid', linewidth = 3 )
	img += plt.plot( x/nm, vx4_real, colors[2], linestyle='solid', linewidth = 3 )
	img += plt.plot( x/nm, vx4_imag, colors[3], linestyle='solid', linewidth = 3 )
	img += plt.plot( x/nm, vx7_real, colors[4], linestyle='solid', linewidth = 3 )
	img += plt.plot( x/nm, vx7_imag, colors[5], linestyle='solid', linewidth = 3 )
	
	#アニメーションに追加
	ims.append( img )

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

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, 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()