今回はプログラミング言語の「Python」を使って作成しました。
下記の本のサンプルの「炎のシミュレーション」をPygameで動くように移植しました。
■ プログラム
import traceback
try:
import random, math, time
import sys
import pygame
from pygame.locals import QUIT
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)]
bitmap=array(200,200)
Xmax = 5; Ymax = 10; MG=1; #DispR = 1
EPSI = 0.5;SIGMA = 0.25;Sig2 = 0.5 * SIGMA
XSig2 = Xmax - Sig2;YSig2 = Ymax - Sig2
Xmax2=Xmax/2
N=2000 # 粒子の数
NCof=math.sqrt(2000/N) # 粒子数の正規化比(2000を標準とする)
DT=0.2 # 時間刻み
timerInt=10 # 表示時間間隔
R=array(2,N,2); F=array(N,2); V=array(2,N,2) #Bname=array(N)
Btemp=array(N); tFact=0.999
E1 = -48 * EPSI * pow(SIGMA,13)
E2 = 24 * EPSI * pow(SIGMA, 7)
# Pygame初期設定
pygame.init()
SURFACE = pygame.display.set_mode((200, 200))
FPSCLOCK = pygame.time.Clock()
fps = 60
src = pygame.Surface((200, 200))
data = pygame.PixelArray(src)
#■分子生成
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
#■表示
Nrend=17; NrendC=int((Nrend-1)/2);NrendC2=NrendC*NrendC
#■レンダリング用テーブル作成
def setRendTab(): #表示中に計算するのは効率が悪いので
global Nrend, NrendC #この表を用意しておく
A=array(Nrend,Nrend)
for i in range(Nrend):
x=i-NrendC
for j in range(Nrend):
y=j-NrendC
r=math.sqrt(x*x+y*y)
if r<NrendC:
h=r-NrendC; A[i][j]=h*h/NrendC2
return A
#■値の範囲チェック
def checkRange(D, D1=0, D2=255):
if D<D1: return D1
if D>D2: return D2
return int(D)
#■表示
def dsp(ID1):
global N, canvas, R, rendTab,Nrend, NrendC,bitmap, RGBFlag, data
Tmap=array(200,200)
#各セルのエネルギー合計を求める
Xmn=200; Xmx=0; Ymn=200; Ymx=0
for i in range(N):
X = int(R[ID1][i][0]*20+50)
Y = int(R[ID1][i][1]*20)
if Btemp[i]>=0.24: #元々エネルギーが小さい粒子については計算しない
Xmn=min(X,Xmn); Ymn=min(Y,Ymn)
Xmx=max(X,Xmx); Ymx=max(Y,Ymx)
for k1 in range(-NrendC,NrendC+1): #粒子が広がりを持っている
YY = Y + k1 #ものとしてエネルギーの
if YY>=0 and YY<200: #合計値を計算する
i1=k1+NrendC
for k2 in range(-NrendC,NrendC+1):
XX = X + k2
if XX>=0 and XX<200:
Tmap[XX][YY] += rendTab[k2+NrendC][i1]*Btemp[i]
#粒子のレンダリング範囲を設定
Xmn-=NrendC; Ymn-=NrendC; Xmx+=NrendC; Ymx+=NrendC
Xmn=checkRange(Xmn,0,200); Ymn=checkRange(Ymn,0,200)
Xmx=checkRange(Xmx,0,200); Ymx=checkRange(Ymx,0,200)
if Xmn<Xmx and Ymn<Ymx:#範囲矛盾のチェック(安全のため)
for i in range(Xmn,Xmx):#レンダリング
for j in range(Ymn,Ymx):
if Tmap[i][j] >=1.0:
TT=Tmap[i][j]*NCof # 2,000個の合計値を標準とする
if TT<=5:
TT1=TT-1
RR=checkRange(63.75*TT1) # 係数 255/4=63.75
G =checkRange(16*TT) # 係数 64/4=16
B=0
elif TT<=20: # 係数 191/25=7.64
RR=255; G=checkRange(64+7.64*(TT-5));B=0
else:# 係数 255/30=85
RR=255; G=255; B = checkRange(85*(TT-30))
if RR !=0 or G!=0 :
data[i][j] = (RR, G, B)
#■力の計算
def force(ID1):
global N, F, R
for i in range(N):
#以下３行、揺らぎを与えるとき
#if R[ID1][i][1]< 5:F[i][0] = (random.random()-0.8*random.random())*3
#else: F[i][0] = (random.random()-1.0*random.random())*3
#F[i][1] = -random.random() *1.5 #上向きのランダムな力
#以下２行、揺らぎを与えないとき
F[i][0] = (random.random()-0.5)*3
F[i][1] = -random.random() *1.5 #上向きのランダムな力
#■速度計算
def velocity(ID1,ID2):
global N, F, V,MG; DMG=DT/MG
for i in range(N):
V[ID2][i][0] = V[ID1][i][0] + DMG * F[i][0]
V[ID2][i][1] = V[ID1][i][1] + DMG * F[i][1]
T = 0 # 温度境界条件
for i in range(N):
VX = V[ID2][i][0]; VY = V[ID2][i][1]
T += VX * VX + VY * VY
if T > 0.001:
T = math.sqrt(T) * 0.1
for i in range(N):
V[ID2][i][0] /= T; V[ID2][i][1] /= T
#■位置計算
def place(ID1,ID2):
global N, V, X, Btemp,tFact
XX2=Xmax/2
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]
XX=R[ID2][i][0]-XX2 # 水平方向のエネルギー低下が大きいとする
YY=(R[ID2][i][1]-Ymax)*0.05
RR=math.sqrt(XX*XX+YY*YY)
Btemp[i] *= tFact*math.exp(-0.05*RR)
#■壁への衝突
def wallHit(ID2):
global Sig2,XSig2, YSig2, Xmax2
for i in range(N): #上壁（天井）の衝突だけを残す
if R[ID2][i][1] < Sig2: # 煙の粒子が天井に届いたら
R[ID2][i][0]=Xmax2 # 下から発生する粒子として
R[ID2][i][1]=YSig2 # 再利用
Btemp[i] =1.0
#■シミュレーション
def simulate(ID1, ID2):
dsp(ID1);force(ID1); velocity(ID1,ID2)
place(ID1,ID2); wallHit(ID2)
#■初期処理用シミュレーション
def simulate2(ID1, ID2):
force(ID1); velocity(ID1,ID2)
place(ID1,ID2); wallHit(ID2)
#■メイン処理
rendTab=setRendTab()
ID1=0; ID2=1-ID1
#■粒子が連続的に発生するよう初期調整
simulate2(ID1,ID2); ID1=ID2
checkN=int(N/100);k=0; NK=1
for i in range(N):
k+=1
if (k % checkN) ==0:
print("準備中です(%2d %%)" % int(i/N*100))
print("■"*NK)
NK+=1
if NK>20: NK=1
R[ID2][i][0]=Xmax/2;R[ID2][i][1]=Ymax-(0.5*SIGMA)
Btemp[i]=1.0 #表示しないのでtime.sleep(0.001000)は不要
simulate2(ID1,ID2)
simulate(ID1,ID2); ID1=ID2
#■シミュレーション開始
while True:
for _ in pygame.event.get(QUIT):
pygame.quit()
sys.exit()
ID2=1-ID1
#----------------------------------------------------
# dataに炎のRGBを格納する
#----------------------------------------------------
data = pygame.PixelArray(src)
simulate(ID1,ID2)
del data
SURFACE.blit(src, (0, 0))
pygame.display.update()
FPSCLOCK.tick(fps)
ID1=ID2
except Exception as e:
print("エラー情報\n" + traceback.format_exc())
input()
■ 参考書
「Python 3.6 による 力学シミュレーション 第４巻: 分子」
■ ゲーム用ライブラリ
「Pygame」
■ プログラミング言語
「Python」
