# Newman-Ziff percolation code from random import * from math import * import Numeric from visual import * from visual.graph import * L=8 N=(L*L) EMPTY=(-N-1) p=0.5 seed(22222) scene=display(x=0,y=400,width=400,height=400,fov=1e-8,background=color.white) pic1=gdisplay(x=0,y=0,width=400,height=400,title="Fraction of sites in largest cluster vs p", xtitle="p",ytitle="Fraction of sites in largest cluster",xmax=1,xmin=0,ymin=0,ymax=1,foreground=color.black,background=color.white) plot1=gcurve(gdisplay=pic1,color=color.blue) def findroot(i): # print "%4d" % (ptr[i]) if ptr[i]<0: return i else: ptr[i] = findroot(ptr[i]) return ptr[i] nn = Numeric.zeros((N,4)) for i in range (N): xpos=i%L ypos=int(i/L) nn[i,0]=i-1+L*(xpos==0) nn[i,1]=i+1-L*(xpos==(L-1)) nn[i,2]=i-L+N*(ypos==0) nn[i,3]=i+L-N*(ypos==(L-1)) order=Numeric.zeros((N)) for i in range (N): order[i] = i for i in range (N): j=int(i+(N-i)*random()) temp = order[i] order[i] = order[j] order[j] = temp # print "%4d\t%4d" % (i,order[i]) big =0 ptr=Numeric.zeros((N)) for i in range (N): ptr[i] = EMPTY for i in range (N): r1 = order[i] s1 = order[i] ptr[s1] = -1 for j in range (4): s2 = nn[s1][j] if ptr[s2]!=EMPTY: r2 = findroot(s2) if r2!=r1: if ptr[r1]>ptr[r2]: ptr[r2] += ptr[r1] ptr[r1] = r2 r1 = r2 else: ptr[r1] += ptr[r2] ptr[r2] = r1 # print "%4d\t%4d" % (ptr[r1],order[i])) if -ptr[r1]>big: big = -ptr[r1] print "%4d\t%4d" % (i+1,big) plot1.plot(pos=( float(i+1)/float(N) , float(big)/float(N) )) for i in range (N): if iptr[r2]): ptr[r2] += ptr[r1] ptr[r1] = r2 r1 = r2 else: ptr[r1] += ptr[r2] ptr[r2] = r1 scene.background=color.white max=0 for i in range(N): if EMPTYmax: max=i box_size=1.1 for i in range (N): # print "%4d\t%4d" % (i,ptr[i]) xpos=i%L ypos=(i/L) # print "%4d\t%4d" % (xpos,ypos) if ptr[i]>=0: box(pos=(box_size*xpos-L/2,box_size*ypos-L/2,0),color=((1.0*ptr[i])/(max+1),(1.0*ptr[i])/(max+1),(1.0*ptr[i])/(max+1))) if EMPTY