#!/usr/bin/python -O 
# -*- encoding: utf-8 -*-
from scipy import * 
from random import *
from pylab import *
import time
# script general para hacer una corrida  a un set de parámetros,
# beta, tamaño de la red,

#beta = 1/T

def En(S):
    energia = 0
    for i in range(0,Lx):
        iright = irightl[i];#mod(i+1,Lx) 
        ileft = ileftl[i];#mod(i+Lx-1,Lx) 
        for j in range (0, Ly):
            jup = jupl[j]; #mod(j+1,Ly) 
            jdown = jdownl[j]; #mod(j+Ly-1,Ly) 
            energia += - S[i,j]*(S[i,jup] + S[i,jdown] + S[iright,j] + S[ileft,j]);

    return energia/2.;

def ising2Dpaso(S,beta):

    DEc = 0
    DM = 0

    for i in range (0,Lx):
        iright = irightl[i];#mod(i+1,Lx)
        ileft = ileftl[i];#mod(i+Lx-1,Lx)
        for j in range (0, Ly):
            jup = jupl[j];#mod(j+1,Ly)
            jdown = jdownl[j];#mod(j+Ly-1,Ly) 
            spin_old = S[i,j]
            spin_new = -S[i,j]

            DE = 2.0*spin_old*(S[iright,j] + S[ileft,j] + S[i,jup] + S[i,jdown] )

	    if( DE < 0) : 
		S[i,j] = spin_new
		DEc += DE
		DM += (spin_new - spin_old)
	    else:
		p = random()
		expbetaE = exp(-beta*DE)
		if (p < expbetaE):
		    S[i,j] = spin_new
		    DEc += DE
		    DM += (spin_new - spin_old)
		
    return DEc, DM


# Tamaño en x,
Lx=8
# Tamaño en y,
Ly=8
irightl=roll(arange(0,Lx),1);
ileftl=roll(arange(0,Lx),-1);
jupl=roll(arange(0,Ly),1);
jdownl=roll(arange(0,Ly),-1);

T=2;
beta=1.0/T;	
npasos = int(1e5)


# propongo un estado inicial al azar
  
# Sij es una matriz de 1 y -1 indicando las dos proyecciones de
# espin
Sij = zeros((Lx,Ly))
for i in range(0,Lx):
    for j in range (0,Ly):
        if(random() > 0.5):
            Sij[i,j] = 1 
        else:
            Sij[i,j] = -1
            
            
energia=zeros(npasos);
magnet=zeros(npasos);
      
      
energia[0] = En(Sij)/(Lx*Ly);
magnet[0] = sum(Sij)/(Lx*Ly);
ion()
gcf()
for n in range(1,npasos):
    DE,dM = ising2Dpaso(Sij,beta);
    energia[n] = energia[n-1] + DE/(Lx*Ly);
    magnet[n] = magnet[n-1] + dM/(Lx*Ly);

    if(mod(n,100) == 0):
            print energia[n], magnet[n];
            pcolor(Sij);
            pause(.5);
            draw();
            time.sleep(0.33);
   

print "# ", T,mean(energia),mean(abs(magnet)),std(energia),std(abs(magnet))