/*  
   	compile using "icpc -c phi4_2D_v2.cpp"
   	compile using "icpc -c twister.cpp"
   	compile using "icpc phi4_2D_v2.o twister.o -o phi4_2D_v2.exe"

   	./phi4_2D_v2.exe 
	-L (L) -beta (beta) -lambda (lambda) 
	-cal_start (cal_start) -cal_stop (cal_stop) -cal_step (cal_step) 
	-run (run)
*/
#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstring>
#include <ctime>
#include <fstream>
#include <sstream>
#include <time.h>
#include "twister.h"

void samp(double ** phi, long L, double beta, double lambda, double *m, double *e);
void config_out(long L, double beta, double lambda, double **phi, long run);

int main (int argc, char *argv[]) {

	long L;     /* The sidelength of the array */
	double beta;
	double lambda;

	long i,j,a; /* loop counters */

	/* Observables */
	double m = 0.0;
	double m_1st = 0.0;    /* average magnetization */
	double m_1st_abs = 0.0;
	double m_2nd = 0.0;
	double m_4th = 0.0;

	double e=0.0;    
	double e_1st=0.0; /* average energy per spin */

	long cal_start;///平衡状態に達したと思われるステップ
	long cal_stop;
	long cal_step;
	long run_start;
	double delta=3.0;
	long trajectory=0;

    for(i=1;i<argc;i++){
        if(!strcmp(argv[i],"-L")){////strcmp関数は文字列argv[i]と"-for"が一致したら0を返す
            sscanf(argv[++i],"%ld",&L);
        }
        else if(!strcmp(argv[i],"-beta")){
            sscanf(argv[++i],"%lf",&beta);
        }
        else if(!strcmp(argv[i],"-lambda")){
            sscanf(argv[++i],"%lf",&lambda);
        }
        else if(!strcmp(argv[i],"-cal_start")){
            sscanf(argv[++i],"%ld",&cal_start);
        }
        else if(!strcmp(argv[i],"-cal_stop")){
            sscanf(argv[++i],"%ld",&cal_stop);
        }
        else if(!strcmp(argv[i],"-cal_step")){
            sscanf(argv[++i],"%ld",&cal_step);
        }
        else if(!strcmp(argv[i],"-run_start")){
            sscanf(argv[++i],"%ld",&run_start);
        }
        else if(!strcmp(argv[i],"-delta")){
            sscanf(argv[++i],"%lf",&delta);
        }
        else if(!strcmp(argv[i],"-trajectory")){
            sscanf(argv[++i],"%ld",&trajectory);
        }
	}

    //long seed=1994;
    //init_genrand(seed);
    init_genrand((unsigned long)(time(NULL)+run_start));///時間を用いた乱数の初期化


	/* Output some initial information */
  	printf("L=%ld, beta=%lf, lambda=%lf, delta=%lf, cal_start=%ld, cal_stop=%ld, cal_step=%ld, run_start=%ld\n",L,beta,lambda,delta,cal_start,cal_stop,cal_step,run_start);

	char file_name[256];
	sprintf(file_name,"phi_vs_time_L%ld_beta%lf_lambda%lf_delta%lf_%ld_%ld_%ld_run_start%ld.txt",L,beta,lambda,delta,cal_start,cal_step,cal_stop,run_start);
  	FILE *fp_output=fopen(file_name,"w");

	char file_name_moments[256];
	sprintf(file_name_moments,"moments_L%ld_%ld_%ld_%ld_run_start%ld.txt",L,cal_start,cal_step,cal_stop,run_start);
  	FILE *fp_output_moments=fopen(file_name_moments,"a");

  	long N=L*L;

  	/* Allocate memory for the system */
  	double **phi=(double**)malloc(L*sizeof(double*));
  	for(i=0;i<L;i++){ 
		phi[i]=(double*)malloc(L*sizeof(double));
	}

  	/* Generate an initial state */
	for(i=0;i<L;i++){
		for(j=0;j<L;j++){
	      	phi[i][j]=0.0;
	    }
	}


	long run=run_start;

  	long num_points = 0;
  	for(long MCstep=0; MCstep<cal_stop; MCstep++){

    	for(a=0;a<N;a++){//one MC step is N-trials
      	/* randomly select a spin */
      		i=genrand_int32()%L;
      		j=genrand_int32()%L;
			double phi_old=phi[i][j];
			phi[i][j]=phi[i][j]+delta*(genrand_real1()-0.5);//trial
			double phi2_new=phi[i][j]*phi[i][j];
			double phi2_old=phi_old*phi_old;

			double delta_u = -beta*(phi[i][j]-phi_old)*( phi[i?(i-1):(L-1)][j] + phi[(i+1)%L][j] + phi[i][j?(j-1):(L-1)] + phi[i][(j+1)%L] )
							 + (phi2_new + lambda*(phi2_new-1)*(phi2_new-1)) - (phi2_old + lambda*(phi2_old-1)*(phi2_old-1));
      		if( genrand_real1() < exp(-delta_u) ){//
				//printf("Accept\n");
      		}else{
				//printf("Reject\n");
				phi[i][j]=phi_old;
			}
    	}
    	/* Sample and accumulate averages */
    	if( (MCstep>=cal_start) && (MCstep%cal_step==0) ){
      		samp(phi,L,beta,lambda,&m,&e);
      		m_1st += m;
      		m_1st_abs += fabs(m);
      		m_2nd += m*m;
      		m_4th += m*m*m*m;
      		e_1st += e;
      		num_points++;
            fprintf(fp_output,"%ld\t%lf\t%lf\n",MCstep,m,e);
            //printf("%ld\t%lf\t%lf\n",MCstep,m,e);
			if(trajectory==1){
				config_out(L,beta,lambda,phi,run);
			}
			run+=1;
    	}
 	}

	m_1st/=num_points;
	m_1st_abs/=num_points;
	m_2nd/=num_points;
	m_4th/=num_points;

	//fprintf(fp_output_moments,"%lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\n",beta,m_1st,m_1st_abs,(L*L)*m_2nd,beta*(L*L)*(m_2nd-m_1st_abs*m_1st_abs),m_4th,0.5*(3.0-m_4th/(m_2nd*m_2nd)));			

	fclose(fp_output);
	fclose(fp_output_moments);


    for(i=0;i<L;i++){
	    free(phi[i]);
    }
	free(phi);
}

void samp(double **phi, long L, double beta, double lambda, double *m, double *e ){

	long i,j;
	*m=0.0;
  	*e=0.0;
  	/* Visit each position (i,j) in the lattice */
  	for(i=0;i<L;i++){
    	for(j=0;j<L;j++){
      		*m+=phi[i][j];
      		*e+= -beta*(phi[i][j])*(phi[i][(j+1)%L]+phi[(i+1)%L][j]) + phi[i][j]*phi[i][j] + lambda*(phi[i][j]*phi[i][j]-1)*(phi[i][j]*phi[i][j]-1);
    	}
  	}
  	*m/=(L*L);
  	*e/=(L*L);
}
void config_out(long L, double beta, double lambda, double **phi, long run){


	char file_name_config[256];
	sprintf(file_name_config,"config_L%ld_beta%lf_lambda%lf_run%ld.txt",L,beta,lambda,run);
  	FILE *fp_output_config=fopen(file_name_config,"w");

	for(long i=0;i<L;i++){
		for(long j=0;j<L;j++){
            fprintf(fp_output_config,"%ld\t%ld\t%lf\n",i,j,phi[i][j]);			
		}
	}
	fclose(fp_output_config);
}
