ブラウン運動のシミュレーション可視化
以前ここで試した単純な粒子の中に大きいサイズの粒子を放り込む。周囲から無数の粒子がぶつかりランダムに動くブラウン運動のような動きをみることができる。
■結果
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 |
import random import math from OpenGL.GL import * from OpenGL.GLU import * from OpenGL.GLUT import * from PIL import Image from PIL import ImageOps a = 2 Dt = 0.01 CUT = 5 #Define xy index PX = 0 PY = 1 VX = 2 VY = 3 FX = 4 FY = 5 pot = [0,0] perticle_num = 1000 xy = [[0 for i in range(6)] for j in range(perticle_num)] BP = [0 for i in range(6)] Radi = 0.1 step = 1 def capture(): global step pad_step = '{0:04d}'.format(step) print pad_step width = glutGet(GLUT_WINDOW_WIDTH) height = glutGet(GLUT_WINDOW_HEIGHT) glReadBuffer(GL_FRONT) glPixelStorei(GL_UNPACK_ALIGNMENT, 1) data = glReadPixels(0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE) image = Image.fromstring("RGBA", (width, height), data) image = ImageOps.flip(image) image.save("test"+pad_step+".png") def find_pair(): global Radi particle_num = len(xy) pnum = 0 w_ratio = 0.1 while pnum < particle_num: dx = xy[pnum][PX] - BP[PX] dy = xy[pnum][PY] - BP[PY] if dx > 0.5*a: dx = dx + a if dy > 0.5*a: dy = dy + a rr = math.sqrt(dx*dx + dy*dy) if rr < Radi: v2 = math.sqrt(xy[pnum][VX]**2 + xy[pnum][VY]**2) cos = (xy[pnum][VX]*dx + xy[pnum][VY]*dy)/(v2*rr) sin = math.sqrt(1 - cos**2) tmpx= xy[pnum][VX]*(cos**2 - sin**2) - xy[pnum][VY]*(2*cos*sin) tmpy= xy[pnum][VX]*(2*cos*sin) + xy[pnum][VY]*(cos**2 - sin**2) xy[pnum][VX] = abs(tmpx) * dx / (abs(dx)) xy[pnum][VY] = abs(tmpy) * dy / (abs(dy)) BP[VX] = BP[VX] - w_ratio * abs(tmpx) * dx / (abs(dx)) BP[VY] = BP[VY] - w_ratio * abs(tmpy) * dy / (abs(dy)) pnum += 1 def integ(): find_pair() particle_num = len(xy) BP[PX] = BP[PX] + Dt*BP[VX] BP[PY] = BP[PY] + Dt*BP[VY] #boundary condition if BP[PX] > a: BP[PX] -= a if BP[PY] > a: BP[PY] -= a if BP[PX] < 0: BP[PX] += a if BP[PY] < 0: BP[PY] += a pnum = 0 while pnum < particle_num: #update position xy[pnum][PX] = xy[pnum][PX] + Dt*xy[pnum][VX] xy[pnum][PY] = xy[pnum][PY] + Dt*xy[pnum][VY] if xy[pnum][PX] > a: xy[pnum][PX] -= a if xy[pnum][PY] > a: xy[pnum][PY] -= a if xy[pnum][PX] < 0: xy[pnum][PX] += a if xy[pnum][PY] < 0: xy[pnum][PY] += a pnum += 1 def calc_ena(): ena = 0.0 particle_num = len(xy) for i in range(0,particle): vx2 = xy[i][VX]*xy[i][VX] vy2 = xy[i][VY]*xy[i][VY] ena = ena + vx2 + vy2 ena = 0.5*ena return ena def update(): global Radi PI = 3.14159263 global step step = step + 1 # If you want to get image... #capture() integ() #print step nnum = 0 while nnum < len(xy): glPointSize(3.0) glColor3f(0.3, 0.3, 1.0) glBegin(GL_POINTS) glVertex2d(xy[nnum][PX], xy[nnum][PY]) glEnd() nnum = nnum + 1 N_t = 30 glColor3f(1.0, 0.0, 0.0) glBegin(GL_POLYGON) for j in range(0,N_t): glVertex2d(Radi*math.cos(2.0*PI*j/N_t)+BP[PX], Radi*math.sin(2.0*PI*j/N_t)+BP[PY]) glEnd() def init_set(): pnum = 0 xv_sum = 0 yv_sum = 0 particle_num = len(xy) while pnum < particle_num: xy[pnum][PX] = 0 # 0 x-position xy[pnum][PY] = 0 # 1 y-position xy[pnum][VX] = random.uniform(-1,1) # 2 x-velocity xy[pnum][VY] = random.uniform(-1,1) # 3 y-velocity xv_sum += xy[pnum][VX] yv_sum += xy[pnum][VY] pnum += 1 xv_sum = xv_sum / particle_num yv_sum = yv_sum / particle_num pnum = 0 while pnum < particle_num: xy[pnum][VX] = xy[pnum][VX] - xv_sum xy[pnum][VY] = xy[pnum][VY] - yv_sum pnum += 1 BP[PX] = a/2 BP[PY] = a/2 BP[VX] = random.uniform(-0.5,0.5) BP[VY] = random.uniform(-0.5,0.5) def draw(): glClearColor(0.0, 0.0, 0.0, 0.0) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) glLoadIdentity() glColor3f(1.0,1.0,0.0) gluLookAt(a/2, a/2, 4.0 ,a/2 ,a/2 ,0 , 0, 1.0, 0.0) glBegin(GL_LINE_LOOP) glVertex2d(0,0) glVertex2d(0,a) glVertex2d(a,a) glVertex2d(a,0) glEnd() update() glutSwapBuffers() def resize(w, h): glViewport(0, 0, w, h) glLoadIdentity() glOrtho(0, a, 0, a, -a, a) def init(): glClearColor(0.7, 0.7, 0.7, 0.7) def idle(): glutPostRedisplay() def reshape(w, h): glViewport(0, 0,w,h) glMatrixMode(GL_PROJECTION) glLoadIdentity() gluPerspective(30.0, w/h, 1.0, 100.0) glMatrixMode (GL_MODELVIEW) if __name__ == "__main__": init_set() glutInitWindowPosition(50, 50); glutInitWindowSize(340, 340); glutInit(sys.argv) glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE ) glutCreateWindow("pyOpenGL TEST") glutDisplayFunc(draw) glutReshapeFunc(reshape) init() glutIdleFunc(idle) glutMainLoop() |
brownian motion simulation visualization by python and openGL
pyOpenGLを使って分子動力学を可視化
以前ここで試したmatplotlibを使った簡単な分子動力学計算の可視化処理ですが、描画が遅いのでOpenGLを使って書き直してみた。python用にはPyOpenGLのパッケージが用意されている。結果だがとても速い！
■結果をgifアニメに落としたものが↓。
gifアニメだと実際の描画より遅いので感動が伝わりにくい。。。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 |
import random from OpenGL.GL import * from OpenGL.GLU import * from OpenGL.GLUT import * from PIL import Image from PIL import ImageOps a = 10 Dt = 0.1 #Define xy index PX = 0 PY = 1 VX = 2 VY = 3 FX = 4 FY = 5 perticle_num = 1000 xy = [[0 for i in range(6)] for j in range(perticle_num)] step = 1 def capture(): global step step = step + 1 pad_step = '{0:04d}'.format(step) print pad_step width = glutGet(GLUT_WINDOW_WIDTH) height = glutGet(GLUT_WINDOW_HEIGHT) glReadBuffer(GL_FRONT) glPixelStorei(GL_UNPACK_ALIGNMENT, 1) data = glReadPixels(0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE) image = Image.fromstring("RGBA", (width, height), data) image = ImageOps.flip(image) image.save("test"+pad_step+".png") def integ(): #verlet integration particle_num = len(xy) pnum = 0 while pnum < particle_num: #update velocity xy[pnum][VX] = xy[pnum][VX] + 0.5*Dt*xy[pnum][FX] xy[pnum][VY] = xy[pnum][VY] + 0.5*Dt*xy[pnum][FY] #update position xy[pnum][PX] = xy[pnum][PX] + Dt*xy[pnum][VX] xy[pnum][PY] = xy[pnum][PY] + Dt*xy[pnum][VY] #boundary condition if xy[pnum][PX] > a: xy[pnum][PX] -= a if xy[pnum][PY] > a: xy[pnum][PY] -= a if xy[pnum][PX] < 0: xy[pnum][PX] += a if xy[pnum][PY] < 0: xy[pnum][PY] += a pnum += 1 def update(): # If you want to get image... # capture() integ() nnum = 0 while nnum < len(xy): glPointSize(3.0) glColor3f(0.3, 0.3, 1.0) glBegin(GL_POINTS) glVertex2d(xy[nnum][PX], xy[nnum][PY]) glEnd() nnum = nnum + 1 def init_set(): pnum = 0 xv_sum = 0 yv_sum = 0 particle_num = len(xy) while pnum < particle_num: xy[pnum][PX] = 0 # 0 x-position xy[pnum][PY] = 0 # 1 y-position xy[pnum][VX] = random.uniform(-1,1) # 2 x-velocity xy[pnum][VY] = random.uniform(-1,1) # 3 y-velocity xv_sum += xy[pnum][VX] yv_sum += xy[pnum][VY] xy[pnum][FX] = 0 # 0 x-force xy[pnum][FY] = 0 # 1 y-force pnum += 1 xv_sum = xv_sum / particle_num yv_sum = yv_sum / particle_num pnum = 0 while pnum < particle_num: xy[pnum][VX] = xy[pnum][VX] - xv_sum xy[pnum][VY] = xy[pnum][VY] - yv_sum pnum += 1 #exit() def draw(): glClearColor(0.0, 0.0, 0.0, 0.0) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) glLoadIdentity() glColor3f(1.0,1.0,0.0) gluLookAt(a/2, a/2, 20.0 ,a/2 ,a/2 ,0 , 0, 1.0, 0.0) glBegin(GL_LINE_LOOP) glVertex2d(0,0) glVertex2d(0,a) glVertex2d(a,a) glVertex2d(a,0) glEnd() update() glutSwapBuffers() def init(): glClearColor(0.7, 0.7, 0.7, 0.7) def idle(): glutPostRedisplay() def reshape(w, h): glViewport(0, 0,w,h) glMatrixMode(GL_PROJECTION) glLoadIdentity() gluPerspective(30.0, w/h, 1.0, 100.0) glMatrixMode (GL_MODELVIEW) if __name__ == "__main__": init_set() glutInitWindowPosition(50, 50); glutInitWindowSize(340, 340); glutInit(sys.argv) glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE ) glutCreateWindow("pyOpenGL TEST") glutDisplayFunc(draw) glutReshapeFunc(reshape) init() glutIdleFunc(idle) glutMainLoop() |
matplotlibでランダムウォーク
ランダムウォークをmatplotlibで表示してみた。
各粒子の動作は乱数で上下左右に動くように設定。
マシンが非力なせいか、粒子数を増やすとgif作成にやたら時間がかかる。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 |
#-*- coding:utf-8 -*- import matplotlib.animation as anm import matplotlib.pyplot as plt import numpy as np import random def pstep(xy): nnum = 0 pnum = 0 while nnum < len(xy): while pnum < 2: rr = random.random() if rr > 0.5: xy[nnum][pnum] = xy[nnum][pnum] + 0.05 else: xy[nnum][pnum] = xy[nnum][pnum] - 0.05 if xy[nnum][pnum] < 0: xy[nnum][pnum] += 10 if xy[nnum][pnum] > 10: xy[nnum][pnum] -= 10 pnum = pnum + 1 pnum = 0 nnum = nnum + 1 def update(i, fig_title, A, xy, total_step): if i != 0: plt.cla() ax.set_xlim(0,10) ax.set_ylim(0,10) pstep(xy) nnum = 0 while nnum < len(xy): plt.plot(xy[nnum][0], xy[nnum][1], "o") nnum = nnum + 1 plt.title(fig_title + ' step = ' + str(i) + '/' + str(total_step)) if __name__ == "__main__": fig = plt.figure(figsize = (10, 6)) fig.set_size_inches(10,10) ax = fig.add_subplot(111) perticle_num = 50 xy = [[0 for i in range(2)] for j in range(perticle_num)] nnum = 0 pnum = 0 while nnum < len(xy): while pnum < 2: xy[nnum][pnum] = random.random() * 10 pnum = pnum + 1 pnum = 0 nnum = nnum + 1 total_step =200 ani = anm.FuncAnimation(fig, \ update, fargs = ('Random Work', 2.0, xy, total_step), \ interval = 100, frames = total_step) ##plt.show() ani.save("Sample_randomwork.gif", writer = 'imagemagick', dpi=50) |
実行した結果が下記。
random work sample by python and matplotlib.