第4章 分　子

４．１　分子間力
　
【リスト4-1】分子の初期位置設定例と共通処理
from tkinter import *
from tkinter import messagebox
import random, math, time
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
　
【リスト4-2】シミュレーション・プログラム
・
・(前略) 【リスト4-1】と同じ
・
Xmax = 10; Ymax = 10; DispR = 10
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] = genMono(canvas, 0, 0, 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
        canvas.coords(Bname[i], X-RD, Y-RD, X+RD, Y+RD)
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)) 
    root.update()
#■画面コピー用マウスイベント  
def leftMouseDown(event):
    global ID1, ID2
    ID2=1-ID1; simulate(ID1, ID2); ID1=ID2
#■メイン処理
root=Tk(); root.title("分子運動のシミュレーション")
canvas=Canvas(root, width=200,height=200); canvas.pack()
init(); genMole();ID1=0
#左ボタン
#canvas.bind("<Button-1>",leftMouseDown)  #画面コピー用マウス
#root.mainloop()
#画面コピーするときは以下をコメントにする
ID1=0
print("計算回数,上壁, 下壁,左壁, 右壁")
print("0,  0.0000,  0.0000,  0.0000,  0.0000")
while True:
    ID2=1-ID1; simulate(ID1,ID2)
    time.sleep(0.001)
    ID1=ID2

４．２　液体もどきの分子運動

【リスト4-3】シミュレーション・プログラム
・
・(前略) 【リスト4-1】と同じ
・
Xmax = 40; Ymax = 20; DispR = 15
EPSI = 1; SIGMA = 1;MG = 1; GC=0.1
N=60                               # 分子の数
DT=0.01                            # 時間刻み
・
・(中略) 【リスト4-2】と同じ
・
def genMono(canvas, X0, Y0,R):#■画面上の分子の生成
    return canvas.create_oval(X0-R,Y0-R, X0+R, Y0+R,fill="blue", outline="blue")
・
・(中略) 【リスト4-2】と同じ
・
def dsp(ID1):           #■表示
    global N, canvas, R, Bname, DispR
    RD=DispR
    for i in range(N):
        X = R[ID1][i][0]*10
        Y = R[ID1][i][1]*10
・
・(中略) 【リスト4-2】と同じ
・
def force(ID1):         #■力の計算
・
・(中略) 【リスト4-2】と同じ
・
                F[i][1] = F[i][1] + FF * dY / DR+GC
def velocity(ID1,ID2):  #■速度の計算
・
・(中略) 【リスト4-2】と同じ
・
def simulate(ID1, ID2):#■シミュレーション
    global sCtr,Vup,Vlf,Vbt, Vrt, root
    dsp(ID1); force(ID1); velocity(ID1,ID2)
    place(ID1,ID2); wallHit(ID2)
    root.update()
・
・(中略) 【リスト4-2】と同じ
・
root=Tk()
root.title("分子運動のシミュレーション")
canvas=Canvas(root, width=400,height=200)
canvas.pack()
init(); genMole()
ID1=0
#左ボタン
#canvas.bind("<Button-1>",leftMouseDown)  #画面コピー用マウスイベント
#for i in range(1200):
#    ID2=1-ID1; simulate(ID1,ID2);ID1=ID2
#root.update();root.mainloop()
while True:
    ID2=1-ID1
    simulate(ID1,ID2)
    time.sleep(0.001)
    ID1=ID2

４．３　ブラウン運動

【リスト4-4】 変更部分のみ(他は【リスト4-2】と同じ）
　・
　・(前略)【リスト4-1】と同じ
　・
Xmax = 10; Ymax = 10; DispR = 1 #ほこりっぽく見せるために小さくする
EPSI = 1.0; SIGMA = 1;MG = 1
N=100                              # 分子の数
　・
　・(中略)【リスト4-2】と同じ
　・
def force(ID1):         #■力の計算
    global N, F, R
    for i in range(N): # 力をランダムな力にする
        F[i][0] = random.random()-0.5
        F[i][1] = random.random()-0.5
　・
　・(中略)【リスト4-2】と同じ
　・
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 * 0.0 # 0.0 以上にすると端に集まりにくい
        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 * 0.0 # 0.0 以上にすると端に集まりにくい
        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 * 0.0 # 0.0 以上にすると端に集まりにくい
        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 * 0.0 # 0.0 以上にすると端に集まりにくい

４．４　煙の運動

【リスト4-5】
・
・(中略)【リスト4-1】と同じ
・
Xmax = 10; Ymax = 10; DispR = 0
EPSI = 0.5; SIGMA = 0.25; MG=1
N=500        # 粒子の数
DT=0.2       # 時間刻み　・
・
・(中略)【リスト4-2】と同じ
・
def genMono(canvas, X0, Y0,R):#■画面上の分子の生成
    return canvas.create_oval(X0-R,Y0-R, X0+R, Y0+R,fill="gray", outline="gray")
・
・(中略)【リスト4-2】と同じ
・
def force(ID1):
    global N, F, R
    for i in range(N):
        if R[ID1][i][1]< 5:
            F[i][0] = (random.random()-0.4)*3
        else:
            F[i][0] = (random.random()-0.6)*3
        F[i][1] = -random.random() *1.5 #上向きのランダムな力
・
・(中略)【リスト4-2】と同じ
・
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 * 0.0 # 0.0 以上にすると端に集まりにくい
        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 * 0.0 # 0.0 以上にすると端に集まりにくい
        if R[ID2][i][1] < (0.5 * SIGMA):  #　煙の粒子が天井に届いたら
            R[ID2][i][0]=Xmax/2           #　下から発生する粒子として
            R[ID2][i][1]=Ymax-(0.5*SIGMA) #　再利用
        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 * 0.0 # 0.0 以上にすると端に集まりにくい
・
・(中略)【リスト4-2】と同じ
・
#■メイン処理
root=Tk(); root.title("煙の運動")
canvas=Canvas(root, width=200,height=200); canvas.pack()
init(); genMole()
ID1=0; ID2=1-ID1
#■煙の粒子が連続的に発生するよう初期調整
simulate(ID1,ID2); ID1=ID2
for i in range(N):
    R[ID2][i][0]=Xmax/2; R[ID2][i][1]=Ymax-(0.5*SIGMA)
    time.sleep(0.001); simulate(ID1,ID2)
#左ボタン
#canvas.bind("<Button-1>",leftMouseDown)  #画面コピー用マウス
#画面コピーするときは以下をコメントにする
while True:
    ID2=1-ID1;simulate(ID1,ID2); time.sleep(0.001); ID1=ID2

４．５　炎

【リスト4-6】
from tkinter import *
from tkinter import messagebox
import random, math, time
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)
#■分子生成
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
    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 :
                        C="#%02X%02X%02X" % (RR, G,B)
                        if C != canvas.itemcget(bitmap[i][j], option='fill'):
                            canvas.itemconfig(bitmap[i][j], fill=C)
    root.update()
#■力の計算
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)
    root.update()
#■初期処理用シミュレーション
def simulate2(ID1, ID2):
    force(ID1); velocity(ID1,ID2)
    place(ID1,ID2);  wallHit(ID2)
#■画面コピー用マウスイベント  
def leftMouseDown(event):
    global ID1, ID2
    ID2=1-ID1; simulate(ID1, ID2);ID1=ID2
#■メイン処理
root=Tk();root.title("炎のシミュレーション")
canvas=Canvas(root, width=200,height=200,bg='black');canvas.pack()
root.update()
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:
        canvas.delete("DT")
        canvas.create_text(100,20,text="準備中です(%2d %%)" % int(i/N*100),
                           font=('Times',10),fill='yellow', tag="DT")
        canvas.create_text(100,40,text="■"*NK,
                           font=('Times',10),fill='yellow', tag="DT")
        NK+=1
        if NK>20: NK=1
        root.update()
    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)
canvas.delete("DT")
for i in range(200):
    for j in range(200):
        bitmap[i][j]= canvas.create_rectangle(i,j,i,j,fill="#000000",
                                              outline="",tag="DT")
simulate(ID1,ID2); ID1=ID2
#■左ボタンクリックしながら画面コピーするとき次の2行のコメントを外す
#canvas.bind("<Button-1>",leftMouseDown)  #画面コピー用マウス
#root.mainloop()
#■シミュレーション開始
while True:
    ID2=1-ID1
    simulate(ID1,ID2)
    time.sleep(0.001)
    ID1=ID2