/*************************    Program 7.2    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

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

void mxwl (n,m,t,v)
/* This subroutine assigns velocities according to the Maxwell
   distribution. N is the total number of velocity components and
   M is the total number of degrees of freedom. T is the system
   temperature in the reduced units.  Copyright (c) Tao Pang 1997. */
int n,m;
double t;
double v[];
{
int i;
double v1,v2,ek,vs;
void grnf();

/* Assign a Gaussian distribution to each velocity component */

for (i = 0; i < n-1; i = i+2)
  {
  grnf(&v1,&v2);
  v[i] = v1;
  v[i+1] = v2;
  }

/* Scale the velocity to satisfy the partition theorem */

ek = 0;
for (i = 0; i < n; ++i)
  {
  ek = ek+v[i]*v[i];
  }
vs = sqrt(ek/(m*t));
for (i = 0; i < n; ++i)
  {
  v[i] = v[i]/vs;
  }
}

void grnf (x,y)
/* Two Gaussian random numbers generated from two uniform
   random numbers.  Copyright (c) Tao Pang 1997. */
double *x,*y;
{
double pi,r1,r2;
double ranf();

pi =  4*atan(1);
r1 = -log(1-ranf());
r2 =  2*pi*ranf();
r1 =  sqrt(2*r1);
*x  = r1*cos(r2);
*y  = r1*sin(r2);
}

double ranf()
/* Uniform random number generator x(n+1)= a*x(n) mod c
   with a = pow(7,5) and c = pow(2,31)-1.
   Copyright (c) Tao Pang 1997. */
{
const int ia=16807,ic=2147483647,iq=127773,ir=2836;
int il,ih,it;
double rc;
extern int iseed;
ih = iseed/iq;
il = iseed%iq;
it = ia*il-ir*ih;
if (it > 0)
  {
  iseed = it;
  }
else
  {
iseed = ic+it;
  }
rc = ic;
return iseed/rc;
}