/*****************  gaussian_ran.c ****************
 *  double gaussian_ran()
 *  Gaussian distributed random number
 *  Probability distribution exp( -x*x/2 ), so < x^2 > = 1
 *  Uses mersenne random number generator
 */
#include <stdlib.h>
#include <math.h>
#include "mersenne.h"

double gaussian_ran()
{
  static int iset=0;
  static double gset;
  register double fac,r,v1,v2;

  if  (iset) {
    iset = 0;
    return(gset);
  } 

  do {
    v1 = 2.0*mersenne() - 1.0;
    v2 = 2.0*mersenne() - 1.0;
    r  = v1*v1 + v2*v2;
  } while (r >= 1.0 || r == 0.0);
  fac  = sqrt( -2.0*log(r)/r );
  gset = v1*fac;
  iset = 1;
  return(v2*fac);
}