def main

  nthermal=20
  naverage=50
  nsweep=10
  betaini  = 0.420
  betafin  = 0.480
  betastep = 0.005

  # creo un array di dimensione NSIZE*NSIZE e con i valori uguali ad 1
  # 1 è lo stesso oggetto per tutti
  spin = Array.new(NSIZE*NSIZE,1)

  # inizializzo il generatore di numeri pseudocasuali
  myrand 108767

  # assegno i valori iniziali agli spin
  startSpin(spin)

  puts "\nReticolo degli spin "+NSIZE.to_s+"x"+NSIZE.to_s
  puts "Passaggi di termalizzazione "+nthermal.to_s
  puts "Media su "+naverage.to_s+" stati"
  puts "calcolata su stati dopo "+nsweep.to_s+" passate"

  # termalizzazione preliminare
  Metropolis(spin, betaini,5*nthermal)
  beta = betaini
  while beta < betafin
    # termalizzazione
    Metropolis(spin, beta,nthermal)

    # azzero i valori medi
    totalEnergy = 0.0
    totalMagnetization = 0.0
    for k in 1..naverage
      Metropolis(spin, beta, nsweep)
      # Calcolo dei valori indipendenti dell'energia
      # e della magnetizzazione
      # Quindi li sommo.
      totalEnergy += Energy(spin)
      totalMagnetization += Magnetization(spin)
    end
    # prendo i valori medi
    totalEnergy /= (naverage.to_f)
    totalMagnetization /= (naverage.to_f)
    # stampo i risultati
    t = ExactEnergy(beta)
    printf("\n beta: %8f   E: %8f  E esatta: %8f  Magnetizzazione %8f",\
           beta,totalEnergy,t,totalMagnetization.abs)
    #incremento beta
    beta =  beta + betastep

  end
  puts " "

end

def startSpin (spin)

  # Scelgo la direzione iniziale degli spin a caso

  spin.each_index do |i|
    if myrand < 0.5
      spin[i] = -1
    end
  end

end

def Energy(spin)
  #
  # contributo di ogni legame all'energia. Ogni legame
  # viene contato una sola volta anziche' due,
  # per cui il risultato finale va moltiplicato per due


  # stampa(spin)

  isum = 0.0
  for i in 0...NSIZE
    for j in 0...(NSIZE-1)
      k = i*NSIZE+j
      isum += spin[k]*spin[k+1]
    end
    isum += spin[i*NSIZE+NSIZE-1]*spin[i*NSIZE]
  end
  for j in 0...NSIZE
    for i in 0...(NSIZE-1)
      k = i*NSIZE + j
      isum += spin[k]*spin[k+NSIZE]
    end
    isum += spin[(NSIZE-1)*NSIZE+j]*spin[j]
  end
  isum=2*isum

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

end

def Magnetization(spin)

  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.abs.to_f/(NSIZE*NSIZE).to_f

end

def Metropolis(spin,beta,nsweep)
  # esegue n "passate" con l'algoritmo di Metropolis su tutto il reticolo

  somma = spin.inject {|sum,q| sum += spin[q]}

  f4=Math.exp(-8*beta)    # k =  4
  f2=Math.exp(-4*beta)    # k =  2

  # loop per NSIZE spin, fatto nsweep volte
  for m in 1..nsweep
    # loop sulle righe
    for i in 0...NSIZE
      ip1=(i+1)%NSIZE
      im1=i-1
      if(i==0)
        im1 = NSIZE-1
      end
      # loop sulle colonne
      for j in 0...NSIZE
        jp1=(j+1)%NSIZE
        jm1=j-1
        if j==0
          jm1=NSIZE-1  # invece di NSIZE-1
        end
        k=spin[ip1*NSIZE+j]+spin[im1*NSIZE+j]+\
          spin[i*NSIZE+jp1]+spin[i*NSIZE+jm1]
        deltaE= 2*spin[i*NSIZE+j]*k
        # Differenza di energia negativa, accetto il cambiamento
        if deltaE <= 0
          spin[i*NSIZE+j] = -spin[i*NSIZE+j]
        # Differenza di energia positiva, il cambiamento
        # viene accettato con probabilità exp(-beta*deltaE)
        else
          case k.abs # valore assoluto di delta spin, può essere solo 2 o 4
            when 4
            if f4 > myrand
              spin[i*NSIZE+j] = -spin[i*NSIZE+j]
            end
            when 2
            if f2 > myrand
              spin[i*NSIZE+j] = -spin[i*NSIZE+j]
            end
          end
        end
      end
    end
  end
end

def ExactEnergy(beta)

   q1=2*Math.sinh(2*beta)/(Math.cosh(2*beta)**2)
   q2=2*(Math.tanh(2*beta))**2-1.0

   h=M_PI/2.0/1000.0

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

   sum=sum-0.5*(1.0+1.0/Math.sqrt(1.0-q1*q1))
   sum=sum*h

   return -2.0/Math.tanh(2*beta)*(1.0+2*q2/M_PI*sum)

end

def myrand(icall=0)

# generatore random predefinito in Ruby
   if icall != 0
      srand(icall)
   end

   return rand

end

def stampa (spin)

  puts " "

  for i in 0...(NSIZE*NSIZE)
    if i%NSIZE == 0
      puts "\n"
    end
    if spin[i] == 1
      printf("+ ")
    else
      printf("- ")
    end
  end
  puts "\n\n"

end

M_PI = 3.14159265
NSIZE = 50

main