/* compile line ; gcc -o ensemble ensemble.c -lm */
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#define SIZE 30000 /* ensemble size */

/* create a macro for generating a sign from two rands */
#define sgn(a,b) (((a)>= (b)) ? (1.0) : (-1.0))

/* parameters for computations */
int BINS =200;
float DIV ,distribution;
int Tmax=210;

/* create a timestep */
float dt=0.02;


/* define a few useful variables */
int n,m,k,l;
long I,J,b_value;
float bin_x;



main()
{
/* create some arrays to store x, v */
float bins[BINS+1], *binptr;
float x[SIZE],v[SIZE];
float *xptr, *vptr, kx1,kx2,kv1,kv2;

/* Set up a data file */
FILE *fpt;
fpt=fopen("data", "w+");


DIV=2.0/BINS;
/* initialize the random number generator */
srand((unsigned long)getpid());

/* initialize the positions and speeds */
for(n=0;n<SIZE;n++){
I=rand();
J=rand();
x[n]=sgn(I,J)*((float)rand()/(float)RAND_MAX); /* random initial x */
I=rand();
J=rand();
v[n]=sgn(I,J)*(1.0-x[n]*x[n]);                 /* corresponding v */
}

xptr=&x[0];
vptr=&v[0];

for (n = 0; n <= BINS; n++)
    {
      bins[n]=0;
    } 
binptr=&bins[0];


/* start a time loop, Runge-Kutta integration */
for(n=0;n<SIZE;n++){
for(m=0;m<Tmax;m++){
kx1=*(vptr+n)*dt;
kv1=-(*(xptr+n))*dt;
kx2=(*(vptr+n)+0.5*kv1)*dt;
kv2=-(*(xptr+n)+0.5*kx1)*dt;
*(xptr+n)=(*(xptr+n))+kx2;
*(vptr+n)=(*(vptr+n))+kv2;
  }
}

/* create a frequency plot from the x-values */
  for (k = 0; k < SIZE; k++)
    {
      b_value=floor((*(xptr+k)+1.0)/DIV);
      *(binptr+b_value)=*(binptr+b_value)+1;
    }
  for (l= 0; l <= BINS; l++)
    {
      bin_x=-1.0+(float)l*(float)DIV;
      distribution=(*binptr)*(float)BINS/(float)SIZE;
      fprintf(fpt,"%f\t%f\n", bin_x, distribution);
      binptr++;
    }                         /*   end of data processing section */

fclose(fpt);

}