Changeset 39242 for trunk/Ohana/src/misc/src/magtomass.c

Timestamp:

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

Author:

eugene

Message:

move multple copies of gaussian deviates to libohana

File:

• trunk/Ohana/src/misc/src/magtomass.c (modified) (4 diffs)

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);
-
-}
-