
import numpy
import numpy.random
import random
import math

class Ising1D:
    def __init__(self,N):
        self.N = N
        self.spin = numpy.ones(N)
        self.E = 0
        self.M = N
                    
    def constantes(self,B,T):
        beta = 1.0/T
        self.dE = [2*B+4,2*B-4,2*B,-2*B+4,-2*B-4,-2*B]
        self.A = []
        for k in range(6):
            self.A.append(math.exp(-beta*self.dE[k]))
                    
    def voisins(self,i):
        if i==0:
            s1 = self.spin[self.N-1]
        else:
            s1 = self.spin[i-1]
        if i==self.N-1:
            s3 = self.spin[0]
        else:
            s3 = self.spin[i+1]
        return (s1,s3)
                    
    def deltaE(self,i):
        s2 = self.spin[i]
        (s1,s3) = self.voisins(i)
        if s2==1:
            p = 0
        else:
            p = 3
        if (s1,s2,s3)==(1,1,1) or (s1,s2,s3)==(-1,-1,-1):
            pass
        elif (s1,s2,s3)==(1,-1,1) or (s1,s2,s3)==(-1,1,-1):
            p += 1
        else:
            p += 2
        return (self.dE[p],self.A[p])
                    
    def _metropolis(self):
        i = random.randint(0,self.N)
        (dE,A)= self.deltaE(i)
        if dE<=0:
            self.spin[i] = -self.spin[i]
            self.E += dE
            self.M += 2*self.spin[i]
        else:
            x = random.random()
            if x<A:
                self.spin[i] = -self.spin[i]
                self.E += dE
                self.M += 2*self.spin[i]
                
            
                    
    def _boucle(self,n):
        m = numpy.zeros(n)
        for k in range(n):
            m[k] = self.M
            for i in range(self.N):
                self.metropolis()
        return (m,numpy.mean(m),numpy.std(m))
                    
    def metropolis(self,i,x):
        (dE,A)= self.deltaE(i)
        if dE<=0:
            self.spin[i] = -self.spin[i]
            self.E += dE
            self.M += 2*self.spin[i]
        else:
            if x<A:
                self.spin[i] = -self.spin[i]
                self.E += dE
                self.M += 2*self.spin[i]
                
    def boucle(self,n):
        m = numpy.zeros(n)
        tab_i = numpy.random.randint(0,self.N,size=n*self.N)
        tab_x = numpy.random.random_sample(size=n*self.N)
        j = 0
        for k in range(n):
            m[k] = self.M
            for i in range(self.N):
                self.metropolis(tab_i[j],tab_x[j])
                j += 1
        return (m,numpy.mean(m),numpy.std(m))
                      