### DA COMPLETARE

#!/usr/bin/ruby
class Rng

    # private:

   # cambiando i due numeri seguenti si può avere un diverso generatore
   # Fibonacci lagged
   Gensize = 250
   Genlag = 103
   A_LCM = 16807
   M_LCM = 2147483647
   Q_LCM = 127773
   R_LCM = 2836

    def initialize(seed = 1)

       #qui sono memorizzati gli ultimi 'Gensize' numeri random interi
       @buffer =  []
       # il seme e' il numero random iniziale
       @k = seed

       # parto con 'Gensize' numeri random generati da LCM
       for i in 0...Gensize
          # moltiplicazione modulo m che evita l'overflow
          @k=A_LCM*(@k % Q_LCM) - R_LCM*(@k/Q_LCM)
          # se il risultato e' negativo aggiungo m
          if(@k < 0)
             @k += M_LCM
          end

          @buffer[i] = @k

       end

       # inizializzazione di r250 presa dalle GSL

       msb = 0x80000000;
       mask = 0xffffffff;

       for i in 0...32
          k = 7 * i + 3;       # Select a word to operate on
          @buffer[k] &= mask;   # Turn off bits left of the diagonal
          @buffer[k] |= msb;    # Turn on the diagonal bit
          mask >>= 1;
          msb >>= 1;
       end

       # fine parte GSL

       @pointer = Gensize - 1 # punta all'ultimo posto riempito

    end

    def irand()

       # lfg prng, run after lfg_init

       # il nuovo numero e' generato dal vecchio sommando
       # due dei 'gensize' numeri con xor (cioe' bit a bit modulo 2)
       # pointer punta all'ultimo elemento trovato
       # rispetto al nuovo elemento n punta a n-'gensize'
       # il secondo indice e' n-103, ovvero (n+147)%'gensize'
       # ma rispetto a pointer va avanzato di uno
       newrand = @buffer[@pointer] ^ @buffer[(@pointer + Gensize - Genlag + 1) % Gensize]
       @pointer = (@pointer + 1 ) % Gensize
       @buffer[@pointer] = newrand
       # il nuovo numero e' inserito nello stack circolare

       return newrand

    end

    def uniform

       # genero un numero in doppia precisione compreso tra 0 e 1
       # a partire da un intero
       return irand.to_f/2147483647.0

    end

    def boxmul

       #
       # algoritmo di Box Mueller
       # genero un numero random con distribuzione gaussiana
       # a partire da due distribuiti uniformemente
       #

       x1 = uniform();
       x2 = uniform();

       return Math.sqrt( - 2*Math.log(x1))*cos(2*Math::PI*x2)

    end

    def esponenziale

       return -Math.log((uniform()).abs)

    end

end # fine classe Rng

class Ising

   # private:

   # dimensioni del reticolo
   NSIZE=50

   Nthermal=50
   Nsweep=10

   def initialize

      myrand = Rng.new(12783)
      # mi serve un generatore di numeri random
      # per inizializzare gli spin
      # alloco memoria per gli spin
      @spin = []
      @cluster = []
      # Scelgo la direzione iniziale degli spin a caso
      # con uguale probabilita' tra 1 e -1
      for j in 0...NSIZE*NSIZE
         if(myrand.uniform()>0.5)
            @spin[j] = 1
         else
            @spin[j] = -1
         end
      end
   end

   # restituisce la dimensione lineare del reticolo
   def getSize
      return NSIZE
   end

   def Energy()

   #
   # contributo di ogni legame all'energia. Ogni legame
   # viene contato una sola volta anziche' due,
   # per cui il risultato finale va moltiplicato per due
   isum=0
   # calcolo prima i legami che stanno sulla stessa riga
   # ciclo sulle righe
   for i in 0...NSIZE
      # ciclo sulle colonne
      for j in 0...NSIZE-1
         isum += @spin[i*NSIZE+j]*@spin[i*NSIZE+j+1]
      end
      isum += @spin[i*NSIZE+NSIZE-1]*@spin[i*NSIZE+0]
   end
   # calcolo ora i legami sulla stessa colonna
   # ciclo sulle colonne
   for j in 0...NSIZE
      # cicolo sulle righe
      for i in 0...NSIZE-1
         isum += @spin[i*NSIZE+j]*@spin[(i+1)*NSIZE+j]
      end
      isum += @spin[(NSIZE-1)*NSIZE+j]*@spin[0*NSIZE+j]
   end
   isum=2*isum

   # Energia senza campo magnetico
   return -isum.to_f / (NSIZE*NSIZE).to_f

   end

   def Magnetization()
   #
   # calcola la magnetizzazione media, data dalla somma degli spin
   # divisa per il loro numero
   #
   isum=0

   #
   # sommo tutti gli spin
   #
   for i in 0...NSIZE*NSIZE
      isum += @spin[i]
   end

   #
   # prendo il valore assoluto della magnetizzazione
   # dato che il sistema e' invariamnte per una inversione di tutti gli spin
   #
   return isum.to_f.abs/(NSIZE*NSIZE).to_f

   end

   def Metropolis(beta, myrand)

   # esegue n "passate" con l'algoritmo di Metropolis su tutto il reticolo
   f4=Math.exp(-8*beta)    # k = 4
   f2=Math.exp(-4*beta)    # k = 2

   # effettuo Nsweep passate su tutti i NSIZE x NSIZE spin del reticolo
   for m in 0...Nsweep
      # loop sulle righe
      for i in 0...NSIZE
         # indici di riga vicini ad 'i'
         ip1=(i+1)%NSIZE
         im1=i-1
         if i==0
            im1=NSIZE-1  # condizioni al contorno periodiche
         end
         # loop sulle colonne
         for j in 0...NSIZE
            # indici di colonna vicini a 'j'
            jp1=(j+1)%NSIZE
            jm1=j-1
            jm1=NSIZE-1 if(j==0)
            # ora calcolo la variazione di energia
            # nel flippare lo spin in (i,j)
            k=@spin[ip1*NSIZE+j]+@spin[im1*NSIZE+j]+\
              @spin[i*NSIZE+jp1]+@spin[i*NSIZE+jm1]
            deltaE= 2*@spin[i*NSIZE+j]*k
            # se l'energia diminuisce accetto la nuova situazione
            if deltaE <= 0
               @spin[i*NSIZE+j] = -@spin[i*NSIZE+j]
            # altrimenti la accetto con probabilita' exp(-beta*Delta_E)
            # genero un numero a caso tra 0 e 1 e accetto la variazione
            # se questo e' minore di Delta_E
            else
               case k.abs
                  when 4
                     @spin[i*NSIZE+j] = -@spin[i*NSIZE+j] if(f4 > myrand.uniform())
                     break
                  when 2
                     @spin[i*NSIZE+j] = -@spin[i*NSIZE+j] if(f2 > myrand.uniform())
                     break
               end
            end # fine if..else per determinare se passare alla nuova configurazione
         end # fine ciclo sulle colonne
      end # fine ciclo sulle righe
    end # fine ciclo sul numero di passaggi
    end # fine funzione

   def ExactEnergy(beta)

      q1 =2*Math.sinh(2*beta)/(Math.cosh(2*beta)**2)
      q2 = 2*(Math.tanh(2*beta)**2)-1.0
      h = Math::PI/2.0/1000.0

      # questo e' un integrale che compare nella
      # formula analitica per l'energia
      # lo calcolo con la regola del trapezio
      sum=0.0
      for j in 0..1000
         x=j*h
         sum += 1.0/Math.sqrt(1.0-(q1*Math.sin(x))**2)
      end
      sum=sum-0.5*(1.0+1.0/Math.sqrt(1.-q1*q1))
      sum=sum*h

      return -2.0/Math.tanh(2*beta)*(1+2*q2/Math::PI*sum)

   end

   def Wolff(beta, myrand)

      # probabilita' di aggiungere un sito al cluster
      pAdd = 1.0-Math.exp(-2*beta)

      for m in 0...Nsweep
         # nessuno spin appartiene al cluster
         for k in 0...NSIZE*NSIZE
            @cluster[k]=0
         end
         # prendo uno spin a caso
         i = (myrand.uniform()*NSIZE).floor
         j = (myrand.uniform()*NSIZE).floor

         # guardo quanto vale lo spin
         s=@spin[i*NSIZE+j]

         # segnalo questo sito come gia' appartenente al cluster
         @cluster[i*NSIZE+j]=1

         # aggiungi elementi al cluster in modo ricorsivo
         addcluster(i,j,pAdd,s,myrand)

         # cambia il segno di tutti gli spin nel cluster
         for k in 0...NSIZE
            for n in 0...NSIZE
               @spin[k*NSIZE+n]= -s if(@cluster[k*NSIZE+n]==1)
            end
         end
      end
   end

   def addcluster(i, j, pAdd, s, myrand)

      # funzione che serve a calcolare il cluster
      # viene chiamata in modo ricorsivo

      # aggiunge al cluster i primi vicini con probabilita' pAdd
      # se gia' non appartengono al cluster e se hanno lo stesso spin
      # determino gli indici di riga e colonna dei primi vicini
      ip1 = i+1
      ip1=0 if ip1 == NSIZE
      jp1 = j+1
      jp1=0 if jp1 == NSIZE
      im1 = i-1
      im1=NSIZE-1 if i == 0
      jm1 = j-1
      jm1=NSIZE-1 if j == 0

      #
      # chiamo addcluster() in modo ricorsivo per ciascuno dei primi vicini
      #
      if(@cluster[ip1*NSIZE+j]==0 && @spin[ip1*NSIZE+j]==s)
         # aggiungo lo spin con probabilità pAdd al cluster
         if(myrand.uniform()< pAdd)
            @cluster[ip1*NSIZE+j]=1
            # ripeto il procedimento per i primi vicini
            addcluster(ip1,j,pAdd,s, myrand)
         end
      end
      if(@cluster[im1*NSIZE+j]==0 && @spin[im1*NSIZE+j]==s)
         if(myrand.uniform()< pAdd)
            @cluster[im1*NSIZE+j]=1
            addcluster(im1,j,pAdd, s, myrand)
         end
      end
      if(@cluster[i*NSIZE+jp1]==0 && @spin[i*NSIZE+jp1]==s)
         if(myrand.uniform()< pAdd)
            @cluster[i*NSIZE+jp1]=1
            addcluster(i,jp1,pAdd, s, myrand)
         end
      end
      if(@cluster[i*NSIZE+jm1]==0 && @spin[i*NSIZE+jm1]==s)
         if(myrand.uniform()< pAdd)
            @cluster[i*NSIZE+jm1]=1
            addcluster(i,jm1,pAdd, s, myrand)
         end
      end
   end

   def StampaReticolo()
      print "Reticolo\n\n"
      for i in 0...NSIZE
         for j in 0...NSIZE
            if(@spin[i*NSIZE+j]==1)
               print "o"
            else
               print"x"
            end
         end
         print "\n"
      end
   end

end

lattice = Ising.new
myrand = Rng.new

nthermal=300
naverage=100
niterazioni=50
betaini=0.439
betafin=0.465
betastep=0.002

# lattice.StampaReticolo()

print "\nReticolo degli spin ",lattice.getSize(),"x",lattice.getSize(),"\n"
print "Passaggi di termalizzazione ",nthermal,"\n"
print "Media su ",naverage," stati","\n"
print "calcolata su stati dopo ",niterazioni," iterazioni","\n"

# termalizzazione preliminare
lattice.Wolff(betaini,myrand)
#lattice.StampaReticolo()
print "\n beta  Energia  E esatta  Magnetizzazione\n"

beta = betaini - betastep
while beta < betafin

   beta += betastep

   # termalizzazione
   lattice.Wolff(beta,myrand)

   # StampaReticolo()
   # azzero i valori medi di energia e magnetizzazione
   totalEnergy=0
   totalMagnetization=0.0
   for k in 1..naverage
      lattice.Wolff(beta,myrand)
      # Calcolo dei valori indipendenti dell'energia e della magnetizzazione
      # Quindi li sommo.
      totalEnergy += lattice.Energy()
      totalMagnetization += lattice.Magnetization.abs
   end

   # divido i valori medi per il numero di stati presi
   # in considerazione
   totalEnergy /= naverage.to_f
   totalMagnetization /= naverage.to_f

   # stampo i risultati
   printf("%.4f %.3f   %.3f     %.3f\n", beta,totalEnergy,lattice.ExactEnergy(beta)\
       ,totalMagnetization.abs)
end
print"\n\n"

