/*       Centre for Complexity Science, University of Warwick -- 2013       */
/* XY model simulation based on code for the Ising model by Charo del Genio */

/* compile with: gcc XY2d.c -o XY2d */
/* run with: ./XY2d [size] [temperature] [sweeps] [seed] [file1] [file2] */

/* BEGIN HEADER */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

/* global variables */
unsigned int STATE[32];
unsigned long int seed;

/* prototypes */
void parameters(int nump, char *param[], int *size, double *temperature, int *sweeps, FILE **p_f_output, FILE **spin_output); 
/* END HEADER */

//=============================
int main(int argc, char *argv[])
{
  char outputfile[80];
  register int i, j;
  int size, sweeps, x, y, square, t_init, t_av;
  double E, M, mx, my, temperature, accept, newspin, **spin;
  FILE *p_f_output, *spin_output;
  double boltzmann, diff;
  
  /* Check parameters */
  parameters(argc, argv, &size, &temperature, &sweeps, &p_f_output, &spin_output); 

  t_init = 10000; // equilibration time -- cheap trick only!
  
  square = size*size; // square lattice only 

  spin = (double **)malloc(sizeof(int *)*size); // allocate memory in a 2d array
  for(i=0; i<size; i++) {
    spin[i] = (double *)malloc(sizeof(int)*size);
  }

  for (i=0; i<size; i++) {
    for (j=0; j<size; j++) {
      spin[i][j] = 2.0*M_PI*drand48(); // Create initial state
    }
  }

  /* initial value of the energy and magnetisation */
  E = 0.0; M = 0.0; mx = 0.0; my = 0.0;
  for (i=0; i<size; i++) {
    for (j=0; j<size; j++) {
      E -= cos(spin[i][j] - spin[i==size-1 ? 0 : i+1][j]);
      E -= cos(spin[i][j] - spin[i][j==size-1 ? 0 : j+1]);
      mx += cos(spin[i][j]); my += sin(spin[i][j]);
    }
  }
  
  /* Simulate */
  for (i=1; i<=sweeps; i++) {  // for each sweep
    
    for (j=1; j<=square; j++) { // run over the lattice size
      x = size * drand48(); // choose a spin at random
      y = size * drand48();
      newspin = 2.0*M_PI*drand48(); // generate a new orientation at random      

      diff = 0.0;
      diff += cos(newspin - spin[x==0 ? size-1 : x-1][y]) - cos(spin[x][y] - spin[x==0 ? size-1 : x-1][y]);
      diff += cos(newspin - spin[x==size-1 ? 0 : x+1][y]) - cos(spin[x][y] - spin[x==size-1 ? 0 : x+1][y]);
      diff += cos(newspin - spin[x][y==0 ? size-1 : y-1]) - cos(spin[x][y] - spin[x][y==0 ? size-1 : y-1]);
      diff += cos(newspin - spin[x][y==size-1 ? 0 : y+1]) - cos(spin[x][y] - spin[x][y==size-1 ? 0 : y+1]);
      
      boltzmann = exp(diff/temperature);  

      if (drand48() < boltzmann) {
	mx += cos(newspin)-cos(spin[x][y]);
	my += sin(newspin)-sin(spin[x][y]);
	E -= diff;
	spin[x][y] = newspin;
      }
     
    } // end run over lattice size 

    // cheap way to calculate the average magnetisation -- the equilibration time should really be determined from the simulation
    if (i>t_init) { 
      t_av = i-t_init;
      M = M*((double) t_av-1)/((double) t_av) + sqrt(mx*mx+my*my)/((double) square)/((double) t_av);
    }
    
  } // end sweep
  
  // output the final spin configuration
  for (x=0; x<size; x++) { 
    for (y=0; y<size; y++) { 
      fprintf(spin_output, "%d\t%d\t%.10f\t%.10f\n",x,y,cos(spin[x][y]),sin(spin[x][y]));
    }
    //    fprintf(spin_output,"\n"); 
  } 
  fclose(spin_output); // close file  

  printf("The average magnetisation is %.6f\n",M);
  //  fprintf(p_f_output,"The average magnetisation is %.10f\n",M);
  
  free(spin[0]); // free memory
  free(spin);
  
  fclose(p_f_output); // close file
  
  return EXIT_SUCCESS;
} // end main
//==============================


/* Check user-provided parameters */
void parameters(int nump, char *param[], int *size, double *temperature, int *sweeps, FILE **p_f_output, FILE **spin_output)  
{	
  if (nump<7) {
    printf("\nUsage: %s [size] [prob] [temp] [sweeps] [seed] [output]\n\n",param[0]);
    printf("\t[size]   is the dimension of each side of the 2d Ising lattice.\n");
    printf("\t[temp]   is the temperature to simulate.\n");
    printf("\t[sweeps] is the number of sweeps to run.\n");
    printf("\t[seed]   is the random number generator seed (positive integer).\n");
    printf("\t[output1] is the output file name of the correlations.\n\n");
    printf("\t[output2] is the output file name of the spin configuration.\n\n");
    exit(EXIT_FAILURE);
  }
  
  *size = atoi(param[1]);
  if (*size<2) {
    printf("\nThe size of the model must be at least 2x2.\n\n");
    exit(EXIT_FAILURE);
  }
  
  *temperature = atof(param[2]);
  
  *sweeps = atoi(param[3]);
  if (*sweeps<1) {
    printf("\nYou must run at least 1 sweep.\n\n");
    exit(EXIT_FAILURE);
  }
  
  seed=strtoul(param[4],NULL,10);
  if (seed<=0) {
    printf("\nThe random number seed must be a positive integer.\n\n");
    exit(EXIT_FAILURE);
  }
  srand48(seed);
    
  if ( (*p_f_output = fopen(param[5],"r") ) != NULL ) {
    printf("\nError!\nOutput file %s already exists!\n\n",param[5]);
    exit(EXIT_FAILURE);
  }
  if ( (*p_f_output = fopen(param[5],"w") ) == NULL ) {
    printf("\nError creating output file %s.\n\n",param[5]);
    exit(EXIT_FAILURE);
  }
  
  if ( (*spin_output = fopen(param[6],"r") ) != NULL ) {
    printf("\nError!\nOutput file %s already exists!\n\n",param[6]);
    exit(EXIT_FAILURE);
  }
  if ( (*spin_output = fopen(param[6],"w") ) == NULL ) {
    printf("\nError creating output file %s.\n\n",param[6]);
    exit(EXIT_FAILURE);
  }
  
  return;
} // end parameters