#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>

using namespace std;


/* Program for 3-state Potts model modified from Ising model program available
 * in course webpage (ising_sim.c). Uses heat bath update and measures
 * energy and magnetization.
 *
 *Changed 29.11.2007: Fixed energy measurement.
 *
 *
*/

extern "C"
{
    #include "mersenne.h"
};

/* Lattice size, adjust */
#define NX 64
#define NY 64

/* How many iterations for each measurements */
#define N_MEASURE 10

/* spin storage */
static int s[NX][NY];
static double beta;

/* Neighbour index arrays, to be filled at the beginning
 */
static int xup[NX],yup[NY],xdn[NX],ydn[NY];


void update(void);
void measure(int iter, FILE *f);

int main(int argc, char* argv[])
{
  int i,x,y,n_loops;
  long seed;
  char fname[200];
  FILE *f;

  /* Read in the input */
  printf(" Beta (value 0: beta_c = log(1+sqrt(3))) : ");
  scanf("%lg",&beta);
  if (beta == 0.0)
    beta = log(1+sqrt(3.0));  /* critical beta for 2D 3-state Potts model*/
  printf(" Number of update sweeps : ");
  scanf("%d",&n_loops);

  printf(" Name of measurement file : ");
  scanf("%s",fname);
  f = fopen(fname,"w");       /* open output file */

  printf(" Random number seed : ");
  scanf("%ld",&seed);
  seed_mersenne( seed );

  /* "Warm up" the rng generator */
  for (i=0; i<543210; i++) mersenne();

  printf(" ++++++++++++++++++++++++++++++++++++++++++\n");
  printf(" Ising model, %d x %d lattice, beta %g\n",NX,NY,beta);
  printf(" Heat bath update, %d updates/measurement, %d updates\n",
	 N_MEASURE, n_loops);
  printf(" Output file %s\n",fname);
  printf(" Random seed %ld\n", seed );


  /* fill up the index array */
  for (i=0; i<NX; i++) {
    xup[i] = (i+1) % NX;
    xdn[i] = (i-1+NX) % NX;
  }

  for (i=0; i<NY; i++) {
    yup[i] = (i+1) % NY;
    ydn[i] = (i-1+NY) % NY;
  }

  /* Initial spins - completely ordered */
  for (x=0; x<NX; x++) for (y=0; y<NY; y++) s[x][y] = 0;

  /* and the update/measure loop */
  for (i=0; i<n_loops; i++) {
    update();
    if (i % N_MEASURE == 0) measure(i,f);
  }

  fclose(f);

  printf(" ** simulation done\n");

  return(1);
}


void update()
//-----------
{
  /* sum over the neighbour sites - typewriter fashion */
    for (int x=0; x<NX; x++)
    {
        for (int y=0; y<NY; y++)
        {
            char M[3] = {0,0,0};
            M[s[xup[x]][y]]++;
            M[s[xdn[x]][y]]++;
            M[s[x][yup[y]]]++;
            M[s[x][ydn[y]]]++;

            const double E0 = -M[0];
            const double E1 = -M[1];
            const double E2 = -M[2];

            const double p0 = exp(beta * -E0);
            const double p1 = exp(beta * -E1);
            const double p2 = exp(beta * -E2);

            const double rand = (p0+p1+p2) * mersenne();
            if(rand < p0)
                s[x][y] = 0;
            else
            {
                if(rand < p0+p1)
                    s[x][y] = 1;
                else
                    s[x][y] = 2;
            }
        }
    }
}



void measure(int iter,FILE *f)
//----------------------------
{
  /* measure energy and magnetization, and print out
   */
    int x,y,e;
    double M[3] = {0,0,0};


    /* sum over the neighbour sites - typewriter fashion */

    e = 0;
    for (x=0; x<NX; x++)
    {
        for (y=0; y<NY; y++)
        {
	        if(s[x][y] != s[xdn[x]][y]) e++;
	        if(s[x][y] != s[x][yup[y]]) e++;
            M[s[x][y]]++;
        }
    }

    const double Mabs = 1.5 * (*max_element(M, M+3)/(NX*NY) - 1.0/3.0);

    fprintf(f,"%d %d %lg\n", iter, e, Mabs);
}