from scipy import *
from pylab import *
import time
from scipy import weave
from scipy.weave import converters

n=64

dx = zeros((4),dtype=int)
dy = zeros((4),dtype=int)


dx[0] = 1
dy[0] = 0

dx[1] = -1
dy[1] = 0

dx[2] = 0
dy[2] = 1

dx[3] = 0
dy[3] = -1 


beta = 1
s = ones((n,n),dtype=int)
r = random((n,n))
for x in range(n):
    for y in range(n):
        if r[x,y] < 0.5:
            s[x,y] = -1

def iterate(s,beta):
   r = random((n,n))
   dir = randint(0,4, (n,n))
   code = """
    #line 40 "conserved.py"
    int x,y,x1,y1;
    int mod = n-1;
    for(x=0; x < n; x++) {
        for(y=0; y < n; y++) {
            x1 = (x + dx(dir(x,y)))&mod;
            y1 =  (y+ dy(dir(x,y)))&mod;

            if (s(x,y) != s(x1,y1)) {
                float field1 = s(((x+1)&mod),y) + s(((x-1)&mod),y) + s(x,((y+1)&mod)) + s(x,((y-1)&mod));
                float field2= s(((x1+1)&mod),y1) + s(((x1-1)&mod),y1) + s(x1,((y1+1)&mod)) + s(x1,((y1-1)&mod));
                float delta_e = s(x,y)*(field1-field2) + 4;
                if (exp(-beta*delta_e) > r(x,y)) {
                    s(x,y) = -s(x,y);
                    s(x1,y1) = -s(x1,y1);
                }
            }
        }
    }
           """
   weave.inline(code, ['s', 'r', 'dir', 'beta', 'dx', 'dy', 'n'], type_converters=converters.blitz, compiler = 'gcc')
   return


ion()
       
for k in range(500000):
       iterate(s,beta)
       if (k%2000 == 0):
	       #imshow(abs(fft2(s)))
	       imshow(s)