#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 13 13:04:10 2017

@author: laro
"""
from scipy import * 
from random import *
from pylab import *
import time
from numpy import cumsum



#este codigo es para ver la convergencia de las cadenas de markov

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)
                else:
                    g=1
                		  
    return DEc, DM, S


# Tamaño en x,
Lx=16
# Tamaño en y,
Ly=16
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 = 10000
#npasos=100


# 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
            
# =============================================================================
# #pretermalizo
# npre=100
# for n in range(1,npre):
#     DE,DM,Sij=ising2Dpaso(Sij,beta)
# 
# =============================================================================
            
            
energia=zeros(npasos)
magnet=zeros(npasos)
      
      
energia[0] = En(Sij)/(Lx*Ly)
magnet[0] = sum(Sij)/(Lx*Ly)
# =============================================================================
# ion()
# gcf()
# =============================================================================
count=0

for n in range(1,npasos):
    count=count+1
    DE,dM,Sij = ising2Dpaso(Sij,beta)
    energia[n] = energia[n-1] + DE/(Lx*Ly)
    magnet[n] = magnet[n-1] + dM/(Lx*Ly)
    if count%1000 ==0 :
        print(count)
        figure(int(count/1000))
        u=linspace(1,count,count)
        v=array(cumsum(magnet[:count]))/array(u)
        plot(u,v)
        
#voy ploteando la magnetizacion media acumulada, es decir, v es un vector cuyo elemento iesimo es
# el promedio de las magnetizaciones de los primeros i elementos de la cadena de markov