# include <ohana.h>
double drand48();
double gaussian();
double rnd_gauss();

# define NGAUSS 2048
double gaussint[NGAUSS];

double rnd_gauss (double mean, double sigma)
{
  
  double y;

  y = gaussint[(int)(NGAUSS*drand48())];
  
  return (y*sigma + mean);

}

gauss_init ()
{

  int i;
  long A, B;
  double val, x, dx, dx1, dx2, dx3, df;
  double mean, sigma;

  A = time(NULL);
  for (B = 0; A == time(NULL); B++);
  srand48(B);

  val = 0;
  dx = 0.001;
  dx1 = dx / 3.0;
  dx2 = 2.0*dx/3.0;
  dx3 = dx;
  mean = 0.0;
  sigma = 1.0;

  for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx)  {
    df = (3.0*gaussian(x    , mean, sigma) + 
	  9.0*gaussian(x+dx1, mean, sigma) +
	  9.0*gaussian(x+dx2, mean, sigma) + 
	  3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    val += df;
    if (val > (i / (double) NGAUSS)) {
      gaussint[i] = x + dx / 2.0;
      i++;
    }
  }
}

double gaussian (double x, double mean, double sigma)
{

  double f, X;
  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
  return (f);

}
