2021年9月 – シミュ君日記

2021年9月23日

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

2021年9月19日

【第８話】サイクロイド曲線に経路が拘束された粒子の運動【ルンゲクッタで行こう！⑧】

############################################################################
#　[008]サイクロイド曲線に経路が拘束された粒子の運動
############################################################################
import os
import numpy as np
import random 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
random.seed(0)

#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
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']
#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
#時間区間
t_min = 0
t_max = 20.0
#時間区間数
NT =   20000
skip =    20

#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#重力加速度
g = np.array([0.0, 0.0, -10.0])

#################################################################
## 経路の設定
#################################################################
class Path:
	#コンストラクタ
	def __init__(self):
		#始点と終点のx座標
		self.x_start = -1
		self.x_end = 1
		#サイクロイド曲線パラメータr
		self.r = ( self.x_end - self.x_start ) / ( 2 * np.pi )
		#始点と終点のz座標
		self.path_z0 = 2 * self.r ; #最下点が0となるように
		#サイクロイド曲線始点の位置ベクトル
		self.start = np.array([self.x_start, 0.0, self.path_z0])

	#経路の位置ベクトルを指定する媒介変数関数
	def position( self, theta ):
		x = self.r * ( theta - np.sin(theta) ) + self.start[0]
		y = 0
		z = - self.r * ( 1 - np.cos(theta) ) + self.start[2]
		return np.array([x, y, z])

	#接線ベクトルを指定する媒介変数関数
	def tangent( self, theta ):
		A = 1.0 / np.sqrt( 2.0 )
		x = A * np.sqrt( 1.0 - np.cos(theta))
		y = 0
		z = - A * np.sqrt( 1.0 + np.cos(theta))
		if( np.sin(theta) < 0 ): z = - z
		return np.array([x, y, z])

	#曲率ベクトルを指定する媒介変数関数
	def curvature( self, theta ):
		A = - 1.0 / ( 4.0 * self.r )
		x = - A * np.sqrt( ( 1.0 + np.cos(theta) ) / ( 1.0 - np.cos(theta) ) )
		y = 0
		z = - A
		if( np.abs(x) == np.inf ):
			x = 0
			z = 0
		if( np.sin(theta) < 0 ): x = - x

		return np.array([x, y, z])
	#媒介変数の取得
	def getTheta( self, r, v ):
		#球体の相対位置ベクトル
		bar_r = r -self.start
		if( np.abs( v[0] ) > 0.1 ):
			_r = np.abs(bar_r[2]) / 2.0 * (  1.0 + (v[2]/v[0])**2 )
		else :
			_r = self.r

		if( bar_r[0] < _r * np.pi ) :
			theta = np.arccos( 1.0 + bar_r[2]/_r )
		else:
			theta = 2.0 * np.pi - np.arccos( 1.0 + bar_r[2]/_r )


		return theta

#経路の生成
path = Path()

#補正パラメータ
compensationK = 100.0;     #補正ばね弾性係数
compensationGamma = 1.0; #補正粘性抵抗係数
compensationFactor = 1;  #補正因子

#粒子に加わる力
def F (self, t, r_t, v_t ):
	#外場（重力）の計算
	fe = m * g

	#経路自体の運動
	v0 = np.array([0.0, 0.0, 0.0])
	a0 = np.array([0.0, 0.0, 0.0])

	#相対速度
	bar_v = v_t - v0

	#媒介変数の取得
	theta = path.getTheta ( r_t, bar_v )

	#媒介変数に対する位置ベクトル、接線ベクトル、曲率ベクトルの計算
	position = path.position( theta )
	tangent = path.tangent( theta )
	curvature = path.curvature( theta )

	#微係数dl/dtとd^2l/dt^2を計算
	dl_dt = np.dot( bar_v, tangent )
	d2l_dt2 = np.dot( fe, tangent ) / m - np.dot( a0, tangent) 

	#加速度を計算
	a = a0 + d2l_dt2 * tangent + dl_dt * dl_dt * curvature
	a_length = np.sqrt(np.sum( a**2 ))

	#補正倍率
	ratio = a_length * compensationFactor
	#ズレ位置ベクトル
	bar_r = r_t - position
	#補正ばね弾性力
	fk = -compensationK * ratio * bar_r
	#法線ベクトル
	curvature_length = np.sqrt(np.sum( curvature**2 ))
	n1 = curvature / curvature_length
	n2 = np.cross(n1, tangent )
	#補正粘性抵抗力
	fgamma1 = -compensationGamma * np.dot( bar_v, n1 ) * ratio * n1 
	fgamma2 = -compensationGamma * np.dot( bar_v, n2 ) * ratio * n2
	#経路補正力を加える
	a += fk + fgamma1 + fgamma2

	return m * a

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		#現時刻の位置と速度
		self.r = r0.copy()
		self.v = v0.copy()
		#１ステップ前の時刻の位置と速度
		self._r = r0.copy()
		self._v = v0.copy()
		#差分
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		#時間間隔
		self.dt = dt
		#内部
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F(self, t, r_t, v_t ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始
#################################################################
#初期条件
particles =[]
for i in range(20):
	#球体の初期位置
	theta0 = np.pi
	#初速度
	x0 = path.r * ( theta0 - np.sin(theta0) ) + path.start[0]
	y0 = 0
	z0 = - path.r * ( 1.0 - np.cos(theta0)) + path.start[2]
	#初期条件
	r0 = np.array([ x0, y0, z0 ])
	v0 = np.array([ 0.170 * (i+1), 0.0, 0.0 ])

	particles.append(  RK4( m, dt, r0, v0 ) )

####################################
## 描画用関数
####################################
#アニメーション作成用
ims=[]

#円座標点配列の生成
thetas = np.linspace(0, 2.0*np.pi, 200)
p = path.position( thetas )
cirsx = p[0]
cirsy = p[2]

def plot( t ):
	#各コマを描画
	img = plt.plot(cirsx, cirsy, color = "black", linestyle='dotted', linewidth = 1 )
	i = 0
	for particle in particles:
		img += plt.plot([particle.r[0]], [particle.r[2]], color = colors[i], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
		i += 1
		if (i>9): i = 0

	time = plt.text(0.8, 0.717, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )
	E = 1.0/2.0 * m * np.sum( particles[0].v**2 ) - m * np.dot( g, particles[0].r )
	print( "t = {:.4f}".format(t)," E =",  E )

####################################
## 各時刻における運動方程式の計算
####################################
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t)

	#ルンゲ・クッタ法による時間発展
	for particle in particles:
		particle.timeEvolution( t )
	#位置と速度を更新
	for particle in particles:
		particle.r += particle.dr
		particle.v += particle.dv

#################################################################
## 分子運動アニメーション
#################################################################
plt.title( u"サイクロイド曲線経路に拘束された球体の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", 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([-1.1, 1.1])
plt.ylim([-0.1, 0.7])
plt.subplots_adjust(left = 0.12, right = 0.97, bottom = 0.10, top = 0.90)

#アニメーションの生成
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()

2021年9月17日2021年9月19日

【第７話】経路が拘束された粒子の運動【ルンゲクッタで行こう！⑦】

############################################################################
#　[007]経路に拘束された剛体球の運動
############################################################################
import os
import numpy as np
import random 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
random.seed(0)

#図全体
fig = plt.figure(figsize=(8, 8))
#全体設定
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']
#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
#時間区間
t_min = 0
t_max = 10.0
#時間区間数
NT =   1000
skip =   10
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#重力加速度
g = np.array([0.0, 0.0, -10.0])

#################################################################
## 経路の設定
#################################################################
class Path:
	#コンストラクタ
	def __init__(self):
		#円の半径
		self.R = 1
		#円の中心位置ベクトル
		self.center = np.array([0.0, 0.0, self.R])

	#経路の位置ベクトルを指定する媒介変数関数
	def position( self, theta ):
		x = self.R * np.cos(theta)
		y = 0
		z = self.R * np.sin(theta) + self.center[2]
		return np.array([x, y, z])

	#接線ベクトルを指定する媒介変数関数
	def tangent( self, theta ):
		x = -np.sin(theta)
		y = 0
		z = np.cos(theta)
		return np.array([x, y, z])
	#曲率ベクトルを指定する媒介変数関数
	def curvature( self, theta ):
		x = -np.cos(theta) / self.R
		y = 0
		z = -np.sin(theta) / self.R
		return np.array([x, y, z])
	#媒介変数の取得
	def getTheta( self, r ):
		#相対位置ベクトル
		bar_r = r - self.center
		bar_r_length = np.sqrt(np.sum( bar_r**2 ))
		sinTheta = bar_r[2] / bar_r_length
		if( sinTheta > 0 ):
			theta = np.arccos( bar_r[0] / bar_r_length )
		else:
			theta = 2.0 * np.pi - np.arccos( bar_r[0] / bar_r_length )

		return theta

#経路の生成
path = Path()

#補正パラメータ
compensationK = 1.0;     #補正ばね弾性係数
compensationGamma = 1.0; #補正粘性抵抗係数
compensationFactor = 0;  #補正因子

#粒子に加わる力
def F (self, t, r_t, v_t ):
	#外場（重力）の計算
	fe = m * g

	#媒介変数の取得
	theta = path.getTheta ( r_t )

	#媒介変数に対する位置ベクトル、接線ベクトル、曲率ベクトルの計算
	position = path.position( theta )
	tangent = path.tangent( theta )
	curvature = path.curvature( theta )

	#経路自体の運動
	v0 = np.array([0.0, 0.0, 0.0])
	a0 = np.array([0.0, 0.0, 0.0])

	#相対速度
	bar_v = v_t - v0

	#微係数dl/dtとd^2l/dt^2を計算
	dl_dt = np.dot( bar_v, tangent )
	d2l_dt2 = np.dot( fe, tangent ) / m - np.dot( a0, tangent) 

	#加速度を計算
	a = a0 + d2l_dt2 * tangent + dl_dt * dl_dt * curvature
	a_length = np.sqrt(np.sum( a**2 ))

	#補正倍率
	ratio = a_length * compensationFactor
	#ズレ位置ベクトル
	bar_r = r_t - position
	#補正ばね弾性力
	fk = -compensationK * ratio * bar_r
	#法線ベクトル
	curvature_length = np.sqrt(np.sum( curvature**2 ))
	n1 = curvature / curvature_length
	n2 = np.cross(n1, tangent )
	#補正粘性抵抗力
	fgamma1 = -compensationGamma * np.dot( bar_v, n1 ) * ratio * n1 
	fgamma2 = -compensationGamma * np.dot( bar_v, n2 ) * ratio * n2
	#経路補正力を加える
	a += fk + fgamma1 + fgamma2

	return m * a

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		#現時刻の位置と速度
		self.r = r0.copy()
		self.v = v0.copy()
		#１ステップ前の時刻の位置と速度
		self._r = r0.copy()
		self._v = v0.copy()
		#差分
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		#時間間隔
		self.dt = dt
		#内部
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F(self, t, r_t, v_t ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始
#################################################################
#初期条件
particles =[]
for i in range(20):
	particles.append(  RK4( m, dt, np.array([ 0.0, 0.0, 2.0 ]), np.array([ 0.001 * (i+1), 0.0, 0.0 ] ) ) )

####################################
## 描画用関数
####################################
#アニメーション作成用
ims=[]

#円座標点配列の生成
thetas = np.linspace(0, 2*np.pi, 200)
cirsx = np.cos( thetas )
cirsy = np.sin( thetas ) + 1.0

def plot( t ):
	#各コマを描画
	img = plt.plot(cirsx, cirsy, color = "black", linestyle='dotted', linewidth = 1 )
	i = 0
	for particle in particles:
		img += plt.plot([particle.r[0]], [particle.r[2]], color = colors[i], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
		i += 1
		if (i>9): i = 0

	time = plt.text(0.8, 2.18, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )
	E = 1.0/2.0 * m * np.sum( particles[0].v**2 ) - m * np.dot( g, particles[0].r )
	print( "t = {:.4f}".format(t)," E =",  E )

####################################
## 各時刻における運動方程式の計算
####################################
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t)

	#ルンゲ・クッタ法による時間発展
	for particle in particles:
		particle.timeEvolution( t )
	#位置と速度を更新
	for particle in particles:
		particle.r += particle.dr
		particle.v += particle.dv

#################################################################
## 分子運動アニメーション
#################################################################
plt.title( u"円経路に拘束された球体の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", 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([-1.1, 1.1])
plt.ylim([-0.1, 2.1])
plt.subplots_adjust(left = 0.12, right = 0.97, bottom = 0.10, top = 0.90)

#アニメーションの生成
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()

2021年9月14日

【第31話】調和ポテンシャル中の波動関数の初期配置の与え方【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 調和ポテンシャル中の電子の固有関数　初期値の与え方
#################################################################
import math
import numpy as np
import scipy.integrate as integrate
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

#################################################################
## 物理系の設定
#################################################################
#量子井戸の幅（m）
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

#状態数
NS = 100

#################################################################
## 初期分布：ガウス分布
#################################################################
sigma = 1.0 * 10**10
x0 = 1.0 * nm
def verphi0(x):
	return ( sigma**2 / np.pi )**(1.0/4.0) * np.exp( - 1.0/2.0 * sigma**2 * (x - x0)**2 )

#被積分関数
def integral_orthonomality(x, n):
	return verphi(n, x) * verphi0(x)

L = 7 * nm

#描画範囲
x_min = -L/2
x_max = L/2
cn = []
for n in range( NS + 1 ):
	#ガウス・ルジャンドル積分
	result = integrate.quad(
		integral_orthonomality, #被積分関数
		x_min, x_max,           #積分区間の下端と上端
		args = ( n )           #被積分関数へ渡す引数
	)
	cn.append( result[0] )	
	#ターミナルへ出力
	print( "(" + str(n) + ")  " + str( result[0]) )

#################################################################
## 波動関数
#################################################################
def psi(x, t):
	_psi = 0.0 + 0.0j
	for n in range(NS):
		_psi += cn[n] * verphi(n, x) * np.exp( - I * ( n + 1.0/2.0) * omega * t )
	return _psi	

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

#基底状態の周期
T0 = 2.0 * np.pi * hbar / Energy(0)
print( "基底状態の周期：" + str(T0) + " [s]" )

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

plt.title( u"調和ポテンシャル中の波動関数", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=20)

L = 7 * 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([-3.5, 3.5])
plt.ylim([-1.5, 6.5])

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

ims = []
#各時刻における波動関数の計算
for tn in range(ts, te):
	#実時間の取得
	t = dt * tn
	#波動関数の計算
	ys = psi(x, t).real / np.sqrt(A) * 2+ yE

	#波動関数の表示
	img  = plt.plot( x/nm, yE, linestyle='dotted', color="black", linewidth = 1 )
	img += plt.plot( x/nm, ys, colors[0], 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.05, 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()

2021年9月9日

【はやくち解説】エルミート多項式の有難みって？【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 規格化エルミート多項式
#################################################################
import math
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
#図全体
fig = plt.figure(figsize=(9, 9))
#全体設定
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']

#インデックス最大数
NS = 6

#################################################################
## エルミート多項式
#################################################################
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

#################################################################
## 直交性の確認
#################################################################
#被積分関数
def integral_orthonomality(x, n, m):
	return NormalizedHermitePolynomial( m, x ) * NormalizedHermitePolynomial(n, x)

for n in range( NS + 1 ):
	for m in range( NS + 1 ):
		#ガウス・ルジャンドル積分
		result = integrate.quad(
			integral_orthonomality, #被積分関数
			-10, 10,           #積分区間の下端と上端
			args=( n, m )           #被積分関数へ渡す引数
		)
		#ターミナルへ出力
		print( "(" + str(n) + ", " + str(m) + ")  " + str( result[0]) )


#########################################################################
# グラフの描画（固有関数）
#########################################################################
plt.title( u"規格化エルミート多項式", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x$", fontsize=30)
#描画区間数
NX = 500
#描画範囲
x_min = -5
x_max = 5
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)

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

#描画範囲を設定
plt.xlim([x_min , x_max])
plt.ylim([-0.5, 10])
plt.yticks(color="None")
plt.tick_params(length=0)

#グラフの描画
ys = [ 0 ] * (NS + 1)
for n in range( NS + 1 ):
	ys[n] = NormalizedHermitePolynomial( n, x ) + n * 1.5

for n in range( NS + 1 ):
	plt.plot(x, ys[n], linestyle='solid', linewidth = 5 )
	plt.text( 5.1, n * 1.5 - 0.12, r"$H_" + str(n) + "$", fontsize = 25)
	plt.hlines([n * 1.5],  -5, 5, "#000000", linestyles='dotted', linewidth = 1 )

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

#グラフの表示
plt.show()

2021年9月7日2021年9月14日

【第29話】調和ポテンシャル中の電子状態【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 調和ポテンシャル中の電子の固有関数
#################################################################
import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
#図全体
fig = plt.figure(figsize=(8, 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

#################################################################
## 物理系の設定
#################################################################
#量子井戸の幅（m）
omega = 1.0 * 1.0E+15
A = np.sqrt( me * omega / hbar)
#エルミート多項式
def HermitePolynomial( n, x ):
	x0 = 1
	x1 = 2.0 * x
	if(n==0): return x0
	if(n==1): return x1

	for m in range(2, n+1):
		x2 = 2.0 * x * x1 - 2.0 * ( m - 1 ) * x0
		x0 = x1
		x1 = x2

	return x2

#固有関数
def verphi(n, x):
	barX = A * x
	return ( A**2 / np.pi) **(1.0/4.0) * 1.0 / np.sqrt( 2**n * math.factorial(n) ) * np.exp( - 1.0/2.0 * barX **2 ) * HermitePolynomial( n, barX )
#固有エネルギー（eV）
def Energy(n):
	return (n + 1.0/2.0) * hbar * omega

#波動関数
def psi(n, x, t):
	return verphi(n, x) * np.exp( - I * ( n + 1.0 / 2.0 ) * omega * t )

#状態数
NS = 4
#固有エネルギーの表示（eV）
for n in range( NS + 1 ):
	if( n == 0 ): text = "基底状態"
	else: text = "第" + str(n) + "励起状態"
	print( text + "：" + str(Energy( n )/eV) + " [eV]" )

#計算時間の幅
ts = 0
te = 100

#基底状態の周期
T0 = 2.0 * np.pi * hbar / Energy(0)
print( "基底状態の周期：" + str(T0) + " [s]" )


#########################################################################
# 波動関数 アニメーション
#########################################################################
#時間間隔
dt = T0 / (te - ts + 1)

plt.title( u"調和ポテンシャル中の電子の固有関数", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\, [{\rm nm}]$", fontsize=20)

L = 6 * 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([-3, 3])
plt.ylim([-0.5, 7.5])

NS = 6

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

#各時刻における波動関数の計算
for tn in range(ts, te):
	#実時間の取得
	t = dt * tn

	img  = plt.plot( x/nm, yE, linestyle='dotted', color="black", linewidth = 1 )

	#各コマを描画
	for n in range( NS + 1 ):
		#波動関数の計算
		y = psi(n, x, t).real / np.sqrt(A) * 0.8 + n + 0.5
		#波動関数の表示
		img += plt.plot( x/nm, y, colors[n], linestyle='solid', linewidth = 3 )

	#アニメーションに追加
	ims.append( img )


#E_0～E_6の表示
for n in range( 0, NS + 1 ):
	plt.text( 3.1, 0.5 + n - 0.12, r"$E_" + str(n) + "$", fontsize = 25)
	plt.hlines([0.5 + n],  -3, 3, "#CCCCCC", linestyles='dotted', linewidth = 1 )

#余白の調整
#plt.subplots_adjust(left=0.15, right=0.90, bottom=0.1, top=0.99)
plt.subplots_adjust(left = 0.05, right = 0.92, 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()

2021年9月6日2021年9月6日

【第６話】剛体球同士の衝突力計算アルゴリズム【ルンゲクッタで行こう！⑥】

############################################################################
#　箱の中のN個の自由粒子（古典力学）
############################################################################
import os
import numpy as np
import random 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
random.seed(0)

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
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']
#################################################################
## 物理系の設定
#################################################################
#質量
m_min = 1.0
m_max = 10.0
#粒子の半径
radius_min = 0.002
radius_max = 0.02
#粒子数
NM = 100
#箱の長さ
L = 1.0
#粒子に加わる力
fex = np.array([0.0, 0.0, 0.0])
#初期位置と速度
v_max = 2.0

#粒子に加わる力
def F (self, t, r_t, v_t ):
	return self.fex + self.fc - self.m * self.beta * v_t

#時間区間
t_min = 0
t_max = 5.0
#時間区間数
NT =    500
skip =    1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#
lw = 0.01
######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		#現時刻の位置と速度
		self.r = r0.copy()
		self.v = v0.copy()
		#１ステップ前の時刻の位置と速度
		self._r = r0.copy()
		self._v = v0.copy()
		#差分
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		#時間間隔
		self.dt = dt
		#外力
		self.fex = np.array([0.0, 0.0, 0.0])
		#衝突力
		self.fc = np.array([0.0, 0.0, 0.0])
		#粘性係数
		self.beta = 0
		#半径
		self.radius = 0
		#内部
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F(self, t, r_t, v_t ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 時間発展の計算
#################################################################

#壁との衝突判定（戻値：法線ベクトル）
def CollisionWallDetection( particle ):
	#衝突リスト
	collisionlist = []
	
	if(particle.r[0] + particle.radius > L):
		collisionlist.append({ "type": 1, "normal" : np.array([-1.0, 0.0, 0.0]) })
	if(particle.r[0] - particle.radius < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([1.0, 0.0, 0.0]) })
	if(particle.r[1] + particle.radius > L):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, -1.0, 0.0]) })
	if(particle.r[1] - particle.radius < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, 1.0, 0.0]) })
	if(particle.r[2] + particle.radius > L):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, 0.0, -1.0]) })
	if(particle.r[2] - particle.radius < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, 0.0, 1.0]) })

	return collisionlist

#粒子の運動エネルギー
def T():
	_T = 0
	for nm in range(NM):
		_T += 1.0/2.0 *  particles[nm].m * np.sum( particles[nm].v**2 )
	return  _T
#粒子の位置エネルギー
def U():
	_R = 0
	for nm in range(NM):
		_R += - np.dot(particles[nm].fex, particles[nm].r)
	return _R

#系の時間発展の計算
def timeEvolution( t, timeScale ):

	#複数衝突
	cylinderCollisionList = []

	for nm in range(NM):
		#時間スケールの指定
		particles[nm].dt = dt / timeScale

		#ルンゲ・クッタ法による時間発展
		particles[nm].fc = np.array([0.0, 0.0, 0.0])
		particles[nm].timeEvolution( t )
		#位置と速度を更新
		particles[nm].r += particles[nm].dr
		particles[nm].v += particles[nm].dv
		#壁との衝突リスト
		CollisionWallList = CollisionWallDetection( particles[nm] )
		#衝突あり
		if( len(CollisionWallList) > 0 ):

			#２つ以上の壁との衝突がある場合　→　高精度計算モードへ
			if( len(CollisionWallList) > 1 ):
				return False

			#以下、衝突が１箇所
			normal = CollisionWallList[ 0 ]["normal"]

			#いったん巻き戻しておく
			particles[nm].r -= particles[nm].dr
			particles[nm].v -= particles[nm].dv

			#衝突力
			if( CollisionWallList[ 0 ]["type"] == 1 ): #壁との衝突
				#print(nm)
				#相対速度
				barV = - particles[nm].v
				particles[nm].fc = (2.0 * particles[nm].m * np.dot(barV, normal) * normal ) / particles[nm].dt - np.dot(particles[nm].fex, normal) * normal
				#位置と速度を更新
				particles[nm].timeEvolution( t )
				particles[nm].r += particles[nm].dr
				particles[nm].v += particles[nm].dv


	#粒子同士の衝突判定
	for nm1 in range(NM):
		for nm2 in range(NM):
			if( nm1 >= nm2 ): continue
			#２点間ベクトル
			r12 = particles[nm2].r - particles[nm1].r
			r12_length = np.sqrt(np.sum( r12**2 ))

			if( r12_length < particles[nm1].radius + particles[nm2].radius ):
				#print( t, nm1, nm2 ) 
				#すでに衝突力が与えられている場合　→　高精度計算モードへ
				if( (particles[nm1].fc[0] != 0 or particles[nm1].fc[1] != 0 or particles[nm1].fc[2] != 0) 
				 or ( particles[nm2].fc[0] != 0 or particles[nm2].fc[1] != 0 or particles[nm2].fc[2] != 0 )):
					#時間発展：失敗
					return False

				n12 = r12 / r12_length
				v12 = particles[nm2].v - particles[nm1].v
				m1 = particles[nm1].m
				m2 = particles[nm2].m
				barf =  particles[nm2].fex / m2 - particles[nm1].fex / m1 
				f21 = 2.0 * m1 * m2 / ( m1 + m2 ) * np.dot(v12, n12) * n12 / particles[nm1].dt
				f21 += m1 * m2 / ( m1 + m2 ) * np.dot(barf, n12) * n12
				f12 = -f21

				#衝突力の設定
				particles[nm1].fc = f21
				particles[nm2].fc = f12
				#いったん巻き戻しておく
				particles[nm1].r -= particles[nm1].dr
				particles[nm1].v -= particles[nm1].dv
				particles[nm2].r -= particles[nm2].dr
				particles[nm2].v -= particles[nm2].dv
				#時間発展を計算
				particles[nm1].timeEvolution( t )
				particles[nm2].timeEvolution( t )
				#位置と速度を更新
				particles[nm1].r += particles[nm1].dr
				particles[nm1].v += particles[nm1].dv
				particles[nm2].r += particles[nm2].dr
				particles[nm2].v += particles[nm2].dv

				#２点間ベクトル
				r12 = particles[nm2].r - particles[nm1].r
				r12_length = np.sqrt(np.sum( r12**2 ))
				#計算後の２粒子が重なっている場合は失敗
				if( r12_length < particles[nm1].radius + particles[nm2].radius ):
					#時間発展：失敗
					return False

	#計算後、壁とあるいは粒子同士が重なっている場合は失敗
	for nm1 in range(NM):
		#壁との衝突リスト
		CollisionWallList = CollisionWallDetection( particles[nm1] )
		#衝突あり
		if( len(CollisionWallList) > 0 ):
			#時間発展：失敗
			return False
		for nm2 in range(NM):
			if( nm1 >= nm2 ): continue
			#２点間ベクトル
			r12 = particles[nm2].r - particles[nm1].r
			r12_length = np.sqrt(np.sum( r12**2 ))
			if( r12_length < particles[nm1].radius + particles[nm2].radius ):
				#時間発展：失敗
				return False


	#時間発展：正常
	return True


#高精度計算モード
def highPrecisionMode(t, level):
	print ( "レベル" , level, "t=", round(t, level+3 ))

	timeScale = 10**level
	for nm in range(NM):
		particles[nm].r = particles[nm]._r.copy()
		particles[nm].v = particles[nm]._v.copy()

	for _n in range(0, 10):
		t_level = t + dt / timeScale * _n

		result_level = timeEvolution( t_level, timeScale )

		if( result_level == False ): 
			#もっと高精度が必要
			level += 1
			highPrecisionMode(t_level, level)
			level -= 1
			
		else:
			#時間発展前の位置と速度を保持
			for nm in range(NM):
				particles[nm]._r = particles[nm].r.copy()
				particles[nm]._v = particles[nm].v.copy()		


#################################################################
## 計算開始
#################################################################
#粒子
particles = []

for nm in range(NM):

	radius = radius_min + ( radius_max - radius_min) * random.random()
	m = m_min + ( m_max - m_min) * random.random()

	theta_v = 2.0 * np.pi * random.random()
	phi_v = 0 * 2.0 * np.pi * random.random()
	_v_max = v_max * random.random()
	v0 = np.array([
		_v_max * np.sin( theta_v ) * np.cos( phi_v ), 
		_v_max * np.sin( theta_v ) * np.sin( phi_v ), 
		_v_max * np.cos( theta_v )
	])

	for i in range(1000):
		r0 = np.array([
			radius + (L - 2.0 * radius) * random.random(), 
			L / 2.0, 
			radius + (L - 2.0 * radius) * random.random()
		])
		flag = True
		for p in particles:
			#２点間ベクトル
			r12 = r0 - p.r
			r12_length =  np.sqrt(np.sum( r12**2 ))
			if( r12_length < radius + p.radius ):
				flag = False
				break
		if(flag == False): continue

		particles.append( RK4( m, dt, r0, v0 ))
		particles[ len( particles ) - 1 ].radius = radius
		break


#粒子外力の設定
for nm in range(NM):
	particles[nm].fex = fex

####################################
## 描画用関数
####################################
#アニメーション作成用
ims=[]
#ピストン位置の配列
zs = []

def plot( t ):
	#各コマを描画
	for nm in range(NM):
		markersize = particles[nm].radius * 1000
		colorRG = 1.0 - particles[nm].m / m_max
		if( nm == 0 ): img  = plt.plot([particles[nm].r[0]], [particles[nm].r[2]], color = (colorRG, colorRG, 1.0), marker = 'o', markersize = markersize, linestyle='solid', linewidth = 0 )
		else: img += plt.plot([particles[nm].r[0]], [particles[nm].r[2]], color = (colorRG, colorRG, 1.0), marker = 'o', markersize = markersize, linestyle='solid', linewidth = 0 )
	time = plt.text(0.9, 1.12, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )
	E = T() + U()
	print( "t = {:.4f}".format(t)," E =",  E )


####################################
## 各時刻における運動方程式の計算
####################################
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t)

	#時間発展前の位置と速度を保持
	for nm in range(NM):
		particles[nm]._r = particles[nm].r.copy()
		particles[nm]._v = particles[nm].v.copy()		

	level = 0

	#時間発展の計算
	result = timeEvolution( t, 1 )
	
	#複数衝突を検知した場合
	if( result == False ): 
		level += 1
		highPrecisionMode(t, level )



#################################################################
## 分子運動アニメーション
#################################################################
plt.title( u"箱の中の粒子の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", 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.1, 1.1])
plt.ylim([-0.1, 1.1])


#箱の概形
plt.vlines([0-lw], -2*lw, L + 2*lw, "black", linestyles='solid', linewidth = 10 )
plt.vlines([L+lw], -2*lw, L + 2*lw, "black", linestyles='solid', linewidth = 10 )
plt.hlines([-lw], 0, L, "black", linestyles='solid', linewidth = 10 )
plt.hlines([L+lw], 0, L, "black", linestyles='solid', linewidth = 10 )

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output100.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
plt.show()

2021年9月3日2021年9月3日

【統計力学001】理想気体で状態方程式シミュレーション！【Pythonで分子動力学入門①】

############################################################################
#　シリンダーの中のN個の粒子
############################################################################
import os
import numpy as np
import random 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
random.seed(0)

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
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']
#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
M = 1000.0
#粒子数
N = 100
#箱の長さ
L = 1.0
#シリンダーに加わる外力：重力
Fex = np.array([0.0, 0.0, -10.0 * M ])
#シリンダーの粘性係数
beta = 0
#粒子に加わる力
fex = np.array([0.0, 0.0, 0.0 ])
#初期位置と速度
v_max = 10.0

#粒子に加わる力
def F (self, t, r_t, v_t ):
	return self.fex + self.fc - self.m * self.beta * v_t

#時間区間
t_min = 0
t_max = 2.0
#時間区間数
NT =   2000
skip =    5
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#
lw = 0.01
######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		#現時刻の位置と速度
		self.r = r0.copy()
		self.v = v0.copy()
		#１ステップ前の時刻の位置と速度
		self._r = r0.copy()
		self._v = v0.copy()
		#差分
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		#時間間隔
		self.dt = dt
		#外力
		self.fex = np.array([0.0, 0.0, 0.0])
		#衝突力
		self.fc = np.array([0.0, 0.0, 0.0])
		#粘性係数
		self.beta = 0
		#内部
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F(self, t, r_t, v_t ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 時間発展の計算
#################################################################

#衝突判定（戻値：法線ベクトル）
def CollisionDetection( r, R ):
	#衝突リスト
	collisionlist = []
	if(r[0] > L):
		collisionlist.append({ "type": 1, "normal" : np.array([-1.0, 0.0, 0.0]) })
	if(r[0] < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([1.0, 0.0, 0.0]) })
	if(r[1] > L):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, -1.0, 0.0]) })
	if(r[1] < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, 1.0, 0.0]) })
	if(r[2] > R[2]):
		collisionlist.append({ "type": 2, "normal" : np.array([0.0, 0.0, -1.0]) })
	if(r[2] < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, 0.0, 1.0]) })
	return collisionlist

#運動エネルギー
def T():
	_T = 1.0/2.0 * M * np.sum( rk4_M.v**2 )
	for n in range(N):
		_T += 1.0/2.0 * m * np.sum( rk4_ms[n].v**2 )
	return  _T
#位置エネルギー
def U():
	_R = - np.dot(rk4_M.fex, rk4_M.r)
	for n in range(N):
		_R += - np.dot(rk4_ms[n].fex, rk4_ms[n].r)
	return _R

#粒子の運動エネルギー
def T_particle():
	_T = 0
	for n in range(N):
		_T += 1.0/2.0 * m * rk4_ms[n].v[2]**2
	return  _T
#シリンダーの位置エネルギー（底を基準）
def U_cylinder():
	_R = - np.dot(rk4_M.fex, rk4_M.r) + M * L/2 
	return _R

#系の時間発展の計算
def timeEvolution( t, timeScale ):
	#時間スケールの指定
	rk4_M.dt = dt / timeScale

	rk4_M.fc = np.array([0.0, 0.0, 0.0])
	rk4_M.timeEvolution( t )
	rk4_M.r += rk4_M.dr
	rk4_M.v += rk4_M.dv

	#複数衝突
	cylinderCollisionList = []

	for mn in updateList:
		#時間スケールの指定
		rk4_ms[mn].dt = dt / timeScale

		#ルンゲ・クッタ法による時間発展
		rk4_ms[mn].fc = np.array([0.0, 0.0, 0.0])
		rk4_ms[mn].timeEvolution( t )
		#位置と速度を更新
		rk4_ms[mn].r += rk4_ms[mn].dr
		rk4_ms[mn].v += rk4_ms[mn].dv
		#衝突リスト
		collisionlist = CollisionDetection( rk4_ms[mn].r, rk4_M.r )
		#衝突あり
		if( len(collisionlist) > 0 ):

			#２つ以上の壁orシリンダーとの衝突がある場合　→　高精度計算モードへ
			if( len(collisionlist) > 1 ):
				return False

			#以下、衝突が１箇所
			normal = collisionlist[ 0 ]["normal"]

			#いったん巻き戻しておく
			rk4_ms[mn].r -= rk4_ms[mn].dr
			rk4_ms[mn].v -= rk4_ms[mn].dv

			#衝突力
			if( collisionlist[ 0 ]["type"] == 1 ): #壁との衝突
				#相対速度
				barV = - rk4_ms[mn].v
				rk4_ms[mn].fc = (2.0 * m * np.dot(barV, normal) * normal ) / rk4_ms[mn].dt - np.dot(rk4_ms[mn].fex, normal) * normal
				#位置と速度を更新
				rk4_ms[mn].timeEvolution( t )
				rk4_ms[mn].r += rk4_ms[mn].dr
				rk4_ms[mn].v += rk4_ms[mn].dv

			if( collisionlist[ 0 ]["type"] == 2 ): #ピストンとの衝突
				cylinderCollisionList.append( { "mn": mn, "normal":normal}  )


	#ピストンとの衝突ありの場合
	if( len(cylinderCollisionList) > 0 ):
		#ピストンをいったん巻き戻して
		rk4_M.r -= rk4_M.dr
		rk4_M.v -= rk4_M.dv
		#衝突数が１個の場合
		if( len(cylinderCollisionList) == 1 ):
			mn = cylinderCollisionList[0]["mn"]
			normal = cylinderCollisionList[0]["normal"]
			#衝突力の計算
			barV = rk4_M.v - rk4_ms[mn].v
			barFex =  rk4_M.fex / M - rk4_ms[mn].fex / m 
			rk4_ms[mn].fc = 2.0 * m * M / ( m + M ) * np.dot(barV, normal) * normal / rk4_ms[mn].dt
			rk4_ms[mn].fc += m * M / ( m + M ) * np.dot(barFex, normal) * normal

			rk4_M.fc = - rk4_ms[mn].fc
			#時間発展を計算
			rk4_ms[mn].timeEvolution( t )
			rk4_M.timeEvolution( t )
			#位置と速度を更新
			rk4_M.r += rk4_M.dr
			rk4_M.v += rk4_M.dv
			rk4_ms[mn].r += rk4_ms[mn].dr
			rk4_ms[mn].v += rk4_ms[mn].dv

		#衝突数が２個以上の場合
		if(len(cylinderCollisionList) > 1):
			#時間発展：失敗
			return False

	#時間発展：正常
	return True

#高精度計算モード
def highPrecisionMode(t, level):
	print ( "レベル" , level, "t=", round(t, level+3 ))

	timeScale = 10**level
	#時間ステップ計算前に巻き戻す
	rk4_M.r = rk4_M._r.copy()
	rk4_M.v = rk4_M._v.copy()
	for mn in range(N):
		rk4_ms[mn].r = rk4_ms[mn]._r.copy()
		rk4_ms[mn].v = rk4_ms[mn]._v.copy()

	for _n in range(0, 10):
		t_level = t + dt / timeScale * _n

		result_level = timeEvolution( t_level, timeScale )

		if( result_level == False ): 
			#もっと高精度が必要
			level += 1
			highPrecisionMode(t_level, level)
			level -= 1
			
		else:
			#時間発展前の位置と速度を保持
			rk4_M._r = rk4_M.r.copy()
			rk4_M._v = rk4_M.v.copy()
			for mn in updateList:
				rk4_ms[mn]._r = rk4_ms[mn].r.copy()
				rk4_ms[mn]._v = rk4_ms[mn].v.copy()		


#################################################################
## 計算開始
#################################################################
#初速度配列
v0s = []
#初期運動エネルギーの計算（z成分のみ）
Tz0 = 0
for n in range(N):
	theta_v = np.pi * random.random()
	phi_v = 2.0 * np.pi * random.random()
	_v_max = v_max * random.random()
	v0s.append( np.array([
		_v_max * np.sin( theta_v ) * np.cos( phi_v ), 
		_v_max * np.sin( theta_v ) * np.sin( phi_v ), 
		_v_max * np.cos( phi_v )
	]))
	Tz0 += 1.0/2.0 * m * v0s[len(v0s)-1][2]**2

#ピストンの初期位置
h0 = L #2*Tz0/(M*10)

rk4_ms = []
for n in range(N):
	r0 = np.array([
		L * random.random(), 
		L * random.random(), 
		h0 * random.random()
	])
	rk4_ms.append( RK4(m, dt, r0, v0s[n]) )
#粒子外力の設定
for n in range(N):
	rk4_ms[n].fex = fex


#シリンダーの初期値

R0 = np.array([0.0, 0.0, h0  ])
V0 = np.array([0.0, 0.0, 0.0])
rk4_M = RK4(M, dt, R0, V0)
rk4_M.beta = beta
rk4_M.fex = Fex

####################################
## 描画用関数
####################################
#アニメーション作成用
ims=[]
#ピストン位置の配列
zs = []

def plot( t ):
	#各コマを描画
	img  = plt.plot([0-lw, L+lw], [rk4_M.r[2] + lw, rk4_M.r[2] + lw], colors[1], linestyle='solid', linewidth = 10 )
	for n in range(N):
		img += plt.plot([rk4_ms[n].r[0]], [rk4_ms[n].r[2]], colors[0], marker = 'o', markersize = 10, linestyle='solid', linewidth = 0 )
	time = plt.text(0.9, 1.12, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )
	E = T() + U()
	print( "{:.2f}".format(t), E )
	zs.append( [t, rk4_M.r[2]] )
	#print( "{:.2f}".format(t), U_cylinder(), T_particle(), 3.0 * U_cylinder()/T_particle() )


####################################
## 各時刻における運動方程式の計算
####################################
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t)

	#計算対象粒子リスト（ここでは全部）
	updateList = list(range(0, N))

	#時間発展前の位置と速度を保持
	rk4_M._r = rk4_M.r.copy()
	rk4_M._v = rk4_M.v.copy()
	for mn in updateList:
		rk4_ms[mn]._r = rk4_ms[mn].r.copy()
		rk4_ms[mn]._v = rk4_ms[mn].v.copy()		

	level = 0

	#時間発展の計算
	result = timeEvolution( t, 1 )

	#複数衝突を検知した場合
	if( result == False ): 
		level += 1
		highPrecisionMode(t, level )

#フォルダ指定
dir_path = "output/"
#フォルダ生成
os.makedirs(dir_path, exist_ok = True)
filepath = dir_path + "z.txt"
#ファイルへ保存
np.savetxt( filepath, zs)

#################################################################
## 分子運動アニメーション
#################################################################
plt.title( u"シリンダーの中の粒子の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", 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.1, 1.1])
plt.ylim([-0.1, 1.1])


#箱の概形
plt.vlines([0-lw], -2*lw, L, "black", linestyles='solid', linewidth = 10 )
plt.vlines([L+lw], -2*lw, L, "black", linestyles='solid', linewidth = 10 )
plt.hlines([-lw], 0, L, "black", linestyles='solid', linewidth = 10 )

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
ani.save("output.gif", writer="imagemagick")
#ani.save("output100.mp4", writer="ffmpeg", dpi=300)

#グラフの表示
#plt.show()

#########################################################################
# グラフの描画（エネルギー固有値）
#########################################################################
fig1 = plt.figure(figsize=(15, 8))
plt.title( u"ピストンの位置", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel( u"時間" + r"$\,[{\rm s}]$", fontsize=30, fontname="Yu Gothic", fontweight=1000  )
plt.ylabel( u"高さ" + r"$\,[{\rm m}]$", fontsize=30, fontname="Yu Gothic", fontweight=1000 )
#余白の調整
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([0, t_max])
plt.ylim([0, L])

xs = []
ys = []
for zt in zs:
	xs.append( zt[0] )
	ys.append( zt[1] )

#グラフを描画	
plt.plot(xs, ys, linewidth = 5)

#グラフの表示
plt.show()

2021年9月1日2021年9月1日

【第28話】量子ビットの作り方５：１量子ビットの完成！！【量子ドットコンピュータの原理⑤】

############################################################################
#　衝立ありの無限に深い量子井戸に静電場を加えた電子状態
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import scipy.integrate as integrate
import numpy as np
import numpy.linalg as LA

#全体設定
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.6260896 * 10**-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.10938215 * 10**-31
#電子ボルト
eV = 1.60217733 * 10**-19
#電気素量
e = 1.60217733 * 10**-19

#複素数
I = 0.0 + 1.0j
######################################
#　物理系の設定
######################################
#量子井戸の幅
L = 1.0 * 10**-9
#計算区間
x_min = -L / 2.0
x_max = L / 2.0
#状態数
n_max = 30
#行列の要素数
DIM = n_max + 1
#空間分割数
NX = 500
#空間刻み間隔
nm = 1.0 * 10**-9
#壁の厚さ
W = L / 5
#壁の高さの最大値
V_max = 30.0 * eV

#電場の強さ
Ex_max = 1.0 * 10**8
#電場強度の分割数
NEx = 10

#基底状態と励起状態
N = 2

######################################
#　衝立なし・静電場無しの解析解
######################################
#固有関数
def varphi(n, x):
	kn = np.pi * (n + 1) / L
	return np.sqrt(2.0 / L) * np.sin(kn * (x + L / 2.0))

#ポテンシャル項
def V(x, Ex):
	if(abs(x) <= W / 2.0): 
		return e * Ex * x + V_max
	else: 
		return e * Ex * x

#固有エネルギー
def Energy(n):
	kn = np.pi * (n + 1) / L
	return hbar * hbar * kn * kn / (2.0 * me)

######################################
#　静電場ごとの固有関数の計算
######################################
#被積分関数（行列要素計算用）
def integral_matrixElement(x, n1, n2, Ex):
	return varphi(n1 ,x) * V(x, Ex) * varphi(n2, x) / eV

#被積分関数（平均計算用）
def average_x(x, a):
	sum = 0
	for n in range(n_max + 1):
		sum += a[n] * varphi(n, x)
	return x * sum**2

#固有値・固有ベクトルの初期化
eigenvalues = [0] * (NEx + 1)
vectors = [0] * (NEx + 1)
for nEx in range(NEx + 1):
	eigenvalues[nEx] = []
	vectors[nEx] = []

#存在確率分布グラフ描画用の配列初期化
xs = []
phi = [0] * (NEx + 1)
for nEx in range(NEx + 1):
	phi[nEx] = [0] * N
	for n in range( len(phi[nEx]) ):
		phi[nEx][n] = [0] * (NX + 1)

#中心の電場依存性グラフ描画用の配列初期化
averageX = [0] * N
for n in range(len(averageX)):
	averageX[n] = [0] * (NEx + 1)


#静電場強度ごとに
for nEx in range(NEx + 1):
	print("電場強度：" + str( nEx * 100 / NEx ) + "%")
	#静電場の強度を設定
	Ex = Ex_max / NEx * nEx
	#エルミート行列（リスト）
	matrix = []

	###行列要素の計算
	for n1 in range(n_max + 1):
		col=[]
		for n2 in range(n_max + 1):
			#ガウス・ルジャンドル積分
			result = integrate.quad(
				integral_matrixElement, #被積分関数
				x_min, x_max,		    #積分区間の下端と上端
				args=(n1, n2, Ex)		#被積分関数へ渡す引数
			)
			real = result[0]
			imag = 0j
			#無静電場のエネルギー固有値（対角成分）
			En = Energy(n1)/eV if (n1 == n2) else 0
			#行の要素を追加
			col.append( En + real )
		#行を追加
		matrix.append( col )

	#リスト → 行列
	matrix = np.array( matrix )

	###固有値と固有ベクトルの計算
	result = LA.eig( matrix )
	eig = result[0] #固有値
	vec = result[1] #固有ベクトル

	#小さい順に並べるためのインデックス（配列）
	index = np.argsort( eig )
	#固有値を小さい順に並び替え
	eigenvalues[nEx] = eig[ index ]

	#転置行列
	vec = vec.T
	#固有ベクトルの並び替え
	vectors[nEx] = vec[ index ]

	### 検算：MA-EA=0 ?
	sum = 0
	for i in range(DIM):
		v = matrix @ vectors[nEx][i] - eigenvalues[nEx][i] * vectors[nEx][i]
		for j in range(DIM):
			sum += abs(v[j])**2
	print("|MA-EA| =" + str(sum))

	###固有関数の空間分布
	for nx in range(NX+1):
		x = x_min + (x_max - x_min) / NX * nx
		if(nEx == 0): xs.append( x/nm )

		for n in range( len(phi[nEx]) ):
			for m in range(n_max+1):
				phi[nEx][n][nx] += vectors[nEx][n][m] * varphi(m, x)

			#描画用データの整形
			phi[nEx][n][nx] = abs(phi[nEx][n][nx])**2 / (1.0 * 10**9)


	for n in range(len(averageX)):
		#ガウス・ルジャンドル積分
		result = integrate.quad(
			average_x,            #被積分関数
			x_min, x_max,		      #積分区間の下端と上端
			args=(vectors[nEx][n])		#被積分関数へ渡す引数
		)
		#計算結果の取得
		averageX[n][nEx] = result[0] * (1.0 * 10**9)


#########################################################################
# グラフの描画（エネルギー固有値）
#########################################################################
fig1 = plt.figure(figsize=(15, 8))
plt.title( u"無限に深い量子井戸（衝立あり）に静電場を加えたときの固有エネルギー", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel( u"静電場の強さ" + r"$\times 10^7 \,[{\rm V/m}]$", fontsize=30, fontname="Yu Gothic", fontweight=1000  )
plt.ylabel( r"$ E\,[{\rm eV}]$", fontsize=30, fontname="Yu Gothic", fontweight=1000 )
#余白の調整
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([0, 10])
#x軸
exs = range( NEx + 1)
#y軸
En_0 = []
En_1 = []
for nEx in range(NEx + 1):
	En_0.append( eigenvalues[nEx][0] )
	En_1.append( eigenvalues[nEx][1] )
	#print( str(nV) + " " + str( eigenvalues[nV][0] ) + " " + str( eigenvalues[nV][1] ))
#基底状態と第１励起状態のグラフを描画	
plt.plot(exs, En_0, marker="o", markersize=15, linewidth = 5)
plt.plot(exs, En_1, marker="o", markersize=15, linewidth = 5)

#グラフの表示
plt.show()

#########################################################################
# グラフの描画（基底状態）
#########################################################################
fig2 = plt.figure(figsize=(15, 8))
plt.title( u"無限に深い量子井戸（衝立あり）に静電場を加えたときの基底状態", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel( r"$x\, [{\rm nm}]$", fontsize=30 )
plt.ylabel( r"$|\varphi^{(m)}(x)|^2$", 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], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.vlines([0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0],  -0.5, 0.5, "black", linestyles='solid', linewidth = 10 )
#衝立の概形
plt.vlines([-0.1], -1, 3, "black", linestyles='solid', linewidth = 10 )
plt.vlines([0.1], -1, 3, "black", linestyles='solid', linewidth = 10 )
plt.hlines([3],  -0.106, 0.106, "black", linestyles='solid', linewidth = 10 )

#描画範囲を設定
plt.xlim([-0.6, 0.6])
plt.ylim([0, 5.0])

#各
for nEx in range(NEx + 1):
	plt.plot(xs, phi[nEx][0] , linewidth = 5)

#グラフの表示
plt.show()

#########################################################################
# グラフの描画（第１励起状態）
#########################################################################
fig3 = plt.figure(figsize=(15, 8))
plt.title( u"無限に深い量子井戸（衝立あり）に静電場を加えたときの第１励起状態", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel( r"$x\, [{\rm nm}]$", fontsize=30 )
plt.ylabel( r"$ |\varphi^{(m)}(x)|^2$", 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([-0.6, 0.6])
plt.ylim([0, 5.0])

#井戸の概形
plt.vlines([-0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.vlines([0.5], -1, 20, "black", linestyles='solid', linewidth = 10 )
plt.hlines([0],  -0.5, 0.5, "black", linestyles='solid', linewidth = 10 )
#衝立の概形
plt.vlines([-0.1], -1, 3, "black", linestyles='solid', linewidth = 10 )
plt.vlines([0.1], -1, 3, "black", linestyles='solid', linewidth = 10 )
plt.hlines([3],  -0.106, 0.106, "black", linestyles='solid', linewidth = 10 )

for nEx in range(NEx + 1):
	plt.plot(xs, phi[nEx][1] , linewidth = 5)	

#グラフの表示
plt.show()