Changeset 39242 for trunk

Timestamp:

Dec 9, 2015, 3:42:26 PM (11 years ago)

Author:

eugene

Message:

move multple copies of gaussian deviates to libohana

Location:

trunk/Ohana/src

Files:

• addstar/src/mkcmf.c (modified) (13 diffs)
• fakeastro/Makefile (modified) (1 diff)
• fakeastro/src/make_2mass_measures.c (modified) (1 diff)
• fakeastro/src/make_fake_stars_catalog.c (modified) (1 diff)
• fakeastro/src/make_fakeqsos.c (modified) (1 diff)
• fakeastro/src/make_fakestars.c (modified) (2 diffs)
• libohana/Makefile (modified) (1 diff)
• libohana/include/ohana.h (modified) (1 diff)
• libohana/src/gaussdev.c (added)
• misc/src/fakedist.c (modified) (5 diffs)
• misc/src/fakepop.c (modified) (5 diffs)
• misc/src/fakestars.c (modified) (3 diffs)
• misc/src/fitdist.c (modified) (4 diffs)
• misc/src/magtoage.c (modified) (4 diffs)
• misc/src/magtomass.c (modified) (4 diffs)
• misc/src/mkgauss.c (modified) (3 diffs)
• opihi/cmd.data/gaussdeviate.c (modified) (3 diffs)
• opihi/lib.data/Makefile (modified) (2 diffs)
• opihi/mana/simsignal.c (modified) (3 diffs)

trunk/Ohana/src/addstar/src/mkcmf.c

@@ -13 +13 @@
 # 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 +280 @@
   }
     
-  gauss_init (2048);
+  gaussdev_init ();
 
   // load test stars from a file:
@@ -564 +562 @@
 
 
-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 +578 @@
 
     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 +624 @@
 
     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 +673 @@
 
     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 +724 @@
 
     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 +781 @@
 
     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 +842 @@
 
     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 +911 @@
 
     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 +1018 @@
 
     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 +1122 @@
 
     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 +1248 @@
 
     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);

trunk/Ohana/src/fakeastro/Makefile

@@ -38 +38 @@
 $(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 \

trunk/Ohana/src/fakeastro/src/make_2mass_measures.c

@@ -80 +80 @@
 
     // 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;

trunk/Ohana/src/fakeastro/src/make_fake_stars_catalog.c

@@ -113 +113 @@
 
     // 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;

trunk/Ohana/src/fakeastro/src/make_fakeqsos.c

@@ -51 +51 @@
 
     // 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;

trunk/Ohana/src/fakeastro/src/make_fakestars.c

@@ -48 +48 @@
       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 +98 @@
     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);
     }
     

trunk/Ohana/src/libohana/Makefile

@@ -35 +35 @@
 $(SRC)/sorts.$(ARCH).o           \
 $(SRC)/bisection.$(ARCH).o       \
+$(SRC)/gaussdev.$(ARCH).o       \
 $(SRC)/string.$(ARCH).o          \
 $(SRC)/findexec.$(ARCH).o        \

trunk/Ohana/src/libohana/include/ohana.h

@@ -486 +486 @@
 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

trunk/Ohana/src/misc/src/fakedist.c

@@ -1 +1 @@
 # include <ohana.h>
-
-# define NGAUSS 2048
-double drand48();
-double gaussint[NGAUSS];
-double gaussian();
-double rnd_gauss();
 
 # define MMIN 1.0
@@ -45 +39 @@
   init_outmatrix (); 
  
-  gauss_init ();
+  gaussdev_init ();
 
   read_datafiles ();
@@ -349 +343 @@
         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 +408 @@
 
 /*****************************************************************************/
-
-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 +436 @@
       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;

trunk/Ohana/src/misc/src/fakepop.c

@@ -2 +2 @@
 # 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 +36 @@
   gfits_scan (&UV_h, "DEC_Y", "%lf", 1, &ldM);
   
-  gauss_init ();
+  gaussdev_init ();
 
   fprintf (stderr, "beginning main loop\n");
@@ -68 +57 @@
         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 +85 @@
 }
 
-
-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 +110 @@
   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);
-
-}
-
 
 /*****************************************************************************/

trunk/Ohana/src/misc/src/fakestars.c

@@ -17 +17 @@
   int found;
 } Stars;
-
-void   gauss_init ();
-double gaussian (double x, double mean, double sigma);
-double rnd_gauss (double mean, double sigma);
 
   /* 
@@ -45 +41 @@
   offset = atof (argv[4]);
 
-  gauss_init ();
+  gaussdev_init ();
 
   /* load header, open file */
@@ -98 +94 @@
     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;
   }

trunk/Ohana/src/misc/src/fitdist.c

@@ -1 +1 @@
 # include <ohana.h>
-
-# define NGAUSS 2048
-double drand48();
-double gaussint[NGAUSS];
-double gaussian();
-double rnd_gauss();
 
 # define MMIN 1.0
@@ -45 +39 @@
   fprintf (stderr, "read in parameters\n");
 
-  gauss_init ();
+  gaussdev_init ();
 
   read_datafiles ();
@@ -360 +354 @@
 /*****************************************************************************/
 
-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 +384 @@
       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;

trunk/Ohana/src/misc/src/magtoage.c

@@ -2 +2 @@
 # 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 +21 @@
   long A, B;
   
+  gaussdev_init();
+
   lAo =   0.0;
   ldA =   1.0;
@@ -118 +112 @@
       }
       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 +161 @@
 }
 
-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);
-
-}
-

trunk/Ohana/src/misc/src/magtomass.c

@@ -1 +1 @@
 # 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 +41 @@
   fprintf (stderr, "using mags in columns %d & %d\n", col1, col2);
 
-  gauss_init ();
+  gaussdev_init ();
 
   lAo = 0.1;
@@ -147 +136 @@
         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 +191 @@
 */
 
-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);
-
-}
-

trunk/Ohana/src/misc/src/mkgauss.c

@@ -1 +1 @@
 # 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 +15 @@
   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 +25 @@
 }
 
-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);
-}

trunk/Ohana/src/opihi/cmd.data/gaussdeviate.c

@@ -17 +17 @@
   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 +28 @@
     
 }
-
-double int_gauss (int i);
 
 int gaussintegral (int argc, char **argv) {
@@ -46 +44 @@
   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);

trunk/Ohana/src/opihi/lib.data/Makefile

@@ -30 +30 @@
 $(SDIR)/starfuncs.$(ARCH).o             \
 $(SDIR)/hermitian.$(ARCH).o             \
-$(SDIR)/gaussian.$(ARCH).o              \
 $(SDIR)/graphtools.$(ARCH).o            \
 $(SDIR)/queues.$(ARCH).o                \
@@ -36 +35 @@
 $(SDIR)/open_kapa.$(ARCH).o             \
 $(SDIR)/style_args.$(ARCH).o
+
+# moved to libohana
+# $(SDIR)/gaussian.$(ARCH).o            \
 
 # fix malloc

trunk/Ohana/src/opihi/mana/simsignal.c

@@ -6 +6 @@
   float *buf;
   double cvalue, dvalue, sigma, SN, period;
-  double rnd_gauss ();
 
   if (argc != 5) {
@@ -28 +27 @@
     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 +43 @@
 8 bit = 2^8
 
-  gauss_init (2*scale);
+  gaussdev_init ();
   sigma = 2.0 / SN;
   range = 1 + 5*sigma;