/* Metropolis Monte Carlo method for generating a Brownian motion
	in 1D with a fixed heatbath temperature and friction constant*/

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

#define  Nt1 100000     /* number of time-steps, transient... */
#define Nt 1000000      /* total number of time-steps, hystogram ...*/
#define h 1.38e-12         /* time-step in the integration of the
									 equation of motion      */
#define  T   300  /* temperature in K */
#define kB  1.38e-23 /* value of the Boltzmann constant */
#define x1 2e-4      /*interval for the random force */
#define mu 1           /* friction constant */
#define M 1e-11       /* mass of the particle */

long int H[24];
long int X[24];

/* ----------------------------------------------------------*/

double randf()
{return( (double)(rand())/(double)(RAND_MAX+1.0));
}

/* ----------------------------------------------------------*/

main()
{long int i;
 int q, s;
 double x, xn, xi, xf, w, r, v, sigma, vmin, vmax, xmin, xmax;
 double vmed, vmed2, xmed, xmed2, sigmav;

 sigma=mu*kB*T/h;
 sigmav=kB*T/M;

 v=0;
 vmin=0;
 vmax=0;
 xmin=0;
 xmax=0;

 for(s=0;s<=20; s++) {H[s]=0; X[s]=0;}

 x=randf()*(2*x1)-x1;
 printf("sigma=%g\n", sigma);
 getchar();
 for(i=1;i<=Nt1; i++)
	{
	 xn=randf()*(2*x1)-x1;
	 if (xn*xn>x*x) w=exp(-(xn*xn-x*x)/2.0/sigma);
	 else w=1;
	 r=randf();
	 if (r<=w) x=xn;
	 v=h/M*(x-mu*v)+v;
	 if (v<vmin) vmin=v;
	 if (v>vmax) vmax=v;
	 if (x<xmin) xmin=x;
	 if (x>xmax) xmax=x;
	}

 printf("xmin=%g     xmax=%g\n", xmin, xmax);
 printf("vmin=%g     vmax=%g\n", vmin, vmax);

  vmed=0;
  vmed2=0;
  xmed=0;
  xmed2=0;

  for(i=1;i<=Nt; i++)
	{
	 xn=randf()*(2*x1)-x1;
	 if (xn*xn>x*x) w=exp(-(xn*xn-x*x)/2.0/sigma);
	 else w=1;
	 r=randf();
	 if (r<=w) x=xn;
	 v=h/M*(x-mu*v)+v;
	 H[1+(int)((v-vmin)/(vmax-vmin)*20.0)]++;
	 X[1+(int)((x-xmin)/(xmax-xmin)*20.0)]++;
	 vmed=vmed+v;
	 vmed2=vmed2+v*v;
	 xmed=xmed+x;
	 xmed2=xmed2+x*x;
	}

	xmed=xmed/(double)(Nt);
	xmed2=xmed2/(double)(Nt);
	vmed=vmed/(double)(Nt);
	vmed2=vmed2/(double)(Nt);

	printf("xmed=%g      sigma=%g    xmed2=%g\n", xmed, sigma, xmed2);
	printf("vmed=%g      sigmav=%g   vmed2=%g\n", vmed, sigmav, vmed2);

	getchar();

	for(s=1; s<=20; s++) printf ("%ld     %ld\n", H[s], X[s]);
	getchar();


}

/* ---------------------------------------------------------*/