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

//./trajectory 55536349 0.01 0 10000. 1.0
// ./trajectory seed dt protential, freq, V0
//potential 0 for barrier potential
//potential 1 for sinus potential
//********************************************************************************************************

//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 pot; 

/*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=2.0*PI;
	d1=PI;
	d2=PI;
 	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);
	//potentials
	v0 = atof(argv[5]); 
	v1 = 0.;

	set_pos_vel();// initialize particles,velocity,forces,...
	FILE * out = fopen("../data/data_trajectory.dat","w");
	MSD_frequency=atol(argv[4]);
	pot=atoi(argv[3]) ;
	// Nsim = atol(argv[3]) ;  //number of timesteps
	Nsim = (long) pow(N_blocks,N_levels)*MSD_frequency;
	printf("Nsim %ld \n ", Nsim);
	Neq=0.;

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

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

	for (long step=1; step<Nsim; step++) {
		fprintf(out,"%d %f %f %f \n",step,x[0],x_nc[0],x[0]-x_nc[0] );
		mc_step(); // perform one mc step	
	} 

	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++) {

		x[i]=L/4.0;
		x_nc[i]=L/4.0;
		calculate_potential(i);

	}
	
	//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) {	
	if (pot==0) {
		px[j]  = v0;
		if(remainder(x[j], L)>=0.&& remainder(x[j], L)<d2){
		  px[j] = v1;	  
		}
    } else {
        px[j]= v0/2* std::cos( 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;
		x_nc[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;
			}
		}
	}

}

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