/*

   grid.cpp

   Description: 

   This code performs a short Monte-Carlo analysis to calculate a typical
   spin configuation for the Ising Model on a torus, and outputs the
   coordinates of all the "up" spins within a portion (100x100 sites) of
   the full lattice. The first file, "outputgrid.txt" is a direct output
   of the spins. For the second file, "outputgrid2.txt", the program
   performs a scale-factor 3 spin-block transformation, and outputs the
   coordinates of the "up" spin-blocks within a 300x300 site region. With
   the parameters set as they are now, it should take roughly 5 minutes to
   run. The time is linear in each of the parameters "N_iterations", "N",
   and "M". If N or M are set to values over 2000, the size of the lattice
   array must be increased.

   Compilation instructions:

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

                        g++ grid.cpp -o grid
                        ./grid

   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=2000; //this is the length in the periodic diection
static const int M=2000; //this is the length in the "infinite" direction
double J=.4406867935; //this is the reduced nearest-neighbor coupling
int N_iterations;


//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 OutfileGrid();
void OutfileGrid2();
void IdentifyPeriodic();


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

//main program

int main(){
  //initialize seed for random number generation
  srand ( time(0));
  InitializeLattice();
  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 IdentifyPeriodic(){
  for(int i=1; i<=N; i++){
    lattice[i][0].s=lattice[i][M].s;
    lattice[i][M+1].s=lattice[i][1].s;
  }
  for(int j=1; j<=M; j++){
    lattice[0][j].s=lattice[N][j].s;
    lattice[N+1][j].s=lattice[1][j].s;
  }
}

void OutfileGrid(){
  FILE *outfileGrid;
  outfileGrid=fopen("outputGrid.txt","w");
  for(int i=1; i<=100; i++){
    for(int j=1; j<=100;j++){
      if(lattice[i][j].s==1){
	fprintf(outfileGrid,"%d \t %d \n",i,j);
      }
    }
  }
}


void OutfileGrid2(){
  FILE *outfileGrid2;
  outfileGrid2=fopen("outputGrid2.txt","w");
  for(int i=1; i<=100; i++){
    for(int j=1; j<=100;j++){
      int k=(2+3*(i-1));
      int l=(2+3*(j-1));
      if((lattice[k-1][l-1].s+lattice[k][l-1].s+lattice[k-1][l].s+lattice[k][l].s+lattice[k][l+1].s+lattice[k+1][l].s+lattice[k+1][l+1].s+lattice[k-1][l+1].s+lattice[k-1][l+1].s)>0){
	fprintf(outfileGrid2,"%d \t %d \n",i,j);
      }
    }
  }
}