/* Monte Carlo simulation of 1-D Boltzmann model */
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#define SEED 171

main ()
{

  /* introduce constants and variables */
  float alpha, beta, betaprime ,div,v_value, largest, smallest, bin_occ;  
  int j,k,l,n,m,t,count,b_value; 
  long rnumber;			/* randum number */
  float vb = 1.0;		/* bath velocity */
  float vsum;
  float v2sum;
  float sigma;
  float sum;
  float a = 0.01;		/* mass ratio */
  int Nmax = 1000000;		/* size of ensemble */
  float v[Nmax];		/*  array of particle velocities */
  int Tmax = 64;			/* maximimum run time */
  int i;			/* parity indicator */
  int Bmax = 100;                  /* number of sort bins */
  long bin[Bmax+1]; 
  long *binptr;   /* bins to store counts */
  /* compute alpha and beta */
  alpha = (1.0 - a) / (1.0 + a);
  betaprime = (2.0 * a * vb) / (1.0 + a);
/* initialize the random number generator */
srand(SEED);

/* intialize the arrays */
  binptr=&bin[0];

  for (n = 0; n < Nmax; n++)
    {
      v[n] = 0.0;
    }
 for (n = 0; n <= Bmax; n++)
    {
      *(binptr+n)=0;
      bin[n]=0;
    } 

/* begin the run */
  for (m = 0; m < Nmax; m++)
    {
      for (t = 1; t <= Tmax; t++)
	{
	  rnumber = rand();
	  if (2 * (rnumber / 2) != rnumber)
	    v[m] = alpha * v[m] - betaprime;
	  else
	    v[m] = alpha * v[m] + betaprime;
	}
    }

/* sort the data and put it into bins */
  smallest = v[0];
  largest = v[0];
  for (j = 0; j < Nmax; j++)
    {
      if (v[j] > largest)
	largest = v[j];
      if (v[j] < smallest)
	smallest = v[j];
    }

/* create 500 bins */
  div = (largest - smallest) / Bmax;
  for (k = 0; k < Nmax; k++)
    {
      b_value=floor((v[k]-smallest)/div);
      *(binptr+b_value)=*(binptr+b_value)+1;
    }
  vsum=0.0;
  v2sum=0.0;
    for (count = 0; count <= Bmax; count++)
    {
      v_value = smallest + (float) count * div;
      bin_occ=(float)(*binptr)/(div*(float)Nmax);
      vsum=vsum+v_value*bin_occ;
      v2sum=v2sum+v_value*v_value*bin_occ;
      printf("%f\t%f\n", v_value , bin_occ);
      binptr++;
    } 				/* end of data processing section */

sigma=sqrt(v2sum-vsum*vsum);
printf("# standard deviation is %f\n",sigma);

}				/* end of program */