#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAX_N 100 // max number of fleas
#define MAX_T 10000 // max number of timesteps
#define h 0.0001 //timestep

#define eps 1e-10 // increment for numerical derivative

// function prototipes
double rgauss(void);
double Vx(double); // potential
double derivative(double (*)(double), double);
double meansquaredisp(double [], double [], int);
double displacement2(double, double); 
void writeonfile(double [], double [], int);


// Parameters for potential
#define Depth 2.
#define sigma 0.5
#define req 0.


int main()

{

  int N;
  int t;
  int i,j;

  double x[MAX_N+1],x0[MAX_N+1];
  double xt,xt0; // temp variables

  double MSD;
  double dVx;
  double kBT,jump,noise;
  double dx;
  double rescape=5.*sigma; //escape distance
  int ninside,escaped;


  srand48(time(0));  
  
  do {
    printf("Input number of particles \n");
    scanf("%d",&N);
  } while (N<=0 || N > MAX_N);



  do {
    printf("Input energy intake (larger than 0 less or equal to 2)\n"); // (temperature)
    scanf("%lf",&kBT);
  } while (kBT < 0. || kBT > 5.);

  noise=sqrt(2.*kBT*h);

  
  // Evolve particles

  ninside=0;

  for (i=0;i<N;i++) {
    
    escaped=0; 
    t=0;

    xt=(double)lrand48()/(RAND_MAX+1.)*0.5;
	  if(lrand48()%2) xt=-xt;  // intial condition in (-0.5,0.5)
    xt0=xt; // store it

    do {
      t++;
      dVx=derivative(Vx,xt);
      jump=rgauss();
      xt = xt-dVx*h +noise*jump;

      //check if particle has escaped
      dx=fabs(xt);      
      if(dx>=rescape) {
	escaped=1;
	printf("Flea %d has escaped at time %d\n",i,t);
      }


    } while (t<MAX_T && escaped==0);

    // printf("Particle %d run time %d\n",i,t);


    if(!escaped) { // If particle has not escaped store its inital and final position
      	ninside++;
	x0[ninside]=xt0;
	x[ninside]=xt;
    }


  } // N particles


  // write on file initial and final positions of particles still inside
  writeonfile(x0,x,ninside);
  

  printf("\n%d particles over %d have escaped\n",N-ninside,N);
  printf("Fraction of escaped particles : %f\n",(double)(N-ninside)/N);
  
  //compute MSD for all particles still inside
  MSD=meansquaredisp(x0,x,ninside);
  printf("\nThe mean squared displacement for the particles stil inside is %f\n",MSD);

  
  return 0;
}


/******** FUNCTIONS DEFINITION  *****************/


double Vx(double r) {
  double vx;
  vx = Depth*( 1.-exp(-0.5*pow((r-req),2)/pow(sigma,2)) );
  return vx;
}


double derivative( double (*f)(double), double x){
  return (f(x+eps)-f(x))/eps;
}


double rgauss(void) {
  /* How to generate a Gaussian dist. 
     1) generate two random numbers r1 and r2 with uniform dist. between 0 and 1
     2) compute gamma=sin(2.*M_PI * r1) * sqrt(-2.*log(r2));  
        gamma is Gaussian with mu=0 and sigma=1
  */
  double g;
  double r1,r2=0.;
  r1=(double)lrand48()/(RAND_MAX+1.);
	while(r2==0){r2=(double)lrand48()/(RAND_MAX+1.);}
  g= sin(2.*M_PI * r1) * sqrt(-2.*log(r2));
  return g;
}


double displacement2(double r0, double r) {
  return pow((r-r0),2.);
}


double meansquaredisp(double r0[], double r[], int ndata) {
  int i;
  double r2;

  r2=0.;
  for (i=0;i<ndata;i++){
    r2+=displacement2(r0[i],r[i]);
  }
  
  r2=r2/ndata;

  return r2;
}



void writeonfile(double r0[], double r[], int ndata) {

  FILE *fp;
  int i;
  
  if((fp = fopen("data.dat", "w")) == NULL) {
    printf("Cannot open file.\n");
    exit(1);
  }

  for (i=0;i<ndata;i++) {
    fprintf(fp,"%4d %10.6f\n",i+1,sqrt(displacement2(r0[i],r[i])));
  }

  if( fclose( fp ))
    printf("File close error.\n");

}