/********************************************************
 * Simulations in statistical physics - practice course *
 *                                                      *
 * Example: 2 dim. Ising model                          *
 *  * lattice size: L x L                               *
 *  * only nearest neighbor interaction (homogeneous)   *
 *  * homogeneous external field possible               *
 *  * periodic b.c.                                     *
 *                                                      *
 * Balogh Laszlo (blaci947@yahoo.co.uk)                 *
 *                                                      *
 * Budapest, Jan, 2011 (mod: Sep, 2011)                 *
 ********************************************************/

#include <math.h>
#include "downloadseed.h"

int main (int argc, char* argv[])
{
  // simulation control parameters:
  const int    L = 40;               // system size: L x L
  const double B = 0.0;              // external magnetic field
  const double T0 = 0.0, T1 = 4.0;   // temperature range: [T0;T1] NOTE: T0 < T1 is required
  const double dT = 0.1;             // temperature range resolution
  const int therm_steps = 1000;      // number of steps, which are ignored (to thermalize the system)
  const int every = 10;              // every ... steps are taken into account in the thermodynamic average
  const int avrg_steps = 10000;      // number of samples to the thermodynamic average
  
  // other parameters, variables:
  gsl_rng*                  rng = gsl_rng_alloc(gsl_rng_mt19937);  // Random Number Generator (Mersenne Twister type)
                                                                   // for more info, see: 
                                                                   //   http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html
  int    S[L][L];        // array for the sample: S[i][j] = 0 or 1 NOTE: this is not \sigma
  double sigma = 0.0;    // store \sum \sigma_i
  double sigma2 = 0.0;   // store ( \sum \sigma_i ) ^ 2
  int    n = 0;          // number of data in the thermodynamic averaging
  double PROB[2][5];     // pre-calculated acceptance ratios (see: documentation)
  double T;              // actual temperature [unit: J = 1]
  double nMC = floor((T1-T0+dT)/dT)*(therm_steps+every*avrg_steps); // number of MC scans (only for information)
  
  int i,j,k,s,sumj;
  //-----------------------------------------------------------------------------------------
  
  // initialize and "warm up" the RNG randomly:
  RNGSeedWarmup(rng);
    
  // initialize the configuration according to coin-flip:
  for (i=0; i<L; i++)
    for (j=0; j<L; j++) {
    S[i][j] = gsl_rng_uniform_int(rng, 2);
  }
  
  // before we start, save the settings to the output (file)
  printf("# system size:                %d x %d\n", L, L);
  printf("# external field:             %.5lE\n", B);
  printf("# temperature range:          %.12lf  -->  %.12lf\n", T1, T0);
  printf("# temperature resolution:     %.5lE\n", dT);
  printf("# simulation parameters:      therm_steps = %d; every = %d; avrg_steps = %d\n", therm_steps, every, avrg_steps);
  printf("# total number of MC scans:   %.3lE\n", nMC);
  printf("# estimated CPU time (@bell): %.2lf h\n", nMC*L*L*1.75E-11);
  fflush(stdout);
  
  // run the simulation over the temperature range:
  for (T=T1; T>=T0; T -= dT) {
    
    // reset thermodynamic averaging:
    n = 0;
    sigma = 0.0;
    sigma2 = 0.0;
    
    // pre-calculate the acceptance ratios for the possible spin-configurations:
    PROB[0][0] = exp(-8.0/T) * exp(2.0*B/T);
    PROB[0][1] = exp(-4.0/T) * exp(2.0*B/T);
    PROB[0][2] = 1.0         * exp(2.0*B/T);
    PROB[0][3] = exp(4.0/T)  * exp(2.0*B/T);
    PROB[0][4] = exp(8.0/T)  * exp(2.0*B/T);
    PROB[1][0] = exp(8.0/T)  * exp(-2.0*B/T);
    PROB[1][1] = exp(4.0/T)  * exp(-2.0*B/T);
    PROB[1][2] = 1.0         * exp(-2.0*B/T);
    PROB[1][3] = exp(-4.0/T) * exp(-2.0*B/T);
    PROB[1][4] = exp(-8.0/T) * exp(-2.0*B/T);
    
    // number of total steps:
    //   + therm_steps: thermalizing the system
    //   + (every * avrg_steps): simulation, but only every "every" steps are taken into account
    //     into the susceptibility calculation
    for (k=0; k < therm_steps + every * avrg_steps; k++) {
    
      // apply a "scan" of M.C. steps:
      
      for (i=0; i<L; i++)
        for (j=0; j<L; j++) {
	  
	  // sum the neighbouring spins of (i,j) 
	  // periodic boundary condition is implemented here:
	  sumj = 0;
	  sumj += S[(i+L-1)%L][ j ] + S[(i+1)%L][j] + S[i][(j+L-1)%L] + S[i][(j+1)%L];
	  
	  // change the (i,j) spin with the Metropolis-probability (pre-calculated):
	  if (gsl_rng_uniform(rng) < PROB[S[i][j]][sumj]) {
	    S[i][j] = 1 - S[i][j];
	  }
      }
    
      // get 1 sample of \sum \sigma_ij:
      //   (magnetization is normalized per site)
      //   (NOTE: \sigma is considered here: +1 or -1)
      if ( ( k>=therm_steps ) && ( (k-therm_steps)%every == 0 ) ) {
        s = 0;
        for (i=0; i<L; i++)
          for (j=0; j<L; j++) {
	    s += 2*S[i][j]-1;
        }
        sigma += abs( s );
        sigma2 += (double)s * (double)s;
        n++;
      }
    }// number of total steps
    
    // print results:
    //   1st col. = temperature
    //   2nd col. = magnetization
    //   3rd col. = susceptibility (calculated on the spot)
    printf("%18.12lf\t%15.10lf\t%15.6lE\n", T, sigma/n/L/L, (sigma2/n-sigma*sigma/n/n)/T/L/L); fflush(stdout);
  }
  
  return 0;
}