Index: trunk/Ohana/src/addstar/src/mkcmf.c
===================================================================
--- trunk/Ohana/src/addstar/src/mkcmf.c	(revision 39241)
+++ trunk/Ohana/src/addstar/src/mkcmf.c	(revision 39242)
@@ -13,6 +13,4 @@
 # define FLAGS 0x1101
 
-void gauss_init (int Nbin);
-double rnd_gauss (double mean, double sigma);
 void writeStars_PS1_V5_Lensing (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
 void writeStars_PS1_V5 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
@@ -282,5 +280,5 @@
   }
     
-  gauss_init (2048);
+  gaussdev_init ();
 
   // load test stars from a file:
@@ -564,78 +562,4 @@
 
 
-static int Ngaussint = 0;
-static double *gaussint;
-
-extern double drand48();
-
-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);
-
-}
-
-/* integrate a gaussian from -5 sigma to +5 sigma */
-void gauss_init (int Nbin) {
- 
-  int i;
-  double val, x, dx, dx1, dx2, dx3, df;
-  double mean, sigma;
- 
-  /* no need to generate this if it already exists */
-  if (Ngaussint == Nbin) return;
-
-  // A = time(NULL);
-  // // XXX this is expensive if called a lot (1 sec min)
-  // // for (B = 0; A == time(NULL); B++);
-  // B = A + 10000;
-  // srand48(B);
- 
-  Ngaussint = Nbin;
-  ALLOCATE (gaussint, double, Ngaussint + 1);
-
-  val = 0;
-  dx = 1.0 / Ngaussint;
-  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 < Ngaussint) && (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 + 0.5) / (double) Ngaussint) {
-      gaussint[i] = x + dx / 2.0;
-      i++;
-    }
-  }
-}
-
-double rnd_gauss (double mean, double sigma) {
- 
-  int i;
-  double y;
- 
-  y = drand48();
-  i = Ngaussint*y;
-  y = gaussint[i]*sigma + mean;
- 
-  return (y);
- 
-}
- 
-double int_gauss (int i) {
-  double y;
-  y = gaussint[i];
-  return (y);
-}
- 
 void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars) {
 
@@ -654,7 +578,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -700,7 +624,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -749,7 +673,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -800,7 +724,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -857,7 +781,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -918,7 +842,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -987,7 +911,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -1094,7 +1018,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -1198,7 +1122,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN * rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN * gaussdev_rnd(0.0, 1.0);
       flux = pow(10.0, -0.4 * M[i]);
       fSN = 1.0 / sqrt(flux);
@@ -1324,7 +1248,7 @@
 
     if (ADDNOISE) {
-      X[i] += FX * fSN * rnd_gauss(0.0, 1.0);
-      Y[i] += FY * fSN * rnd_gauss(0.0, 1.0);
-      M[i] += fSN*rnd_gauss(0.0, 1.0);
+      X[i] += FX * fSN * gaussdev_rnd(0.0, 1.0);
+      Y[i] += FY * fSN * gaussdev_rnd(0.0, 1.0);
+      M[i] += fSN*gaussdev_rnd(0.0, 1.0);
       flux = pow (10.0, -0.4*M[i]);
       fSN = 1.0 / sqrt(flux);
Index: trunk/Ohana/src/fakeastro/Makefile
===================================================================
--- trunk/Ohana/src/fakeastro/Makefile	(revision 39241)
+++ trunk/Ohana/src/fakeastro/Makefile	(revision 39242)
@@ -38,5 +38,4 @@
 $(SRC)/fakestar_io.$(ARCH).o \
 $(SRC)/insert_fakestar.$(ARCH).o \
-$(SRC)/gaussian.$(ARCH).o \
 $(SRC)/fakeastro_images.$(ARCH).o \
 $(SRC)/fakeastro_images_region.$(ARCH).o \
Index: trunk/Ohana/src/fakeastro/src/make_2mass_measures.c
===================================================================
--- trunk/Ohana/src/fakeastro/src/make_2mass_measures.c	(revision 39241)
+++ trunk/Ohana/src/fakeastro/src/make_2mass_measures.c	(revision 39242)
@@ -80,6 +80,6 @@
 
     // uR,uD in linear arcsec
-    double dRsee = rnd_gauss (0.0, 1.0 / SN);
-    double dDsee = rnd_gauss (0.0, 1.0 / SN);
+    double dRsee = gaussdev_rnd (0.0, 1.0 / SN);
+    double dDsee = gaussdev_rnd (0.0, 1.0 / SN);
 
     double dRoff = (dRpm + dRsee) / 3600.0;
Index: trunk/Ohana/src/fakeastro/src/make_fake_stars_catalog.c
===================================================================
--- trunk/Ohana/src/fakeastro/src/make_fake_stars_catalog.c	(revision 39241)
+++ trunk/Ohana/src/fakeastro/src/make_fake_stars_catalog.c	(revision 39242)
@@ -113,6 +113,6 @@
 
     // uR,uD in linear arcsec
-    double dRsee = rnd_gauss (0.0, 1.0 / SN);
-    double dDsee = rnd_gauss (0.0, 1.0 / SN);
+    double dRsee = gaussdev_rnd (0.0, 1.0 / SN);
+    double dDsee = gaussdev_rnd (0.0, 1.0 / SN);
 
     double dRoff = (dRpm + dRsee) / 3600.0;
Index: trunk/Ohana/src/fakeastro/src/make_fakeqsos.c
===================================================================
--- trunk/Ohana/src/fakeastro/src/make_fakeqsos.c	(revision 39241)
+++ trunk/Ohana/src/fakeastro/src/make_fakeqsos.c	(revision 39242)
@@ -51,5 +51,5 @@
 
     // Mr will be interpretted as m_r
-    double Mr = rnd_gauss (18.0, 1.5);
+    double Mr = gaussdev_rnd (18.0, 1.5);
     
     stars[i].R = R;
Index: trunk/Ohana/src/fakeastro/src/make_fakestars.c
===================================================================
--- trunk/Ohana/src/fakeastro/src/make_fakestars.c	(revision 39241)
+++ trunk/Ohana/src/fakeastro/src/make_fakestars.c	(revision 39242)
@@ -48,5 +48,5 @@
       int inPatch = FALSE;
       while (!inPatch) {
-	z = rnd_gauss (0.0, FAKEASTRO_ZGAL);
+	z = gaussdev_rnd (0.0, FAKEASTRO_ZGAL);
 	r = sqrt(drand48()) * FAKEASTRO_RGAL;
 	Lrad = drand48() * 2 * M_PI;
@@ -98,7 +98,7 @@
     double Mr;
     if (bigPeak) {
-      Mr = rnd_gauss (11.25, 1.0);
+      Mr = gaussdev_rnd (11.25, 1.0);
     } else {
-      Mr = rnd_gauss (10.00, 2.5);
+      Mr = gaussdev_rnd (10.00, 2.5);
     }
     
Index: trunk/Ohana/src/libohana/Makefile
===================================================================
--- trunk/Ohana/src/libohana/Makefile	(revision 39241)
+++ trunk/Ohana/src/libohana/Makefile	(revision 39242)
@@ -35,4 +35,5 @@
 $(SRC)/sorts.$(ARCH).o		 \
 $(SRC)/bisection.$(ARCH).o       \
+$(SRC)/gaussdev.$(ARCH).o       \
 $(SRC)/string.$(ARCH).o		 \
 $(SRC)/findexec.$(ARCH).o	 \
Index: trunk/Ohana/src/libohana/include/ohana.h
===================================================================
--- trunk/Ohana/src/libohana/include/ohana.h	(revision 39241)
+++ trunk/Ohana/src/libohana/include/ohana.h	(revision 39242)
@@ -486,3 +486,9 @@
 int sprintf_double (char *output, double value);
 
-# endif
+/* in gaussdev.c */
+double gaussian (double x, double mean, double sigma);
+void gaussdev_init (void);
+void gaussdev_free (void);
+double gaussdev_rnd (double mean, double sigma);
+
+# endif
Index: trunk/Ohana/src/libohana/src/gaussdev.c
===================================================================
--- trunk/Ohana/src/libohana/src/gaussdev.c	(revision 39242)
+++ trunk/Ohana/src/libohana/src/gaussdev.c	(revision 39242)
@@ -0,0 +1,83 @@
+# include <ohana.h>
+
+// NOTE: this code uses a local static array which must be freed on exit 
+static int Ngaussint = 0;
+static double *gaussint = NULL;
+
+// extern double drand48();
+
+// return the value of a gaussian at position x with mean,sigma
+double gaussian (double x, double mean, double sigma) {
+
+  double f;
+
+  double S2 = sigma*sigma;
+
+  f = exp (-0.5 * SQ(x - mean) / S2) / sqrt(2 * M_PI * S2);
+
+  return (f);
+}
+
+/* integrate a gaussian from -10 sigma to +10 sigma */
+// how long does this take with a reasonably high resolution?
+void gaussdev_init (void) {
+ 
+  int i;
+  double val, x, dx, dx1, dx2, dx3, df;
+  double mean, sigma;
+ 
+  if (gaussint) return;
+
+  // 16k bins takes ~0.1 sec to generate on a ~3GHz core and gives good resolution
+  Ngaussint = 0x4000;
+  ALLOCATE (gaussint, double, Ngaussint + 1);
+
+  val = 0;
+  dx = 1.0 / Ngaussint;
+  dx1 = dx / 3.0;
+  dx2 = 2.0*dx/3.0;
+  dx3 = dx;
+  mean = 0.0;
+  sigma = 1.0;
+ 
+  for (i = 0, x = -10.0; (i < Ngaussint) && (x <= 10.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 + 0.5) / (double) Ngaussint) {
+      gaussint[i] = x + dx / 2.0;
+      myAssert (i < Ngaussint + 1, "oops");
+      i++;
+    }
+  }
+}
+
+void gaussdev_free (void) {
+  free (gaussint);
+}
+
+// return a number drawn from a gaussian deviate with mean, sigma
+// must first call gaussdev_init()
+double gaussdev_rnd (double mean, double sigma) {
+ 
+  int i;
+  double y;
+ 
+  myAssert (gaussint, "need to call gaussdev_init before calling gaussdev");
+
+  y = drand48();
+  i = Ngaussint*y;
+  y = gaussint[i]*sigma + mean;
+ 
+  return (y);
+}
+ 
+// return the value of \int gaussian(x,mean,sigma) 
+// double gaussian_integral (double x, double mean, double sigma) {
+double gaussian_integral (int i) {
+
+  // int i = (x - mean) / (sigma * Ngaussint);
+  return gaussint[i];
+}
Index: trunk/Ohana/src/misc/src/fakedist.c
===================================================================
--- trunk/Ohana/src/misc/src/fakedist.c	(revision 39241)
+++ trunk/Ohana/src/misc/src/fakedist.c	(revision 39242)
@@ -1,9 +1,3 @@
 # include <ohana.h>
-
-# define NGAUSS 2048
-double drand48();
-double gaussint[NGAUSS];
-double gaussian();
-double rnd_gauss();
 
 # define MMIN 1.0
@@ -45,5 +39,5 @@
   init_outmatrix (); 
  
-  gauss_init ();
+  gaussdev_init ();
 
   read_datafiles ();
@@ -349,6 +343,6 @@
 	if (noise < dMo) {
 	  *Noise = noise;
-	  *mag = rnd_gauss (v, noise);
-	  *color = rnd_gauss (uv, 1.4*noise);
+	  *mag = gaussdev_rnd (v, noise);
+	  *color = gaussdev_rnd (uv, 1.4*noise);
 	  return;
 	}
@@ -414,67 +408,4 @@
 
 /*****************************************************************************/
-
-double
-rnd_gauss (mean, sigma)
-double mean, sigma; 
-{
-
-  int i;
-  double y;
-
-  y = drand48();
-  i = NGAUSS*y;
-  y = gaussint[i]*sigma + mean;
-
-  return (y);
-
-}
-
-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++;
-    }
-  }
-  return (x + dx / 2.0);
-}
-
-double 
-gaussian (x, mean, sigma) 
-double x, mean, sigma; 
-{
-
-  double f, X;
-
-  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
-
-  return (f);
-
-}
 
 /*****************************************************************************/
@@ -505,6 +436,6 @@
       Ngood = 1.0;
       for (k = 0; k < NTRY; k++) {
-	v = rnd_gauss (mag, noise);
-	uv = rnd_gauss (color, 1.4*noise);
+	v = gaussdev_rnd (mag, noise);
+	uv = gaussdev_rnd (color, 1.4*noise);
 	x = (uv - UV0 - 0.7*A_V) / DUV;
 	y = (v - V0 - Dist - A_V) / DV;
Index: trunk/Ohana/src/misc/src/fakepop.c
===================================================================
--- trunk/Ohana/src/misc/src/fakepop.c	(revision 39241)
+++ trunk/Ohana/src/misc/src/fakepop.c	(revision 39242)
@@ -2,18 +2,7 @@
 # define MMIN 1.0
 # define MMAX 120.0
-extern double drand48();
 extern double rnd_mass();
 double rnd_mass();
 double term;
-extern double rnd_gauss();
-extern double rnd_integrate ();
-extern double gaussian ();
-
-double gaussian();
-double rnd_gauss();
-double rnd_integrate();
-
-# define NGAUSS 2048
-double gaussint[NGAUSS];
 
 double AgeS, AgeE, alpha, d, Av, dMo, dVo, dVref;
@@ -47,5 +36,5 @@
   gfits_scan (&UV_h, "DEC_Y", "%lf", 1, &ldM);
   
-  gauss_init ();
+  gaussdev_init ();
 
   fprintf (stderr, "beginning main loop\n");
@@ -68,6 +57,6 @@
 	v = gfits_get_matrix_value (&V_i, X, Y) + d + Av; 
 	noise = dVo*sqrt(1.0 + pow (10.0, (0.4*(v - dVref)))); 
-	V = rnd_gauss (v, noise);
-	UV = rnd_gauss (uv, 1.4*noise);
+	V = gaussdev_rnd (v, noise);
+	UV = gaussdev_rnd (uv, 1.4*noise);
 	if (noise < dMo) {
 	  fprintf (stdout, "%f %f   %f %f   %f  %f\n", V, noise, V+UV, noise, mass, age);
@@ -96,55 +85,4 @@
 }
 
-
-double
-rnd_gauss (mean, sigma)
-double mean, sigma; 
-{
-
-  int i;
-  double y;
-
-  y = drand48();
-  i = NGAUSS*y;
-  y = gaussint[i]*sigma + mean;
-
-  return (y);
-
-}
-
-
-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++;
-    }
-  }
-  return (x + dx / 2.0);
-}
 
 double 
@@ -172,18 +110,4 @@
   return (x + dx / 2.0);
 }
-
-double 
-gaussian (x, mean, sigma) 
-double x, mean, sigma; 
-{
-
-  double f, X;
-
-  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
-
-  return (f);
-
-}
-
 
 /*****************************************************************************/
Index: trunk/Ohana/src/misc/src/fakestars.c
===================================================================
--- trunk/Ohana/src/misc/src/fakestars.c	(revision 39241)
+++ trunk/Ohana/src/misc/src/fakestars.c	(revision 39242)
@@ -17,8 +17,4 @@
   int found;
 } Stars;
-
-void   gauss_init ();
-double gaussian (double x, double mean, double sigma);
-double rnd_gauss (double mean, double sigma);
 
   /* 
@@ -45,5 +41,5 @@
   offset = atof (argv[4]);
 
-  gauss_init ();
+  gaussdev_init ();
 
   /* load header, open file */
@@ -98,5 +94,5 @@
     dMr = 0.001 * stars[i].dM;
     dMo = hypot (dMs, dMr);
-    dM  = rnd_gauss (0.0, dMo);
+    dM  = gaussdev_rnd (0.0, dMo);
     stars[i].M += dM + offset;
   }
Index: trunk/Ohana/src/misc/src/fitdist.c
===================================================================
--- trunk/Ohana/src/misc/src/fitdist.c	(revision 39241)
+++ trunk/Ohana/src/misc/src/fitdist.c	(revision 39242)
@@ -1,9 +1,3 @@
 # include <ohana.h>
-
-# define NGAUSS 2048
-double drand48();
-double gaussint[NGAUSS];
-double gaussian();
-double rnd_gauss();
 
 # define MMIN 1.0
@@ -45,5 +39,5 @@
   fprintf (stderr, "read in parameters\n");
 
-  gauss_init ();
+  gaussdev_init ();
 
   read_datafiles ();
@@ -360,67 +354,4 @@
 /*****************************************************************************/
 
-double
-rnd_gauss (mean, sigma)
-double mean, sigma; 
-{
-
-  int i;
-  double y;
-
-  y = drand48();
-  i = NGAUSS*y;
-  y = gaussint[i]*sigma + mean;
-
-  return (y);
-
-}
-
-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++;
-    }
-  }
-  return (x + dx / 2.0);
-}
-
-double 
-gaussian (x, mean, sigma) 
-double x, mean, sigma; 
-{
-
-  double f, X;
-
-  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
-
-  return (f);
-
-}
-
 /*****************************************************************************/
 
@@ -453,6 +384,6 @@
       Ngood = 1.0;
       for (k = 0; k < NTRY; k++) {
-	v = rnd_gauss (Mag, Noise);
-	uv = rnd_gauss (Color, 1.4*Noise);
+	v = gaussdev_rnd (Mag, Noise);
+	uv = gaussdev_rnd (Color, 1.4*Noise);
 	x = (uv - UV0 - 0.7*A_V) / DUV;
 	y = (v - V0 - Dist - A_V) / DV;
Index: trunk/Ohana/src/misc/src/magtoage.c
===================================================================
--- trunk/Ohana/src/misc/src/magtoage.c	(revision 39241)
+++ trunk/Ohana/src/misc/src/magtoage.c	(revision 39242)
@@ -2,12 +2,4 @@
 # define MMIN 1.0
 # define MMAX 120.0
-extern double drand48();
-extern double rnd_gauss();
-extern double rnd_integrate ();
-extern double gaussian ();
-
-double gaussian();
-double rnd_gauss();
-double rnd_integrate();
 
 void main (argc, argv)
@@ -29,4 +21,6 @@
   long A, B;
   
+  gaussdev_init();
+
   lAo =   0.0;
   ldA =   1.0;
@@ -118,6 +112,6 @@
       }
       else {
-	v = rnd_gauss (V, dV);
-	uv = rnd_gauss ((U-V), dUV);
+	v = gaussdev_rnd (V, dV);
+	uv = gaussdev_rnd ((U-V), dUV);
       }
       x = (uv - UV0 - 0.7*Av) / DUV;
@@ -167,55 +161,2 @@
 }
 
-double
-rnd_gauss (mean, sigma)
-double mean, sigma; 
-{
-
-  double range, x;
-
-  range = drand48();
-  x = rnd_integrate (*gaussian, range, mean, sigma);
-
-  return (x);
-
-}
-
-
-double 
-rnd_integrate (function, range, mean, sigma) 
-double (*function) ();
-double range, mean, sigma;
-{
-
-  double val, x, dx, dx1, dx2, dx3, df;
-
-  range += 0.0001;
-  val = 0;
-  dx = sigma / 10.0;
-  dx1 = dx / 3.0;
-  dx2 = 2.0*dx/3.0;
-  dx3 = dx;
-
-  for (x = mean - 7*sigma; (val < range) && (x < mean + 7*sigma); x += dx)  {
-    df = (3.0*function(x    , mean, sigma) + 
-	  9.0*function(x+dx1, mean, sigma) +
-	  9.0*function(x+dx2, mean, sigma) + 
-	  3.0*function(x+dx3, mean, sigma)) * (dx1/8.0);
-    val += df;
-  }
-  return (x + dx / 2.0);
-}
-
-double 
-gaussian (x, mean, sigma) 
-double x, mean, sigma; 
-{
-
-  double f, X;
-
-  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
-
-  return (f);
-
-}
-
Index: trunk/Ohana/src/misc/src/magtomass.c
===================================================================
--- trunk/Ohana/src/misc/src/magtomass.c	(revision 39241)
+++ trunk/Ohana/src/misc/src/magtomass.c	(revision 39242)
@@ -1,15 +1,4 @@
 # include <ohana.h>
-extern double drand48();
-extern double rnd_gauss();
-extern double rnd_integrate ();
-extern double gaussian ();
 # define NEWWAY 1
-
-double gaussian();
-double rnd_gauss();
-double rnd_integrate();
-
-# define NGAUSS 2048
-double gaussint[NGAUSS];
 
 void main (argc, argv)
@@ -52,5 +41,5 @@
   fprintf (stderr, "using mags in columns %d & %d\n", col1, col2);
 
-  gauss_init ();
+  gaussdev_init ();
 
   lAo = 0.1;
@@ -147,6 +136,6 @@
 	Ngood = 1.0;
 	for (k = 0; k < Ntry; k++) {
-	  v = rnd_gauss (V, dV);
-	  uv = rnd_gauss ((U-V), dUV);
+	  v = gaussdev_rnd (V, dV);
+	  uv = gaussdev_rnd ((U-V), dUV);
 	  x = (uv - UV0 - 0.7*Av) / DUV;
 	  y = (v - d - V0 - Av) / DV;
@@ -202,91 +191,2 @@
 */
 
-double
-rnd_gauss (mean, sigma)
-double mean, sigma; 
-{
-
-  int i;
-  double y;
-
-  y = drand48();
-  i = NGAUSS*y;
-  y = gaussint[i]*sigma + mean;
-
-  return (y);
-
-}
-
-
-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++;
-    }
-  }
-  return (x + dx / 2.0);
-}
-
-double 
-rnd_integrate (function, range, mean, sigma) 
-double (*function) ();
-double range, mean, sigma;
-{
-
-  double val, x, dx, dx1, dx2, dx3, df;
-
-  range += 0.0001;
-  val = 0;
-  dx = sigma / 100.0;
-  dx1 = dx / 3.0;
-  dx2 = 2.0*dx/3.0;
-  dx3 = dx;
-
-  for (x = mean - 7*sigma; (val < range) && (x < mean + 7*sigma); x += dx)  {
-    df = (3.0*function(x    , mean, sigma) + 
-	  9.0*function(x+dx1, mean, sigma) +
-	  9.0*function(x+dx2, mean, sigma) + 
-	  3.0*function(x+dx3, mean, sigma)) * (dx1/8.0);
-    val += df;
-  }
-  return (x + dx / 2.0);
-}
-
-double 
-gaussian (x, mean, sigma) 
-double x, mean, sigma; 
-{
-
-  double f, X;
-
-  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
-
-  return (f);
-
-}
-
Index: trunk/Ohana/src/misc/src/mkgauss.c
===================================================================
--- trunk/Ohana/src/misc/src/mkgauss.c	(revision 39241)
+++ trunk/Ohana/src/misc/src/mkgauss.c	(revision 39242)
@@ -1,11 +1,3 @@
 # 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) {
@@ -23,8 +15,8 @@
   Npts = atoi (argv[3]);
 
-  gauss_init ();
+  gaussdev_init ();
 
   for (i = 0; i < Npts; i++) {
-    f = rnd_gauss (mean, sigma);
+    f = gaussdev_rnd (mean, sigma);
     fprintf (stdout, "%f\n", f);
   }
@@ -33,53 +25,2 @@
 }
 
-void gauss_init () {
-
-  int i;
-  long A, B;
-  double val, x, dx, dx1, dx2, dx3, df;
-
-  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;
-
-  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);
-}
Index: trunk/Ohana/src/opihi/cmd.data/gaussdeviate.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/gaussdeviate.c	(revision 39241)
+++ trunk/Ohana/src/opihi/cmd.data/gaussdeviate.c	(revision 39242)
@@ -17,7 +17,7 @@
   ResetVector (vec, OPIHI_FLT, Npts);
 
-  gauss_init (2048);
+  gaussdev_init ();
   for (i = 0; i < Npts; i++) {
-    vec[0].elements.Flt[i] = rnd_gauss (mean, sigma);
+    vec[0].elements.Flt[i] = gaussdev_rnd (mean, sigma);
   }
   return (TRUE);
@@ -28,6 +28,4 @@
     
 }
-
-double int_gauss (int i);
 
 int gaussintegral (int argc, char **argv) {
@@ -46,7 +44,7 @@
   ResetVector (vec, OPIHI_FLT, Npts);
 
-  gauss_init (Npts);
+  gaussdev_init ();
   for (i = 0; i < Npts; i++) {
-    vec[0].elements.Flt[i] = int_gauss (i);
+    vec[0].elements.Flt[i] = gaussian_int (i);
   }
   return (TRUE);
Index: trunk/Ohana/src/opihi/lib.data/Makefile
===================================================================
--- trunk/Ohana/src/opihi/lib.data/Makefile	(revision 39241)
+++ trunk/Ohana/src/opihi/lib.data/Makefile	(revision 39242)
@@ -30,5 +30,4 @@
 $(SDIR)/starfuncs.$(ARCH).o		\
 $(SDIR)/hermitian.$(ARCH).o		\
-$(SDIR)/gaussian.$(ARCH).o		\
 $(SDIR)/graphtools.$(ARCH).o            \
 $(SDIR)/queues.$(ARCH).o		\
@@ -36,4 +35,7 @@
 $(SDIR)/open_kapa.$(ARCH).o             \
 $(SDIR)/style_args.$(ARCH).o
+
+# moved to libohana
+# $(SDIR)/gaussian.$(ARCH).o		\
 
 # fix malloc
Index: trunk/Ohana/src/opihi/mana/simsignal.c
===================================================================
--- trunk/Ohana/src/opihi/mana/simsignal.c	(revision 39241)
+++ trunk/Ohana/src/opihi/mana/simsignal.c	(revision 39242)
@@ -6,5 +6,4 @@
   float *buf;
   double cvalue, dvalue, sigma, SN, period;
-  double rnd_gauss ();
 
   if (argc != 5) {
@@ -28,5 +27,5 @@
     ivalue = scale * 0.5 * (sin (i*2*M_PI/period) + 1) + 0.5;
     /*
-    dvalue = rnd_gauss (cvalue, sigma);
+    dvalue = gaussdev_rnd (cvalue, sigma);
     cvalue = (dvalue + range) / (2.0*range);
     dvalue = MAX (0, MIN (0.99999, cvalue));
@@ -44,5 +43,5 @@
 8 bit = 2^8
 
-  gauss_init (2*scale);
+  gaussdev_init ();
   sigma = 2.0 / SN;
   range = 1 + 5*sigma;