/* Metropolis Monte Carlo update of 2-d Ising ferromagnet */

#include <math.h>

extern int size;

void
metropolis(sweeps)
int sweeps;
{
   int imet,i,j,ipos,jpos,ineg,jneg;
   int old_spin,new_spin,spin_sum;
   int old_energy,new_energy;
   double energy_diff;
   double ranf();  /* random number generator */

   /* Loop over sweeps */
   for (imet=0;imet<sweeps;imet++) {

      /* Loop over sites */
      for (i=0;i<size;i++) {
      for (j=0;j<size;j++) {

         /* periodic boundary conditions */
         ipos = (i+1) % size;
         jpos = (j+1) % size;
         ineg = (i+size-1) % size;
         jneg = (j+size-1) % size;

	 old_spin = spin[i][j];
	 new_spin = - old_spin;

	 /* Sum of neighboring spins */
	 spin_sum = spin[i][jpos] + spin[ipos][j] + 
		    spin[i][jneg] + spin[ineg][j];

	 old_energy = - old_spin * spin_sum;   
	 new_energy = - new_spin * spin_sum;
	 energy_diff = beta * (new_energy - old_energy);   /* beta = J/kT */

	 /* Do the update */
	 if ( (energy_diff ><= 0.0) || (exp(-energy_diff) > ranf() ) ) {
	    /* Accept the change */
	    spin[i][j] = new_spin;
	 }

      }}  /* end of loop over sites */

   }  /* end of loop over sweeps */

}