ブラウン運動のシミュレーション可視化

以前ここで試した単純な粒子の中に大きいサイズの粒子を放り込む。周囲から無数の粒子がぶつかりランダムに動くブラウン運動のような動きをみることができる。

■結果

 
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