#################################################################
## ガウス波束の運動
#################################################################
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.62607015 * 1.0E-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.1093837015 * 1.0E-31
#電子ボルト
eV = 1.602176634 * 1.0E-19
#################################################################
## 物理系の設定
#################################################################
#波束の中心エネルギー
E_bar = 1.0 * eV
#重ね合わせの数
NK = 200
#ガウス分布の幅
sigma = 1.0 / ( 1.25 * 10**-9 )
#波数の間隔
dk = 10.0 * sigma / NK
#ガウス分布の中心
k_bar = np.sqrt( 2.0 * me * E_bar / hbar**2)
#空間分割サイズ
dx = 1.0 * 10**-9
#空間分割数
NX = 500
#計算区間
x_min = -10.0 * dx
x_max = 10.0 * dx
#計算時間の幅
ts = -100
te = 100
#時間間隔
dt = 4.0 * 10**-16
#ガウス波束の値を計算する関数
def Psi( x, t ):
#波動関数値の初期化
psi = x * (0.0 + 0.0j)
#各波数ごとの寄与を足し合わせる
for kn in range(NK):
#各波数を取得
k = k_bar + dk * (kn - NK/2)
#波数から各振動数を取得
omega = hbar / (2.0 * me) * k**2
#平面波を足し合わせる
psi += np.exp(I * k * x ) * np.exp(- I * omega * t) * np.exp( -( (k - k_bar) / (2.0 * sigma) )**2 )
return psi
#基準となる振幅を取得
psi_abs = abs(Psi( 0, 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"ガウス波束の運動", 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.xlim([x_min/dx, x_max/dx])
plt.ylim([-1.1, 1.1])
#アニメーションの生成
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 = k_bar - dk * NK / 2.0
k_max = k_bar + dk * NK / 2.0
#波数座標点配列の生成
k = np.linspace(k_min, k_max, NK)
#正規分布
c_n = np.exp( -((k - k_bar) / (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- k_bar)/sigma, (k_max- k_bar)/sigma])
plt.ylim([0, 1.05])
#波数分布グラフの描画
plt.plot((k-k_bar)/sigma, c_n, linestyle='solid', linewidth = 5)
#グラフの表示
plt.show()