# Newman-Ziff percolation code

from random import *
from math import *
import Numeric 
from visual import *
from visual.graph import *

L=32
N=(L*L)
EMPTY=(-N-1)
p=0.5
kcore=2
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))
   
occupation=Numeric.zeros((N))
ptr=Numeric.zeros((N))
for i in range (N):
    if random()<p:
        occupation[i]=1
    else:
        occupation[i]=0

change=1
while change==1:
    update=0
    for i in range (N): 
        xpos=i%L
        ypos=(i/L)
        if occupation[i]==1:
            sum=0
            if occupation[i-1+L*(xpos==0)]==1:
                sum=sum+1
            if occupation[i+1-L*(xpos==(L-1))]==1:
                sum=sum+1
            if occupation[i-L+N*(ypos==0)]==1:
                sum=sum+1
            if occupation[i+L-N*(ypos==(L-1))]==1:
                sum=sum+1
            if sum<kcore:
                occupation[i]=0
                update=update+1
    if update==0:
        change=0

for i in range (N):
    if (occupation[i]==1):
        ptr[i] = -1
    else:
        ptr[i]=EMPTY


        
for k in range (2):    
    for i in range (N):
        if ptr[i]!=EMPTY: 
            s1=i
            r1=findroot(s1);
            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

scene.background=color.white                            

max=0
for i in range(N):
    if EMPTY<ptr[i]<0:
        if i>max:
            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<ptr[i]<0:
        box(pos=(box_size*xpos-L/2,box_size*ypos-L/2,0),color=((1.0*i)/(max+1),(1.0*i)/(max+1),(1.0*i)/(max+1)))