/* compile with:   g++ observables_periodic_sinus_potential.cpp nrutil.cpp random_mars.cpp -o observables_periodic_sinus_potential -O3 -w */

//./observables_periodic_sinus_potential 55536349 0.0001 100 100. 10.0
// ./observables_periodic_sinus_potential seed dt Nsim oder Multiple sim, freq, V0, test
//Header ********************************************************************************************************

//parameters
double d1; 	
double d2; 	 
double L; 		//Boxlenght
double Tinit;
double mu;
long N;
double N_double;
class RanMars *random_mars;

double dx; 
double dt; 		//timestep for integration
int Neq;
long Nsim;
double B;

int seed;
using namespace std; //for the execution time measurement

double *x;		
double *x_nc;		
double *px;

double beta_T;

//functions
void set_pos_vel();
void calculate_potential(int j);
void mc_step();
void pbc(double &x);

double x_before;
double px_before;
double px_after;

double v0; 
double v1;

//MSD Variables
double x_save; //save positions
double x_save_nc; 
double delta_x; 
double delta_x0;
double  delta_x_nc;

int N_blocks; 
int N_levels; 
int count_MSD;
long MSD_frequency;
int multiple_Nsim; 

/*END Header************************************************************************ */

#include <stdio.h> 
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string>
#include <sstream>
#include <complex>
#include <iostream>
#include "nrutil.h"
#include "random_mars.h"

#include <bits/stdc++.h>  //for the execution time measurement

#define PI 3.1415926535897
const std::complex<double> i(0, 1);
std::complex<double> isf;

/*********************************************************************************************************************************/

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

	printf("Main working \n");
	clock_t start_eq, end_eq,start_sim, end_sim; //for the execution time measurement

	seed=atoi(argv[1]);
	random_mars = new RanMars(seed);	// initialize random number generator
	srand (time(NULL));// initialize random number generator
		
	//N=atoi(argv[4]); 	//numer of particles I average 
	N=1;
	x = dvector(0,N-1);
	x_nc = dvector(0,N-1);
	px = dvector(0,N-1);
	L=1;
	d1=0.5;
	d2=L-d1;
 	dt = atof(argv[2]);//intgration time step	

	//for the bd_step
	Tinit = 1.0; //temperature in units of epsilon/kB.
	mu = 1.0;  //mobility
	beta_T= 1/(double) Tinit;
	double D=Tinit* mu;
	B = sqrt(2.0*D*dt);

	double Q=2.0;
	//potentials
	v0 = atof(argv[5]); 
	v1 = 0.;
		
	set_pos_vel();// initialize particles,velocity,forces,...

	//FILE * out = fopen("data_MC_position.dat","w");
	N_blocks = 10; // 10; // number of blocks -> how many steps I want to make for each level
	N_levels = 4; // number of levels -> how many orders of magnitude 
	MSD_frequency=atol(argv[4]);

	//Nsim = (int) pow(N_blocks,N_levels)*MSD_frequency;  //minimum number of timesteps
	multiple_Nsim=atoi(argv[3]) ;
	Nsim = (long) pow(N_blocks,N_levels)*MSD_frequency*multiple_Nsim;
	printf("Nsim %ld \n ", Nsim);
	Neq=100.;


	start_eq = clock();
	for (long step=0; step<Neq; step++) {
		mc_step();
	} 
	end_eq = clock(); 

	start_sim = clock(); 
	// declare MSD array for blockings
	double blocking_sumx[N_blocks][N_levels]; 
	double blocking_sumx_nc[N_blocks][N_levels];	 
	double count_block[N_blocks][N_levels]; //for counting how often sth is added to the msdblocking and devide tho get the averade
	double msd_blocking[N_blocks][N_levels];
    //noice cancelling
	double msd_blocking_nc[N_blocks][N_levels];
	double msd_blocking_nc_cc[N_blocks][N_levels];
	double msd_blocking_isf[N_blocks][N_levels];
	// initialize array with zeros
	for (int i=0; i<N_blocks; i++){
		for (int j=0; j<N_levels; j++){
			blocking_sumx[i][j]=0.0;
			blocking_sumx_nc[i][j]=0.0;
	
			msd_blocking[i][j] =0.0;				
			msd_blocking_nc[i][j]=0.0;	
			msd_blocking_nc_cc[i][j]=0.0;
			msd_blocking_isf[i][j] =0.0;	

			count_block[i][j] =0.0;
		}
	}
	long step_msd;
	x_save = x[0];
	x_save_nc= x_nc[0];

	for (long step=1; step<Nsim; step++) {

		mc_step(); // perform one mc step	

		if (step%MSD_frequency==0){ 
			step_msd=step/MSD_frequency;
			
			delta_x = (x[0]-x_save); 
			//delta_x_nc = delta_x-(x_nc[0]-x_save_nc); //old version
			delta_x_nc = (x_nc[0]-x_save_nc); 

			x_save = x[0];
			x_save_nc = x_nc[0];


			double delx;
			double delx_nc;

			int j0; // = step % N_blocks;
			int max_levels = (int) round(log(step_msd)/log(N_blocks));
			
			if (max_levels>=N_levels) {
				max_levels=N_levels-1;
				} //Important, as otherwise, when k == 4, and so on, data will be written to memory where none exists, leading to the overwriting of count_msd, for example.

			for (int k=0; k<(max_levels+1); k++) {
				
				if (k==0) { //delta... what to add to the previous blocksum to get the new one. for k=0 its always the length of 1 timestep, for k=1 its length of 10 timesteps, for k=2 its the length of 100 timesteps
					delx = delta_x; //the delta I add for k==0, so for the level 0 is always just vx=x[t]-x[t-1] /dt so for J0=0/k=0==v. for J0=1/k=0==J0=0/k=0+v for J0=2/k=0==J0=1/k=0+v ....
					delx_nc = delta_x_nc; 

				} else {
					delx = blocking_sumx[N_blocks-1][k-1];  //the delta I add for level k>0, so for level 1,2,3... its always the last element (N_blocks-1) of the level before [k-1]. for k=1 J0=0/k=1==J0=9/k=0, for J0=1/k=1==J0=0/k=1+J0=9/k=0,.....		
					delx_nc = blocking_sumx_nc[N_blocks-1][k-1];  //the delta I add for level k>0, so for level 1,2,3... its always the last element (N_blocks-1) of the level before [k-1]. for k=1 J0=0/k=1==J0=9/k=0, for J0=1/k=1==J0=0/k=1+J0=9/k=0,.....

				}

				int nblocks_to_k = 1 ; // (int) round(pow(N_blocks,k)); //for calculating j0
				for (int kk=0; kk<k; kk++){
					nblocks_to_k *= N_blocks;
				}
					
				if (step_msd % nblocks_to_k == 0 ) {

					j0 =  ((step_msd) / nblocks_to_k-1 ) % N_blocks;  // was -1

					if (j0==0) {
						blocking_sumx[j0][k] = delx; //blocking_sum[N_blocks-1][k-1][n] ; for the first element in the block take just the last element from the level before	
						blocking_sumx_nc[j0][k] = delx_nc; //blocking_sum[N_blocks-1][k-1][n] ; for the first element in the block take just the last element from the level before

					} else {
						blocking_sumx[j0][k] = blocking_sumx[j0-1][k] + delx; //for elements j0>0 in the block take value from previous element of the block +delta
						blocking_sumx_nc[j0][k] = blocking_sumx_nc[j0-1][k] + delx_nc; //for elements j0>0 in the block take value from previous element of the block +delta
					}

					double dr2x = blocking_sumx[j0][k]*blocking_sumx[j0][k]; //calculate mean-square-displacement
					double dr2x_nc = blocking_sumx_nc[j0][k]*blocking_sumx_nc[j0][k];
					double dr2x_nc_cc  = blocking_sumx_nc[j0][k]*blocking_sumx[j0][k];
					
					msd_blocking[j0][k] +=  dr2x;
					msd_blocking_nc[j0][k] +=  dr2x_nc;
					msd_blocking_nc_cc[j0][k] += dr2x_nc_cc;
					
					isf=exp(i*Q*(double)blocking_sumx[j0][k]);  
					msd_blocking_isf[j0][k] += std::real(isf);
					count_block[j0][k] += 1.0; //for taking the average (over particles and number of timesteps
				} 
			}
		}
	} 
	//fclose(out);
	end_sim = clock(); 

	FILE * out_msd = fopen("data_observables_periodic_sinus_potential.dat","w");
	fprintf(out_msd, "Montecarlo_stufenpotential_periodic Parameters: dt:%f Tinit:%f mu:%f N_blocks:%d N_levels:%d MSD_frequency:%d multiple_Nsim:%d Nsim:%d Neq:%d seed:%d  d2:%f d1:%f\n", dt, Tinit, mu,N_blocks,N_levels,MSD_frequency,multiple_Nsim,Nsim, Neq, seed, d2, d1);
	
	//Fehler/Varianz des MSD irgendwann benötigt?
	for (int i=0; i<N_levels; i++){
		for (int j=0; j<N_blocks; j++){ //-1 because else the end block and the beginning of the new block has the same time sampled
			double time_here = (j+1) * (pow(N_blocks,i))*dt*MSD_frequency;
			double msd_here = msd_blocking[j][i]/count_block[j][i]; 
			double msd_here_nc = msd_blocking_nc[j][i]/count_block[j][i];
			double msd_here_nc_cc = msd_blocking_nc_cc[j][i]/count_block[j][i];
			double msd_here_isf = msd_blocking_isf[j][i]/count_block[j][i]; 
			fprintf(out_msd, "%d %d %lf %lf %.15lf %.15lf %.15lf %.15lf \n",i, j, count_block[j][i], time_here, msd_here, msd_here_nc, msd_here_nc_cc, msd_here_isf);	
		}
	}
	fclose(out_msd);
	//fclose(out);
	end_sim = clock();     // Calculating total time taken by the program. 
    printf("time needed for equilibration %f \n", double(end_eq - start_eq) / double(CLOCKS_PER_SEC));
    printf("time needed for main simulation %f \n", double(end_sim - start_sim) / double(CLOCKS_PER_SEC));	

}
/*********************************************************************************************************************************/

void set_pos_vel(){
	for (long i=0; i<N; i++) {
		//double randi_xi=random_mars->gaussian();
		x[i]=L/4.0;
		x_nc[i]=0.0;
		calculate_potential(i);
		//px[i]=v0;
 	}

	//monte carlo 
	x_before=0.0;
	px_before=0.0;
	px_after=0.0;

	//for the MSD
	delta_x = 0.0;
	delta_x_nc = 0.0;

	x_save = 0.0;
	x_save_nc= 0.0;

}

/*********************************************************************************************************************************/

void calculate_potential(int j) {	

		px[j]= 0.5*v0* std::cos(2*PI*x[j]);

}

/*********************************************************************************************************************************/

void mc_step(){  
	 for (long i=0; i<N; i++) {

		double randi_x=random_mars->gaussian();
		x_before=x[i];
		px_before=px[i];
		x[i]+=B*randi_x;
		//pbc(x[i]);

		calculate_potential(i);
		px_after=px[i];

		double deltaE= px_after-px_before;
		if (deltaE> 0.0){
			double randi_exp= random_mars->uniform();

			if  (exp ( - beta_T *deltaE) < randi_exp) {
				x[i]=x_before;
				px[i]=px_before;
				x_nc[i]+=-B*randi_x;		
			}
		}
	}

}

void pbc(double &x){
  if (x >= L/2.0) x -= L;
  if (x < -L/2.0) x += L;
}