#################################################################
## 透過率と反射率の入射波のエネルギー依存性のグラフ描画（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
######################################
# 物理系の設定
######################################
#空間刻み間隔
nm = 1E-9
#壁の厚さ
d = 2.0 * nm
#壁の高さ
V1 = 0 * eV
V2 = 1.0 * eV
V3 = 0 * eV
######################################
# 反射係数と透過係数を計算する関数定義
######################################
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, 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 / 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
M32[0][0] = (1.0 + 0.0j + k2 / k3) / 2.0
M32[0][1] = (1.0 + 0.0j - k2 / k3) / 2.0
M32[1][0] = (1.0 + 0.0j - k2 / k3) / 2.0
M32[1][1] = (1.0 + 0.0j + k2 / k3) / 2.0
T2[0][0] = np.exp( I * k2 * d )
T2[0][1] = 0
T2[1][0] = 0
T2[1][1] = np.exp( - I * k2 * d )
#転送行列の計算
M31 = M32 @ T2 @ M21;
#反射係数と透過係数の計算
rc = ReflectionCoefficient(M31);
tc = TransmissionCoefficient(M31);
return rc, tc
#################################################################
## 透過率と反射率の入射波のエネルギー依存性のグラフ描画
#################################################################
#エネルギーの計算範囲
E_min = 0.001
E_max = 3.0
#計算点数
NK = 500
#エネルギーの行列
Es = np.linspace(E_min, E_max, NK)
#反射率と透過率の配列
Rs = []
Ts = []
for E in Es:
E = E * 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 ))
k3 = np.sqrt(np.complex128( 2.0 * me * (E - V3) / hbar**2 ))
rc, tc = calculateCoefficient( k1, k2, k3, d)
Rs.append( abs(rc)**2 )
Ts.append( abs(tc)**2 )
#グラフの描画（固有関数）
plt.title( u"透過率と反射率の入射波のエネルギー依存性のグラフ（壁の厚さ：" + r"$2.0[{\rm nm}]$" + ", 壁の高さ：" + r"$1.0[{\rm eV}]$" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(u"入射波のエネルギー" + r"$[{\rm eV}]$" , 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, 3])
plt.ylim([0, 1])
#波数分布グラフの描画
plt.plot(Es, Rs, linestyle='solid', linewidth = 5)
plt.plot(Es, Ts, linestyle='solid', linewidth = 5)
#グラフの表示
plt.show()
【第11話】Q．壁の高さがマイナスの場合、波束の運動はどうなる？【Pythonコピペで量子力学完全攻略マニュアル】
【第10話】Q．中心エネルギーが20[eV]のガウス波束を、中心エネルギーと同じ、20[eV]の壁に衝突させたらどうなるだろう？【Pythonコピペで量子力学完全攻略マニュアル】
【第９話】有限の高さの障壁に向けて照射した電子パルスのアニメーション【Pythonコピペで量子力学完全攻略マニュアル】
#############################################################
## 有限の高さの障壁に向けてガウス波束を照射したときの時間発展（E = 10.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 = 20.0 * eV
#重ね合わせの数
NK = 1000
#ガウス分布の幅
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 = 10.0 * dx
#壁の高さ
V = 10.0 * eV
#ガウス波束の初期位置
x0 = -7.5 * dx
#計算時間の幅
ts = 0
te = 300
#時間間隔
dt = 0.3 * 10**-16
#ガウス波束の値を計算する関数
def Psi( x, t, region ):
#波動関数値の初期化
psi = x * (0.0 + 0.0j)
#各波数ごとの寄与を足し合わせる
for kn in range(NK):
#各波数を取得
k1 = k0 + dk * (kn - NK/2)
#波数から各振動数を取得
omega = hbar / (2.0 * me) * k1**2
#平面波のエネルギー
E = (hbar * k1)**2 / (2.0 * me)
#領域IIの波数
k2 = np.sqrt(np.complex128( 2.0 * me * (E - V) / hbar**2 ))
#反射係数
rc = ( k1 - k2 ) / ( k1 + k2 )
#透過係数
tc = 2.0 * k1 / ( k1 + k2 )
#入射波と反射波と透過波
phi_I = np.exp( I * k1 * x - I * omega * t )
phi_R = rc * np.exp( - I * k1 * x - I * omega * t )
phi_T = tc * np.exp( I * k2 * x - I * omega * t )
#ガウス分布
Ck = np.exp( - I * k1 * x0 ) * np.exp( -( (k1 - k0) / (2.0 * sigma) )**2 )
#平面波を足し合わせる
if(region == 1):
psi += (phi_I + phi_R) * Ck
elif(region == 2):
psi += phi_T * Ck
return psi
#基準となる振幅を取得
psi_abs = abs(Psi( x0, 0, 1 ))
#座標点配列の生成
x1 = np.linspace(x_min, 0, NX)
x2 = np.linspace(0, x_max, NX)
#アニメーション作成用
ims=[]
### 波束の空間分布の計算
for tn in range(ts, te + 1):
#実時間の取得
t = dt * tn
#時刻tの波動関数を取得
psi1 = Psi( x1, t, 1 )
psi2 = Psi( x2, t, 2 )
#波動関数の規格化
psi1 /= psi_abs
psi2 /= psi_abs
#各コマを描画
img = plt.plot(x1/dx, psi1.real, colors[0], linestyle='solid', linewidth = 3.0)
img += plt.plot(x1/dx, psi1.imag, colors[1], linestyle='solid', linewidth = 3.0)
img += plt.plot(x1/dx, abs(psi1), colors[2], linestyle='solid', linewidth = 3.0)
img += plt.plot(x2/dx, psi2.real, colors[3], linestyle='solid', linewidth = 3.0)
img += plt.plot(x2/dx, psi2.imag, colors[4], linestyle='solid', linewidth = 3.0)
img += plt.plot(x2/dx, abs(psi2), colors[5], linestyle='solid', linewidth = 3.0)
#time = plt.text( 8.5, 1.57, "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_0 = 20.0 [{\rm eV}] , V = -20.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], +0.025, -1.025, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0], -10, 0, "black", linestyles='solid', linewidth = 10 )
plt.hlines([-1], 0, 10, "black", linestyles='solid', linewidth = 10 )
#描画範囲を設定
plt.xlim([x_min/dx, x_max/dx])
plt.ylim([-1.1, 1.5])
#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=50)
#アニメーションの保存
ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
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()
【第６話】無限に高い障壁へ照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】
#################################################################
## 無限に高い障壁へ照射した平面波の時間発展（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
######################################
# 物理系の設定
######################################
#電子のエネルギー
E1 = 1.0 * eV
#波数
k1 = np.sqrt(2.0 * me * E1 / hbar**2 )
#角振動数
omega1 = E1/hbar
#１周期
T1 = 2.0 * np.pi / omega1
#時間間隔
#dt = 0.2 * 10**-16
dt = T1 / 100
#時間刻み数
NT = 1000
#空間刻み間隔
dx = 1E-9
#描画範囲
x_min = -5.0 * dx
x_max = 0.0 * dx
#描画区間数
NX = 500
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#アニメーション作成用
ims=[]
#各時刻における計算
for tn in range(NT):
t = dt * tn
#波動関数の計算
phi_I = np.exp( I * k1 * x - I * omega1 * t )
phi_R = np.exp( - I * k1 * x - I * omega1 * t )
phi = phi_I - phi_R
#各コマを描画
img = plt.plot(x/dx, phi_I.real, colors[0], linestyle='solid', linewidth = 5.0)
img += plt.plot(x/dx, phi_R.real, colors[1], linestyle='solid', linewidth = 5.0)
img += plt.plot(x/dx, phi.real, colors[2], 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], -5, 5, "black", linestyles='solid', linewidth = 10 )
#描画範囲を設定
plt.xlim([x_min/dx, x_max/dx + 0.5])
plt.ylim([-3.1, 3.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()
【第５話】電子パルスの運動アニメーション【Pythonコピペで量子力学完全攻略マニュアル】
#################################################################
## ガウス波束の運動
#################################################################
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()
【第４話】電子パルスの作り方【Pythonコピペで量子力学完全攻略マニュアル】
#################################################################
## ガウス波束の空間分布
#################################################################
import numpy as np
import matplotlib.pyplot as plt
#図全体
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
#################################################################
## 物理系の設定
#################################################################
#波束の中心エネルギー
E0 = 10.0 * eV
#重ね合わせの数
NK = 500
#ガウス分布の幅
sigma = 1.0 / ( 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 = -5.0 * dx
x_max = 5.0 * dx
#ガウス波束の値を計算する関数
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
#平面波を足し合わせる
psi += np.exp(I * k * x ) * np.exp(- I * omega * t) * np.exp( -( (k - k0) / (2.0 * sigma) )**2 )
return psi
#基準となる振幅を取得
psi_abs = abs(Psi( 0, 0 ))
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#時刻t=0の波動関数を取得
psi = Psi( x, 0 )
#波動関数の規格化
psi /= psi_abs
#グラフの描画（波動関数）
plt.title( u"ガウス波束", fontsize=20, fontname="Yu Gothic", fontweight=1000 )
plt.xlabel( r"$ x [{\rm nm} ] $", fontsize=30 )
plt.ylabel( r"$ \psi(x,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([x_min/dx, x_max/dx])
plt.ylim([-1.1, 1.1])
#グラフの描画
plt.plot(x/dx, psi.real, linestyle='solid', linewidth = 5)
plt.plot(x/dx, psi.imag, linestyle='solid', linewidth = 5)
plt.plot(x/dx, abs(psi), linestyle='solid', linewidth = 5)
#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.15, top = 0.95)
#################################################################
## ガウス分布の描画
#################################################################
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_k = np.exp( -((k - k0) / (2.0 * sigma))**2 )
#グラフの描画（固有関数）
plt.title( u"ガウス分布（正規分布）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$ (k - k_0)/\sigma $", fontsize=30)
plt.ylabel(r"$ C(k)/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.subplots_adjust(left = 0.1, right = 0.98, bottom=0.15, top = 0.95)
#描画範囲を設定
plt.xlim([(k_min- k0)/sigma, (k_max- k0)/sigma])
plt.ylim([0, 1.05])
#波数分布グラフの描画
plt.plot((k-k0)/sigma, C_k, linestyle='solid', linewidth = 5)
#グラフの表示
plt.show()
【第３話】自由粒子の運動アニメーション【Pythonコピペで量子力学完全攻略マニュアル】
#################################################################
## 自由粒子の波動関数の時間発展（E = 1.0[eV], 0.25[eV], 4.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
######################################
# 物理系の設定
######################################
#電子のエネルギー
E1 = 1.0 * eV
E2 = 0.25 * eV
E3 = 4.0 * eV
#波数
k1 = np.sqrt(2.0 * me * E1 / hbar**2 )
k2 = np.sqrt(2.0 * me * E2 / hbar**2 )
k3 = np.sqrt(2.0 * me * E3 / hbar**2 )
#角振動数
omega1 = E1/hbar
omega2 = E2/hbar
omega3 = E3/hbar
#時間間隔
dt = 1 * 10**-16
#時間刻み数
NT = 1000
#空間刻み間隔
dx = 1E-9
#描画範囲
x_min = -2.0 * dx
x_max = 2.0 * dx
#描画区間数
NX = 500
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#アニメーション作成用
ims=[]
#各時刻における計算
for tn in range(NT):
t = dt * tn
#波動関数の計算
phi1 = np.exp( I * k1 * x - I * omega1 * t )
phi2 = np.exp( I * k2 * x - I * omega2 * t )
phi3 = np.exp( I * k3 * x - I * omega3 * t )
#各コマを描画
img = plt.plot(x/dx, phi1.real, colors[0], linestyle='solid', linewidth = 5.0)
img += plt.plot(x/dx, phi2.real, colors[1], linestyle='solid', linewidth = 5.0)
img += plt.plot(x/dx, phi3.real, colors[2], linestyle='solid', linewidth = 5.0)
ims.append( img )
#グラフの描画（波動関数）
plt.title( u"自由粒子の波動関数（" + r"$ E = 0.25, 1.0, 4.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.xlim([x_min/dx, x_max/dx])
plt.ylim([-1.1, 1.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()
【第２話】自由粒子のスナップショット【Pythonコピペで量子力学完全攻略マニュアル】
#################################################################
## 自由粒子の固有関数（E = 0.25[eV], 1.0[eV], 4.0[eV]）
#################################################################
import numpy as np
import matplotlib.pyplot as plt
#図全体
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
#################################################################
## 物理系の設定
#################################################################
#電子のエネルギー
E1 = 1.0 * eV
E2 = 0.25 * eV
E3 = 4.0 * eV
#波数
k1 = np.sqrt(2.0 * me * E1 / hbar**2)
k2 = np.sqrt(2.0 * me * E2 / hbar**2)
k3 = np.sqrt(2.0 * me * E3 / hbar**2)
#空間刻み間隔
dx = 1E-9
#描画範囲
x_min = -2.0 * dx
x_max = 2.0 * dx
#描画区間数
NX = 500
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#平面波の計算
phi1 = np.exp( I * k1 * x )
phi2 = np.exp( I * k2 * x )
phi3 = np.exp( I * k3 * x )
#グラフの描画（固有関数）
plt.title( u"自由粒子の固有関数（" + r"$ E = 0.25, 1.0, 4.0[{\rm eV}] $" + u"）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$ x\, [{\rm nm}]$ ", fontsize=30)
plt.ylabel(r"$ {\rm Re}\{\varphi(x)\} $", 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([x_min/dx, x_max/dx])
plt.ylim([-1.1, 1.1])
#グラフの描画
plt.plot(x/dx, phi1.real, linestyle='solid', linewidth = 5)
plt.plot(x/dx, phi2.real, linestyle='solid', linewidth = 5)
plt.plot(x/dx, phi3.real, linestyle='solid', linewidth = 5)
#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.15, top = 0.95)
#グラフの表示
plt.show()
【第１話】プログラムの動作確認【Pythonコピペで量子力学完全攻略マニュアル】
#################################################################
## Pythonでグラフアニメーション
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
#図全体
fig = plt.figure(figsize=(10, 6))
#描画点数
NX = 100
#描画範囲
x_min = -1.0
x_max = 1.0
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#アニメーション分割数
AN = 30
#アニメーション作成用
ims = []
#アニメーションの各グラフを生成
for a in range(AN):
phi = 2.0 * np.pi * a / AN
#sin関数値を計算
y = np.sin( np.pi * x + phi )
#各コマのグラフの描画
img = plt.plot(x, y, color="blue", linewidth=3.0, linestyle="solid" )
ims.append( img )
#グラフタイトルの設定
plt.title("sin function")
#x軸ラベル・y軸ラベルの設定
plt.xlabel("x-axis")
plt.ylabel("y-axis")
#描画範囲を設定
plt.xlim([-1.0,1.0])
plt.ylim([-1.0,1.0])
#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=50)
#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)
#ani.save("output.mp4", writer=matplotlib.animation.AVConvWriter, dpi=300)
#グラフの表示
plt.show()