#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*
 * Simple Ising model simulation program
 * Uses heat bath update
 * The program simulates Ising model at the given coupling
 * Prints out measurements of energy and magnetization
 * 
 * Needs files: ising_sim.c mersenne.h mersenne_inline.c
 * To compile:
 *   cc -O2 -o ising_sim ising_sim.c mersenne_inline.c -lm
 *
 * Asks parameters to be used
 */

#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(2))) : ");
  scanf("%lg",&beta);
  if (beta == 0.0)  
    beta = log(1+sqrt(2.0));  /* critical beta */
  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] = 1;

  /* 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()
{
  int sum,x,y;
  double p_plus,p_minus;


  /* sum over the neighbour sites - typewriter fashion */
  for (x=0; x<NX; x++) for (y=0; y<NY; y++) {
    sum = s[xup[x]][y] + s[xdn[x]][y] 
        + s[x][yup[y]] + s[x][ydn[y]];

    /* Heat bath update - calculate probability of +-1 */
    p_plus  = exp( (beta/2.0) * sum );   /* prob. of +1, unnormalized */
    p_minus = 1.0/p_plus;                 /* and -1 */
    p_plus  = p_plus/(p_plus + p_minus);  /* normalized prob of +1 */

    /* and choose the state appropriately */
    if (mersenne() < p_plus) s[x][y] = 1; else s[x][y] = -1;
  }
}


void measure(int iter,FILE *f)
{
  /* measure energy and magnetization, and print out
   */
  int sum,x,y,e,m;

  /* sum over the neighbour sites - typewriter fashion */
  e = m = 0;
  for (x=0; x<NX; x++) for (y=0; y<NY; y++) {
    e += 2 - s[x][y] * (s[xup[x]][y] + s[x][yup[y]]); 
    m += s[x][y];
  }
  
  fprintf(f,"%d %d %d\n", iter, e, m);
}