import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt

N = 100    # taille de la grille (on suppose N pair pour simplifier les interactions au bord)
betaCritique = np.log(1+np.sqrt(2))/2

def Proba(i,j,M,beta): # proba que le site (i,j) = 1 sachant les voisins 
    D = M[i-1,j] + M[(i+1)%N,j] + M[i,(j+1)%N] + M[i,j-1]
    return 1/(1+np.exp(-2*beta*D))

def EtapeBlanche(M,beta):
    for i in range(N):
        for j in range(N):
            if (i-j)%2 == 0:
                if npr.random() < Proba(i,j,M,beta):
                    M[i,j] = 1
                else:
                    M[i,j] = -1
    return(M)

            
def EtapeNoire(M,beta):
    for i in range(N):
        for j in range(N):
            if (i-j)%2 == 1:
                if npr.random() < Proba(i,j,M,beta):
                    M[i,j] = 1
                else:
                    M[i,j] = -1
    return(M)


def Gibbs(Nstep,M0,beta):
    M = M0
    for i in range(Nstep):
        if i%2 == 0: 
            M = EtapeBlanche(M,beta)
        else:
            M = EtapeNoire(M,beta)
    return(M)

# Conditions initiales 

def M0rand(N): return(2*np.random.randint(2,size=(N,N))-1) # aléatoire
def M0align(N): return(np.ones((N,N))) # aligné

# Illustration 

Nstep = 500
M0 = M0align(N)
beta = betaCritique-0.1

IsingGibbs = Gibbs(Nstep,M0,beta)

plt.imshow(IsingGibbs)
plt.show()