/*************************************************************************/
/* Program file name: maxwell.c                                          */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: July 7, 2006                                           */
/*                                                                       */
/* (1) This C program is created for the book, "An Introduction to       */
/*     Computational Physics, 2nd Edition," written by Tao Pang and      */
/*     published by Cambridge University Press on January 19, 2006.      */
/*                                                                       */
/* (2) No warranties, express or implied, are made for this program.     */
/*                                                                       */
/*************************************************************************/

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

// Function to draw initial velocities fron the Maxwell
// distribution for a given mass m, temperature T, and
// particle number N.

void maxwell (N, m, T, v)
int N;
double m, T;
double v[nv];
{
  int i, nv = 3*N, ng = nv-6;
  double v1, v2, ek, vs;
  void rang();

// Assign a Gaussian number to each velocity component
  for (i=0; i< nv-1; i+=2) {
    rang(&v1, &v2);
    v[i] = v1;
    v[i+1] = v2;
  }

// Scale the velocity to satisfy the partition theorem
  ek = 0;
  for (i=0; i<nv; ++i) ek += v[i]*v[i];
  vs = sqrt(m*ek*nv/(ng*T));
  for (i=0; i<nv; ++i) v[i] = v[i]/vs;
}
void rang (x, y)
// Two Gaussian random numbers generated from two uniform
// random numbers.
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() {
// Function to generate a uniform random number in [0,1]
// following x(i+1)=a*x(i) mod c with a=pow(7,5) and
// c=pow(2,31)-1.  Here the seed is a global variable.

  const int a = 16807,  c = 2147483647, q = 127773, r = 2836;
  int l, h, t;
  double cd = c;
  extern int seed;

  h = seed/q;
  l = seed%q;
  t = a*l - r*h;
  if (t > 0) seed = t;
  else seed = c + t;
  return seed/cd;
}