Index: trunk/Ohana/src/addstar/src/mkcmf.c
===================================================================
--- trunk/Ohana/src/addstar/src/mkcmf.c	(revision 39225)
+++ 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);