///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Maxwell.java                                       //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: July 7, 2006                                           //
//                                                                       //
// (1) This Java program is part of 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.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

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

  public static double[] maxwell(double m, double T,
    int N) {
    int nv = 3*N;
    int ng = nv-6;
    double v[] = new double[nv];
    double r[] = new double[2];

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

// Scale the velocity to satisfy the partition theorem
    double ek = 0;
    for (int i=0; i<nv; ++i) ek += v[i]*v[i];
    double vs = Math.sqrt(m*ek*nv/(ng*T));
    for (int i=0; i<nv; ++i) v[i] /= vs;
    return v;
  }

// Method to create two Gaussian random numbers from two
// uniform random numbers in [0,1].

  public static double[] rang() {
    double x[] = new double[2];
    double r1, r2;
    r1 = - Math.log(1-ranf());
    r2 = 2*Math.PI*ranf();
    r1 = Math.sqrt(2*r1);
    x[0] = r1*Math.cos(r2);
    x[1] = r1*Math.sin(r2);
    return x;
  }

// Method 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.

  public static double ranf() {
    final int a = 16807, c = 2147483647, q = 127773,
      r = 2836;
    final double cd = c;
    int h = seed/q;
    int l = seed%q;
    int t = a*l-r*h;
    if (t > 0) seed = t;
    else seed = c+t;
    return seed/cd;
  }