/*

monte_carlo.cpp

Description:

 This code performs a Monte-Carlo analysis to calculate a
typical spin configuration for the Ising Model on a torus (5000x500 spins)
and then computes correlation functions in the longitudinal direction. It
outputs the resulting correlation function (as a function of lattice
distance) in the file "correlations.txt". Currently, the number of Monte
Carlo iterations is set to 50,000. This is enough to produce an accurate
critical correlation function over 200+ lattice site spacings. It should
take roughly 5 hours of CPU time on 2008-grade PCs. Run time is roughly
linear in N, M and NUM_ITER. If N or M is increased above its present
value, the size of the spin array must be increased.

Compilations instructions:

Strip off the html commands and store the source in the file
monte_carlo.cpp. Under gcc then execute

                        g++ monte_carlo.cpp -o monte_carlo
                        ./monte_carlo

Copyright 2008 Douglas Stanford. The author grants permission to copy,
distribute and display this work in unaltered form, with attribution to
the author, for noncommercial purposes only. All other rights, including
commercial rights, are reserved to the author.

*/

#include 
#include 
#include 
#include 
#include 


//variables
static const int N=500; //this is the length in the periodic diection
static const int M=5000; //this is the length in the "infinite" direction
static const int NUM_ITER=50000; //This is the number of Monte Carlo
			     //iterations performed
double J=.4406867935; //this is the critical reduced nearest-neighbor coupling
time_t start_time;
time_t end_time;


//functions
double H(); //this function returns the energy
void InitializeLattice(); //this function initializes the spins randomly
void MonteCarlo(); //this function applies the Monte Carlo process to each spin, in turn
void ComputeCorrelations(int height, int width);
double Distance(int i1,int i2,int j1,int j2);
void IdentifyPeriodic();


//structs
struct lattice_point{
  double s;
}lattice[510][5010];
double correlations[1010];
int counter[1010];

//main program

int main(){
  //initialize seed for random number generation
  srand ( time(0));
 
  //open output files
  
  FILE *outfile2;
  outfile2=fopen("correlations.txt","w"); 
  
  //set up the lattice of spins randomly
  InitializeLattice();

  //Iterate Monte Carlo to realize a typical configuration
   for(int i=1; i0.5){
	lattice[i][j].s=1;
      }
      else{
	lattice[i][j].s=-1;
      }
    }
  }
}


void MonteCarlo(){
  //NOTE: this MonteCarlo iteration assumes toroidal topology as well
  double delta=0;
  int j1=0,j2=0,i1=0,i2=0;
  for(int i=1; i<=N; i++){
    i1=i-1;
    i2=i+1;
    if(i1==0){i1=N;}
    if(i2==M+1){i2=1;}
    for(int j=1; j<=M; j++){
      j1=j-1;
      j2=j+1;
      if(j1==0){j1=M;}
      if(j2==M+1){j2=1;}
      delta=J*lattice[i][j].s*(lattice[i][j1].s+lattice[i][j2].s+lattice[i1][j].s+lattice[i2][j].s);
      if(delta<0){
	lattice[i][j].s=-lattice[i][j].s;
      }
      else if(-2*delta>log((double)rand()/((double)RAND_MAX))){
	lattice[i][j].s=-lattice[i][j].s;
      }
    }
  }
}

void ComputeCorrelations(int height, int width){
  for(int i=1; i<=N; i++){
    for(int j1=501; j1<=M; j1++){
      for(int j2=(j1-500); j2