//    Simple program to compute the energy for the Potts class of models
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include "lib.h"
using namespace std;

ofstream ofile;


void ReadData(int&,double&,double&,double&,int&,int&);
void Initialize(int, int **, double *);
void Energy(double,double *);
void Thermalization(int,double *,int **,long&,double&,int);
void Metropolis(int,double *,int **,long&,double&);
void WriteAverages(double *, double, int);


const int J = 1;
int N          ;

inline int periodic(int i, int limit, int add) { 
  return (i+limit+add) % (limit);
}

int main(void){
  
  int MCs, nTermal, **Spin, q;
  double T0 , T1, dT,  E , *Boltzmann, *Average ;    
  long seed = -1;
  
  ofile.open("potts.txt");
  ReadData( q, T0, T1, dT, MCs,nTermal );
  
  // Vectors
  Average = new(nothrow) double [5]          ;
  Boltzmann = new(nothrow) double [4]      ;
  Spin = (int**) matrix(N, N, sizeof(int));  
  if(!Average||!Boltzmann){
    cout <<"Error in Boltzmann........."       ;
    exit(EXIT_FAILURE)                     ;
  }
  
  Initialize(q, Spin, Average );   

  for( double T = T0; T <= T1; T += dT){    
    Energy( T, Boltzmann)                               ; 
    Thermalization(q, Boltzmann, Spin, seed, E, nTermal); 
    for(int i = 1; i <= MCs; i++){
      Metropolis(q,Boltzmann,Spin,seed, E); 
      Average[0] += E      ; 
      Average[2] += E*E    ;
    }
    // write results 
    WriteAverages(Average, T, MCs );    
    //Initialize for new temperature
    for(int i=0; i<5 ;i++){
      Average[i]=0;
    }
  }
 ofile.close();
 free_matrix((void **) Spin);

}// main


void ReadData(int& q,double& T0,double& T1,double& dT,int& MCs,int& nTermal){

    cout << "-------Potts model      ------------- " << endl   ;
    cout << "-------------------------------------- " << endl   ;
    cout << "Start temperature (KT [eV]) :"                      ;
    cin  >> T0                                                  ;
    cout << "End temperature (KT [eV]) :"                      ;
    cin  >> T1                                                  ;
    cout << "Temperature interval (KT [eV]) :"                   ;
    cin  >> dT                                                  ;
    cout << " Monte Carlo cycles :"                     ;
    cin  >> MCs                                                 ;
    cout << "q for the Potts model :"                            ;   
    cin  >> q                                                   ;
    cout << "Lattice size N (total square lattice N x N):"                            ;
    cin  >> N                                                   ;
    cout << "Monte Carlo termalisation cycles" << endl;
    cin  >> nTermal                                             ;  

}//ReadData


void Initialize(int q, int **Spin, double *Average){
  
  long seed = -4;
    for(int i = 0; i < N; i++){
      for(int j = 0; j < N; j++){
        //Spin[i][j] = 1;
        Spin[i][j] = (int)(ran1(&seed)*(q)+1); 
      }
    } 
    for(int i=0; i<5 ;i++){
      Average[i]=0;
    }

}//Initialize

void Energy(double T,double *Boltzmann){

  Boltzmann[0] = exp(-J/T)  ;
  Boltzmann[1] = exp(-2*J/T);
  Boltzmann[2] = exp(-3*J/T);
  Boltzmann[3] = exp(-4*J/T);
    
}//Energy  
 

void Thermalization(int q,double *Boltzmann,int **Spin,long& seed,double& E, int nTermal){ 

  int  SpinFlip, LocalEnergy0, LocalEnergy, x, y, dE;
    
    for(int i = 0; i < N; i++){
      for(int j = 0; j < N; j++){
        LocalEnergy0 = 0;
        LocalEnergy = 0;
        dE = 0;
        x = (int) (ran1(&seed)*N);
        y = (int) (ran1(&seed)*N);
        if(Spin[x][y] == Spin[x][periodic(y,N,-1)])
          LocalEnergy0 --;
        if(Spin[x][y] == Spin[periodic(x,N,-1)][y])
          LocalEnergy0 --;
        if(Spin[x][y] == Spin[x][periodic(y,N,1)])
          LocalEnergy0 --;
        if(Spin[x][y] == Spin[periodic(x,N,1)][y]) 
          LocalEnergy0 --;
        
        do{                                                        
        SpinFlip = (int)(ran1(&seed)*(q)+1);               
        }while(SpinFlip == Spin[x][y]);

        if(SpinFlip == Spin[x][periodic(y,N,-1)])
          LocalEnergy --;
        if(SpinFlip == Spin[periodic(x,N,-1)][y])
          LocalEnergy --;
        if(SpinFlip == Spin[x][periodic(y,N,1)])
          LocalEnergy --;
        if(SpinFlip == Spin[periodic(x,N,1)][y]) 
          LocalEnergy --;
        
        dE =  (LocalEnergy - LocalEnergy0); 

        if(dE<=0){
          Spin[x][y] = SpinFlip;
        }
        else if(ran1(&seed)<Boltzmann[dE-1]){
          Spin[x][y] = SpinFlip;
        }
      }
    }
  E=0;
  for(int x =0; x < N; x++) {
    for (int y= 0; y < N; y++){
      if(Spin[x][y] == Spin[periodic(x,N,-1)][y])
        E -= J;
      if(Spin[x][y] == Spin[x][periodic(y,N,-1)])
        E -= J;
        }
  }
 
}//Thermalization    


void Metropolis(int q,double *Boltzmann,int **Spin,long& seed,double& E){ 

  int  SpinFlip, LocalEnergy0, LocalEnergy, x, y, dE;
   
    for(int i = 0; i < N; i++){
      for(int j = 0; j < N; j++){
        x = (int) (ran1(&seed)*N);
        y = (int) (ran1(&seed)*N);
        LocalEnergy0 = 0;
        LocalEnergy = 0;
        dE = 0;
        if(Spin[x][y] == Spin[x][periodic(y,N,-1)])
          LocalEnergy0 --;
        if(Spin[x][y] == Spin[periodic(x,N,-1)][y])
          LocalEnergy0 --;
        if(Spin[x][y] == Spin[x][periodic(y,N,1)])
          LocalEnergy0 --;
        if(Spin[x][y] == Spin[periodic(x,N,1)][y]) 
          LocalEnergy0 --;
        
        do{                                                        
        SpinFlip = (int)(ran1(&seed)*(q)+1);               
        }while(SpinFlip == Spin[x][y]);

        if(SpinFlip == Spin[x][periodic(y,N,-1)])
          LocalEnergy --;
        if(SpinFlip == Spin[periodic(x,N,-1)][y])
          LocalEnergy --;
        if(SpinFlip == Spin[x][periodic(y,N,1)])
          LocalEnergy --;
        if(SpinFlip == Spin[periodic(x,N,1)][y]) 
          LocalEnergy --;
        
        dE =  LocalEnergy - LocalEnergy0; 

        if(dE<=0){
          Spin[x][y] = SpinFlip;
          E += J*dE;
        }
        else if(ran1(&seed)<Boltzmann[dE-1]){
          Spin[x][y] = SpinFlip;
          E += J*dE;
        }
      }
    }
 
}//Metropolis

void WriteAverages(double *Average , double T, int MCs ){

  double norm = 1/(double)(MCs);
  double EAverage    = Average[0]*norm;
  double E2Average   = Average[2]*norm;
  double EVariance   = (E2Average - EAverage*EAverage)/N/N;
    
  ofile << setiosflags(ios::showpoint | ios::uppercase);
  ofile << setw(15) << setprecision (8) << T;
  ofile << setw(15) << setprecision (8) << EAverage/N/N;
  ofile << setw(15) << setprecision (8) << EVariance/T/T;
  ofile << setw(15) << setprecision (8) << E2Average<<endl;
 
} // WriteAverage