今回はプログラミング言語の「Python」を使って作成しました。
下記の本のサンプルの「分子運動のシミュレーション」をPygameで動くように移植しました。
下画像の式で動いています。
■ プログラム
import traceback
try:
import sys
import math
import random
import pygame
from pygame.locals import QUIT, KEYDOWN, K_LEFT, K_RIGHT, K_SPACE, Rect
#-----------------------------------------------------------
# クラス定義
#-----------------------------------------------------------
class Block:
""" ボールオブジェクト """
def __init__(self, col, rect, W):
#--- 初期値設定 ---
self.col = col
self.rect = rect
self.W = W
self.x = rect.x
self.y = rect.y
def move(self):
""" ボールを動かす """
self.rect = Rect(self.x,self.y, self.W, self.W)
def draw(self):
""" ボールを描画する """
pygame.draw.ellipse(SURFACE, self.col, self.rect)
def array(N1, N2=0, N3=0):
if N2 == 0 and N3==0: return [0 for i in range(N1)]
elif N3 == 0: return [[0 for i in range(N2)] for j in range(N1)]
else: return [[[0 for i in range(N3)] for j in range(N2)]
for k in range(N1)]
def genMole():
global N, R
for i in range(N):
R[0][i][0] = random.random() * (Xmax - 2.0) + 1.0
R[0][i][1] = random.random() * (Ymax - 2.0) + 1.0
Xmax = 10; Ymax = 10; DispR = 20
EPSI = 1.0; SIGMA = 1;MG = 1
N=40 # 分子の数
DT=0.05 # 時間刻み
timerInt=10 # 表示時間間隔
R=array(2,N,2); F=array(N,2) # R:位置，F:力
V=array(2,N,2); Bname=array(N) # V:速度, Bname:表示円形
E1 = -48 * EPSI * pow(SIGMA,13) # 計算用係数
E2 = 24 * EPSI * pow(SIGMA, 7)
alfa= 0.1 # alfa = M/(3kT) とみなす(M: 分子質量合計)
ID1 = 0; ID2 = 1 # 実行制御用
Nup=0;Vup=0;Nlf=0; Vlf=0 #衝突回数と衝突速度の絶対値合計
Nbt=0;Vbt=0;Nrt=0; Vrt=0
sCtr=0 #シミュレーション回数
#--------------------------------------------------------------
# 画面上の分子の生成
#--------------------------------------------------------------
def genMono(canvas, X0, Y0,R):
return canvas.create_oval(X0-R,Y0-R, X0+R, Y0+R,fill="red")
def init(): #■初期設定
global N, F, V,Bname,DispR
for i in range(N):
Bname[i] = Block((random.randint(0, 255), random.randint(0, 255), random.randint(0, 255)), Rect(random.randint(0, 150), random.randint(-150, 0), DispR, DispR), DispR)
#--------------------------------------------------------------
# 画面上の分子の表示
#--------------------------------------------------------------
def dsp(ID1):
global N, canvas, R, Bname, DispR
RD=DispR
for i in range(N):
X = R[ID1][i][0]*20
Y = R[ID1][i][1]*20
""" ボールを描画する """
Bname[i].rect.x = X-RD
Bname[i].rect.y = Y-RD
Bname[i].rect.w = RD
Bname[i].rect.h = RD
pygame.draw.ellipse(SURFACE, Bname[i].col, Bname[i].rect)
def force(ID1): #■力の計算
global N, F, R, alfa
for i in range(N):
F[i][0] = 0; F[i][1] = 0
for j in range(N):
if i != j :
dX = R[ID1][j][0] - R[ID1][i][0]
dY = R[ID1][j][1] - R[ID1][i][1]
DR = math.sqrt(dX * dX + dY * dY)
if DR < 1 : DR = 1
FF = E1 / math.pow(DR,13) + E2 / math.pow(DR,7)
F[i][0] = F[i][0] + FF * dX / DR
F[i][1] = F[i][1] + FF * dY / DR
def velocity(ID1,ID2): #■速度の計算
global N, F, V,MG
for i in range(N):
V[ID2][i][0] = V[ID1][i][0] + DT * F[i][0]/ MG
V[ID2][i][1] = V[ID1][i][1] + DT * F[i][1]/ MG
T = 0 # 温度境界条件
for i in range(N):
VX = V[ID2][i][0]; VY = V[ID2][i][1]
T += VX * VX + VY * VY
T = math.sqrt(T)* alfa
if T > 0.00001:
for i in range(N):
V[ID2][i][0] /= T; V[ID2][i][1] /= T
def place(ID1,ID2): #■位置の計算
global N, V, X
for i in range(N):
R[ID2][i][0] = R[ID1][i][0] + DT * V[ID2][i][0]
R[ID2][i][1] = R[ID1][i][1] + DT * V[ID2][i][1]
def wallHit(ID2):#■壁への衝突
global Nup,Vup,Nlf, Vlf, Nbt, Vbt,Nrt, Vrt
for i in range(N):
if R[ID2][i][0] < (0.5 * SIGMA):
R[ID2][i][0] = 0.5 * SIGMA
AV=abs(V[ID2][i][0])
V[ID2][i][0] = AV; Nlf+=1; Vlf+=AV
if R[ID2][i][0] > (Xmax - 0.5 * SIGMA):
R[ID2][i][0] = Xmax - 0.5 * SIGMA
AV=abs(V[ID2][i][0])
V[ID2][i][0] = -AV; Nrt+=1; Vrt+=AV
if R[ID2][i][1] < (0.5 * SIGMA):
R[ID2][i][1] = 0.5 * SIGMA
AV=abs(V[ID2][i][1])
V[ID2][i][1] = AV; Nup+=1; Vup+=AV
if R[ID2][i][1] > (Ymax - 0.5 * SIGMA):
R[ID2][i][1] = Ymax - 0.5 * SIGMA
AV=abs(V[ID2][i][1])
V[ID2][i][1] = -AV; Nbt+=1; Vbt+=AV
#--------------------------------------------------------------
# 画面上の分子のシミュレーション
#--------------------------------------------------------------
def simulate(ID1, ID2):
global sCtr,Vup,Vlf,Vbt, Vrt, root
dsp(ID1); force(ID1); velocity(ID1,ID2)
place(ID1,ID2); wallHit(ID2)
sCtr+=1;
#print("%d, %7.4f, %7.4f, %7.4f, %7.4f" % (sCtr, Vup/sCtr,Vbt/sCtr, Vlf/sCtr,Vrt/sCtr))
#--------------------------------------------------------------
# メイン処理
#--------------------------------------------------------------
init(); genMole();ID1=0
ID1=0
# Pygame初期設定
pygame.init()
SURFACE = pygame.display.set_mode((200, 200))
FPSCLOCK = pygame.time.Clock()
fps = 60
#--------------------------------------------------------------
# ループ処理
#--------------------------------------------------------------
while True:
#-----------------------------------------------------------
# 描画処理
#-----------------------------------------------------------
SURFACE.fill((0, 0, 0))
ID2=1-ID1;
simulate(ID1,ID2)
ID1=ID2
pygame.display.update()
FPSCLOCK.tick(fps)
except Exception as e:
print("エラー情報\n" + traceback.format_exc())
input()
■ 参考書
「Python 3.6 による 力学シミュレーション 第４巻: 分子」
■ ゲーム用ライブラリ
「Pygame」
■ プログラミング言語
「Python」
