# include <ohana.h>

void   gauss_init ();
double gaussian (double x, double mean, double sigma);
double rnd_gauss (double mean, double sigma);

/* NGAUSS defines the resolution of the distribution */
# define NGAUSS 2048
static double gaussint[NGAUSS];

int main (int argc, char **argv) {

  int i, Npts;
  double mean, sigma, f;

  if (argc != 4) {
    fprintf (stderr, "USAGE: mkgauss (mean) (sigma) (npts)\n");
    exit (2);
  }

  mean = atof (argv[1]);
  sigma = atof (argv[2]);
  Npts = atoi (argv[3]);

  gauss_init ();

  for (i = 0; i < Npts; i++) {
    f = rnd_gauss (mean, sigma);
    fprintf (stdout, "%f\n", f);
  }

  exit (0);
}

void 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    , 0.0, 1.0) + 
	  9.0*gaussian(x+dx1, 0.0, 1.0) +
	  9.0*gaussian(x+dx2, 0.0, 1.0) + 
	  3.0*gaussian(x+dx3, 0.0, 1.0)) * (dx1/8.0);
    val += df;
    if (val > (i / (double) NGAUSS)) {
      gaussint[i] = x + dx / 2.0;
      i++;
    }
  }
  return;
}

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

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

double rnd_gauss (double mean, double sigma) {

  int i;
  double y;

  i = 0;
  while ((i == 0) || (i == NGAUSS - 1)) {
    y = drand48();
    i = MIN (NGAUSS - 1, MAX (0, NGAUSS*y));
  }
  y = gaussint[i]*sigma + mean;

  return (y);
}
