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

############################################################################
#　[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()
コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル
