Move gaussian_noise() to utils.c and synergise it with poisson_noise()

3 files changed, 12 insertions, 26 deletions

diff --git a/src/partial_sim.c b/src/partial_sim.c
index 55e8521d..c0c44212 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -31,20 +31,6 @@
 #include "stream.h"
 
 
-static double gaussian_noise(double expected, double variance)
-{
-	double x1, x2, noise;
-
-	/* A uniformly distributed random number between 0 and 1 */
-	x1 = ((double)random()/RAND_MAX);
-	x2 = ((double)random()/RAND_MAX);
-
-	noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2);
-
-	return expected + noise*variance;
-}
-
-
 static void mess_up_cell(UnitCell *cell)
 {
 	double ax, ay, az;
diff --git a/src/utils.c b/src/utils.c
index 7e61f250..bfa5c7d4 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -120,24 +120,23 @@ void progress_bar(int val, int total, const char *text)
 }
 
 
-static int fake_poisson_noise(double expected)
+double gaussian_noise(double expected, double stddev)
 {
-	double x1, x2, w;
-	double noise, rf;
-
-	do {
+	double x1, x2, noise;
 
-		x1 = 2.0 * ((double)random()/RAND_MAX) - 1.0;
-		x2 = 2.0 * ((double)random()/RAND_MAX) - 1.0;
-		w = pow(x1, 2.0) + pow(x2, 2.0);
+	/* A uniformly distributed random number between 0 and 1 */
+	x1 = ((double)random()/RAND_MAX);
+	x2 = ((double)random()/RAND_MAX);
 
-	} while ( w >= 1.0 );
+	noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2);
 
-	w = sqrt((-2.0*log(w))/w);
-	noise = w * x1;
+	return expected + noise*stddev;
+}
 
-	rf = expected + noise*sqrt(expected);
 
+static int fake_poisson_noise(double expected)
+{
+	double rf = gaussian_noise(expected, sqrt(expected));
 	return (int)rf;
 }
 
diff --git a/src/utils.h b/src/utils.h
index e3c78d6c..4433023f 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -91,6 +91,7 @@ extern int assplode(const char *a, const char *delims, char ***pbits,
                     AssplodeFlag flags);
 
 extern void progress_bar(int val, int total, const char *text);
+extern double gaussian_noise(double expected, double stddev);
 extern int poisson_noise(double expected);
 
 /* Keep these ones inline, to avoid function call overhead */