from scipy import *
import random as r
from pylab import *
import matplotlib.animation as animation
from numba import jit
@jit(nopython=True)
def GlueSingleWalker(lattice):
""" Generates a new random walker on a boundary, and walks it around
until it comes close to the grown cluster. It then glues to the cluster.
"""
N = len(lattice)
# Random walker is released from the boundary"
while True:
istart = int(r.random()*4*N)
icase = int(istart/N) # icase can be either 0,1,2,3
if int(icase/2): # can be 0 or 1
x0, y0 = istart%N, (0,N-1)[icase%2] # either top or button of the figure, and x is random number.
else:
x0, y0 = (0,N-1)[icase%2], istart%N # either left or right, and y is random number
if lattice[x0,y0]==0: break # this is empty site, hence walker can start
while True: # Random walker is walking until it glues with the cloud
x,y = x0,y0 # Random walker before being moved
ipos = int(r.random()*4) # one of the four moves is choosen randomly
x_or_y = int(ipos/2) # x_or_y is either 0 or 1
dr = 2*int(ipos%2)-1 # dr can be -1 or +1
if x_or_y:
x = (x+dr)%N # moves left one step. notice periodic boundary conditions
else:
y = (y+dr)%N # moves up on step
# Checking if particle can be glued!
# Here we do not take periodic boundary conditions since they put too many particles on the boundary
r_neighbor = x0 and lattice[x-1,y] # neighbor left
d_neighbor = y>0 and lattice[x,y-1] # neighbor below
if r_neighbor or u_neighbor or l_neighbor or d_neighbor: # Particle must be glued
lattice[x,y]=1
#print 'Glued', x, y
break # exit the infinite loop once the walker is glued.
x0,y0=x,y
def Mass(lattice):
dr=20
Rmax = sqrt(2.)/2*len(lattice) # maximum radius the particle could be
Nr = int(Rmax/dr)+1 # number of rings
mass = zeros(Nr)
N=len(lattice)
for i in range(N):
for j in range(N):
if lattice[i][j]==1: # this site is occupied
rx = i-N/2 # the distance of the point from the center in x-direction
ry = j-N/2 # the distance of the point from the center in y-direction
r = sqrt(rx**2+ry**2) # distance from the origin
mass[int(r/dr)]+=1 # the ring at r/dr containe one more particle
for i in range(1,Nr):
mass[i] += mass[i-1] # cumulative sum -> mass inside the circle of radius dr*(i+1)
fi = open('mass.dat', 'w')
for i in range(Nr):
print(dr*(i+1), mass[i], file=fi)
fi.close()
N=400
lattice = zeros((N,N),dtype=int)
lattice[int(N/2),int(N/2)]=1 # We put first point in the middle of the system
fig = figure(figsize=(8, 8)) # make figure
# make axesimage object
# the vmin and vmax here are very important to get the color map correct
im = imshow(lattice, cmap=cm.Greys, vmin=0, vmax=1, interpolation='nearest')
def updatefig(i):
for j in range(50):
GlueSingleWalker(lattice)
# set the data in the axesimage object
im.set_array(lattice)
if i%50==1:
print('Computing mass')
Mass(lattice)
# return the artists set
return im,
ani = animation.FuncAnimation(fig, updatefig, frames=100, blit=True)
show()