RNG overhaul

Previously, we were using random(), which is really really bad.

10 files changed, 151 insertions, 82 deletions

diff --git a/tests/cell_check.c b/tests/cell_check.c
index f62cd7be..009b4c1f 100644
--- a/tests/cell_check.c
+++ b/tests/cell_check.c
@@ -3,7 +3,11 @@
  *
  * Check that unit cells work correctly
  *
- * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2012-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -49,6 +53,9 @@ int main(int argc, char *argv[])
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	cell = cell_new_from_parameters(27.155e-9, 28.155e-9, 10.987e-9,
 	                                deg2rad(90.0),
@@ -56,7 +63,7 @@ int main(int argc, char *argv[])
 	                                deg2rad(120.0));
 	if ( cell == NULL ) return 1;
 
-	orientation = random_quaternion();
+	orientation = random_quaternion(rng);
 	cell = cell_rotate(cell, orientation);
 
 	cell_get_reciprocal(cell, &asx, &asy, &asz,
@@ -112,6 +119,7 @@ int main(int argc, char *argv[])
 	                         cx, cy, cz);
 	cell_print(cell);
 
+	gsl_rng_free(rng);
 
 	return fail;
 }
diff --git a/tests/centering_check.c b/tests/centering_check.c
index f4073072..887311dd 100644
--- a/tests/centering_check.c
+++ b/tests/centering_check.c
@@ -3,7 +3,11 @@
  *
  * Check that centering of cells works
  *
- * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2012-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -53,7 +57,7 @@ static int check_cell(UnitCell *cell, const char *text)
 
 static int check_centering(double a, double b, double c,
                            double al, double be, double ga,
-                           LatticeType latt, char cen, char ua)
+                           LatticeType latt, char cen, char ua, gsl_rng *rng)
 {
 	UnitCell *cell, *cref;
 	UnitCell *n;
@@ -77,7 +81,7 @@ static int check_centering(double a, double b, double c,
 	cell_set_centering(cref, cen);
 	cell_set_unique_axis(cref, ua);
 
-	cell = cell_rotate(cref, random_quaternion());
+	cell = cell_rotate(cref, random_quaternion(rng));
 	if ( cell == NULL ) return 1;
 	cell_free(cref);
 
@@ -110,9 +114,9 @@ static int check_centering(double a, double b, double c,
 
 		do {
 
-			h = flat_noise(0, 30);
-			k = flat_noise(0, 30);
-			l = flat_noise(0, 30);
+			h = gsl_rng_uniform_int(rng, 30);
+			k = gsl_rng_uniform_int(rng, 30);
+			l = gsl_rng_uniform_int(rng, 30);
 
 		} while ( forbidden_reflection(cell, h, k, l) );
 
@@ -156,9 +160,9 @@ static int check_centering(double a, double b, double c,
 		int f = 0;
 		long int ih, ik, il;
 
-		h = flat_noise(0, 30);
-		k = flat_noise(0, 30);
-		l = flat_noise(0, 30);
+		h = gsl_rng_uniform_int(rng, 30);
+		k = gsl_rng_uniform_int(rng, 30);
+		l = gsl_rng_uniform_int(rng, 30);
 
 		x = h*asx + k*bsx + l*csx;
 		y = h*asy + k*bsy + l*csy;
@@ -199,110 +203,115 @@ static int check_centering(double a, double b, double c,
 int main(int argc, char *argv[])
 {
 	int fail = 0;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	/* Triclinic P */
 	fail += check_centering(50e-10, 55e-10, 70e-10, 67.0, 70.0, 77.0,
-	                        L_TRICLINIC, 'P', '*');
+	                        L_TRICLINIC, 'P', '*', rng);
 
 	/* Monoclinic P */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
-	                        L_MONOCLINIC, 'P', 'a');
+	                        L_MONOCLINIC, 'P', 'a', rng);
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
-	                        L_MONOCLINIC, 'P', 'b');
+	                        L_MONOCLINIC, 'P', 'b', rng);
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
-	                        L_MONOCLINIC, 'P', 'c');
+	                        L_MONOCLINIC, 'P', 'c', rng);
 
 	/* Monoclinic "C"-centered, unique axis a, three cell choices */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
-	                        L_MONOCLINIC, 'B', 'a');
+	                        L_MONOCLINIC, 'B', 'a', rng);
 	fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
-	                        L_MONOCLINIC, 'C', 'a');
+	                        L_MONOCLINIC, 'C', 'a', rng);
 	fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
-	                        L_MONOCLINIC, 'I', 'a');
+	                        L_MONOCLINIC, 'I', 'a', rng);
 
 	/* Monoclinic "C"-centered, unique axis b, three cell choices */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
-	                        L_MONOCLINIC, 'C', 'b');
+	                        L_MONOCLINIC, 'C', 'b', rng);
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
-	                        L_MONOCLINIC, 'A', 'b');
+	                        L_MONOCLINIC, 'A', 'b', rng);
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
-	                        L_MONOCLINIC, 'I', 'b');
+	                        L_MONOCLINIC, 'I', 'b', rng);
 
 	/* Monoclinic "C"-centered, unique axis c, three cell choices */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
-	                        L_MONOCLINIC, 'A', 'c');
+	                        L_MONOCLINIC, 'A', 'c', rng);
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
-	                        L_MONOCLINIC, 'B', 'c');
+	                        L_MONOCLINIC, 'B', 'c', rng);
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
-	                        L_MONOCLINIC, 'I', 'c');
+	                        L_MONOCLINIC, 'I', 'c', rng);
 
 	/* Orthorhombic P */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'P', '*');
+	                        L_ORTHORHOMBIC, 'P', '*', rng);
 
 	/* Orthorhombic A */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'A', '*');
+	                        L_ORTHORHOMBIC, 'A', '*', rng);
 
 	/* Orthorhombic B */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'B', '*');
+	                        L_ORTHORHOMBIC, 'B', '*', rng);
 
 	/* Orthorhombic C */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'C', '*');
+	                        L_ORTHORHOMBIC, 'C', '*', rng);
 
 	/* Orthorhombic I */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'I', '*');
+	                        L_ORTHORHOMBIC, 'I', '*', rng);
 
 	/* Orthorhombic F */
 	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_ORTHORHOMBIC, 'F', '*');
+	                        L_ORTHORHOMBIC, 'F', '*', rng);
 
 	/* Tetragonal P */
 	fail += check_centering(10e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_TETRAGONAL, 'P', 'a');
+	                        L_TETRAGONAL, 'P', 'a', rng);
 	fail += check_centering(30e-10, 10e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_TETRAGONAL, 'P', 'b');
+	                        L_TETRAGONAL, 'P', 'b', rng);
 	fail += check_centering(30e-10, 30e-10, 10e-10, 90.0, 90.0, 90.0,
-	                        L_TETRAGONAL, 'P', 'c');
+	                        L_TETRAGONAL, 'P', 'c', rng);
 
 	/* Tetragonal I */
 	fail += check_centering(10e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_TETRAGONAL, 'I', 'a');
+	                        L_TETRAGONAL, 'I', 'a', rng);
 	fail += check_centering(30e-10, 10e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_TETRAGONAL, 'I', 'b');
+	                        L_TETRAGONAL, 'I', 'b', rng);
 	fail += check_centering(30e-10, 30e-10, 10e-10, 90.0, 90.0, 90.0,
-	                        L_TETRAGONAL, 'I', 'c');
+	                        L_TETRAGONAL, 'I', 'c', rng);
 
 	/* Rhombohedral R */
 	fail += check_centering(10e-10, 10e-10, 10e-10, 60.0, 60.0, 60.0,
-	                        L_RHOMBOHEDRAL, 'R', '*');
+	                        L_RHOMBOHEDRAL, 'R', '*', rng);
 
 	/* Hexagonal P */
 	fail += check_centering(30e-10, 10e-10, 10e-10, 120.0, 90.0, 90.0,
-	                        L_HEXAGONAL, 'P', 'a');
+	                        L_HEXAGONAL, 'P', 'a', rng);
 	fail += check_centering(10e-10, 30e-10, 10e-10, 90.0, 120.0, 90.0,
-	                        L_HEXAGONAL, 'P', 'b');
+	                        L_HEXAGONAL, 'P', 'b', rng);
 	fail += check_centering(10e-10, 10e-10, 30e-10, 90.0, 90.0, 120.0,
-	                        L_HEXAGONAL, 'P', 'c');
+	                        L_HEXAGONAL, 'P', 'c', rng);
 
 	/* Hexagonal H (PDB-speak for rhombohedral) */
 	fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0,
-	                        L_HEXAGONAL, 'H', 'c');
+	                        L_HEXAGONAL, 'H', 'c', rng);
 
 	/* Cubic P */
 	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_CUBIC, 'P', '*');
+	                        L_CUBIC, 'P', '*', rng);
 
 	/* Cubic I */
 	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_CUBIC, 'I', '*');
+	                        L_CUBIC, 'I', '*', rng);
 
 	/* Cubic F */
 	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
-	                        L_CUBIC, 'F', '*');
+	                        L_CUBIC, 'F', '*', rng);
+
+	gsl_rng_free(rng);
 
 	return fail;
 }
diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c
index c183a2a2..d2475213 100644
--- a/tests/gpu_sim_check.c
+++ b/tests/gpu_sim_check.c
@@ -3,7 +3,11 @@
  *
  * Check that GPU simulation agrees with CPU version
  *
- * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2012-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -80,6 +84,9 @@ int main(int argc, char *argv[])
 	double start, end;
 	double gpu_time, cpu_time;
 	SymOpList *sym;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	gctx = setup_gpu(1, NULL, NULL, NULL, 0);
 	if ( gctx == NULL ) {
@@ -90,7 +97,7 @@ int main(int argc, char *argv[])
 	cell_raw = cell_new_from_parameters(28.1e-9, 28.1e-9, 16.5e-9,
 	                          deg2rad(90.0), deg2rad(90.0), deg2rad(120.0));
 
-	cell = cell_rotate(cell_raw, random_quaternion());
+	cell = cell_rotate(cell_raw, random_quaternion(rng));
 
 	gpu_image.width = 1024;
 	gpu_image.height = 1024;
@@ -217,5 +224,7 @@ int main(int argc, char *argv[])
 
 	}
 
+	gsl_rng_free(rng);
+
 	return 0;
 }
diff --git a/tests/integration_check.c b/tests/integration_check.c
index f1b5384b..0196742f 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -43,7 +43,7 @@ int main(int argc, char *argv[])
 	struct image image;
 	int fs, ss;
 	FILE *fh;
-	unsigned int seed;
+	unsigned long int seed;
 	int fail = 0;
 	const int w = 128;
 	const int h = 128;
@@ -57,11 +57,14 @@ int main(int argc, char *argv[])
 	int i;
 	Histogram *hi;
 	double esd_sum = 0.0;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	fh = fopen("/dev/urandom", "r");
 	fread(&seed, sizeof(seed), 1, fh);
 	fclose(fh);
-	srand(seed);
+	gsl_rng_set(rng, seed);
 
 	image.flags = NULL;
 	image.beam = NULL;
@@ -110,7 +113,7 @@ int main(int argc, char *argv[])
 
 		for ( fs=0; fs<w; fs++ ) {
 		for ( ss=0; ss<h; ss++ ) {
-			image.dp[0][fs+w*ss] = 10.0*poisson_noise(40);
+			image.dp[0][fs+w*ss] = 10.0*poisson_noise(rng, 40);
 			if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 2 ) continue;
 			//image.dp[0][fs+w*ss] += 10.0*poisson_noise(10);
 		}
@@ -124,7 +127,7 @@ int main(int argc, char *argv[])
 		cell_set_centering(cell, 'P');
 		cell_set_parameters(cell, 800.0e-10, 800.0e-10, 800.0e-10,
 		                    deg2rad(90.0), deg2rad(90.0), deg2rad(90.0));
-		cell = cell_rotate(cell, random_quaternion());
+		cell = cell_rotate(cell, random_quaternion(rng));
 
 		ic.halfw = ir_out;
 		ic.image = &image;
diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c
index a4fd4526..4b3894ba 100644
--- a/tests/pr_l_gradient_check.c
+++ b/tests/pr_l_gradient_check.c
@@ -3,11 +3,11 @@
  *
  * Check Lorentz factor gradients for post refinement
  *
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2012-2013 Thomas White <taw@physics.org>
+ *   2012-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -284,6 +284,7 @@ int main(int argc, char *argv[])
 	int quiet = 0;
 	int plot = 0;
 	int c;
+	gsl_rng *rng;
 
 	const struct option longopts[] = {
 		{"quiet",       0, &quiet,        1},
@@ -334,12 +335,14 @@ int main(int argc, char *argv[])
 	                                deg2rad(90.0),
 	                                deg2rad(90.0));
 
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+
 	for ( i=0; i<1; i++ ) {
 
 		UnitCell *rot;
 		double val;
 
-		orientation = random_quaternion();
+		orientation = random_quaternion(rng);
 		rot = cell_rotate(cell, orientation);
 		crystal_set_cell(cr, rot);
 
@@ -354,5 +357,7 @@ int main(int argc, char *argv[])
 
 	}
 
+	gsl_rng_free(rng);
+
 	return fail;
 }
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index 9adf4780..c0582036 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -3,11 +3,11 @@
  *
  * Check partiality gradients for post refinement
  *
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2012-2013 Thomas White <taw@physics.org>
+ *   2012-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -397,6 +397,7 @@ int main(int argc, char *argv[])
 	int quiet = 0;
 	int plot = 0;
 	int c;
+	gsl_rng *rng;
 
 	const struct option longopts[] = {
 		{"quiet",       0, &quiet,        1},
@@ -447,12 +448,14 @@ int main(int argc, char *argv[])
 	                                deg2rad(90.0),
 	                                deg2rad(90.0));
 
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+
 	for ( i=0; i<1; i++ ) {
 
 		UnitCell *rot;
 		double val;
 
-		orientation = random_quaternion();
+		orientation = random_quaternion(rng);
 		rot = cell_rotate(cell, orientation);
 		crystal_set_cell(cr, rot);
 
@@ -511,5 +514,7 @@ int main(int argc, char *argv[])
 
 	}
 
+	gsl_rng_free(rng);
+
 	return fail;
 }
diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c
index f603fff0..5ef01bf5 100644
--- a/tests/pr_pl_gradient_check.c
+++ b/tests/pr_pl_gradient_check.c
@@ -3,11 +3,11 @@
  *
  * Check partiality gradients for post refinement
  *
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2012-2013 Thomas White <taw@physics.org>
+ *   2012-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -399,6 +399,7 @@ int main(int argc, char *argv[])
 	int quiet = 0;
 	int plot = 0;
 	int c;
+	gsl_rng *rng;
 
 	const struct option longopts[] = {
 		{"quiet",       0, &quiet,        1},
@@ -449,12 +450,14 @@ int main(int argc, char *argv[])
 	                                deg2rad(90.0),
 	                                deg2rad(90.0));
 
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+
 	for ( i=0; i<1; i++ ) {
 
 		UnitCell *rot;
 		double val;
 
-		orientation = random_quaternion();
+		orientation = random_quaternion(rng);
 		rot = cell_rotate(cell, orientation);
 		crystal_set_cell(cr, rot);
 
@@ -514,5 +517,7 @@ int main(int argc, char *argv[])
 
 	}
 
+	gsl_rng_free(rng);
+
 	return fail;
 }
diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c
index c1718a92..f1665a6c 100644
--- a/tests/prof2d_check.c
+++ b/tests/prof2d_check.c
@@ -3,7 +3,11 @@
  *
  * Check 2D profile fitting
  *
- * Copyright © 2013 Thomas White <taw@physics.org>
+ * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2013-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -48,7 +52,7 @@ int main(int argc, char *argv[])
 	struct image image;
 	int fs, ss;
 	FILE *fh;
-	unsigned int seed;
+	unsigned long int seed;
 	int fail = 0;
 	const int w = 1024;
 	const int h = 1024;
@@ -65,11 +69,14 @@ int main(int argc, char *argv[])
 	int n = 0;
 	int n_strong = 0;
 	int n_weak = 0;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	fh = fopen("/dev/urandom", "r");
 	fread(&seed, sizeof(seed), 1, fh);
 	fclose(fh);
-	srand(seed);
+	gsl_rng_set(rng, seed);
 
 	image.flags = NULL;
 	image.beam = NULL;
@@ -121,7 +128,7 @@ int main(int argc, char *argv[])
 	cell_set_centering(cell, 'P');
 	cell_set_parameters(cell, 800.0e-10, 800.0e-10, 800.0e-10,
 	                    deg2rad(90.0), deg2rad(90.0), deg2rad(90.0));
-	cell = cell_rotate(cell, random_quaternion());
+	cell = cell_rotate(cell, random_quaternion(rng));
 
 	cr = crystal_new();
 	crystal_set_profile_radius(cr, 0.001e9);
@@ -136,7 +143,7 @@ int main(int argc, char *argv[])
 
 	for ( fs=0; fs<w; fs++ ) {
 	for ( ss=0; ss<h; ss++ ) {
-		image.dp[0][fs+w*ss] = 10.0*poisson_noise(40);
+		image.dp[0][fs+w*ss] = 10.0*poisson_noise(rng, 40);
 	}
 	}
 
@@ -154,11 +161,11 @@ int main(int argc, char *argv[])
 		get_indices(refl, &hj, &kj, &lj);
 		if ( lj % 2 ) {
 			const int pk_ph = 1000;
-			ADD_PX(fs, ss, 10.0*poisson_noise(pk_ph));
-			ADD_PX(fs-1, ss, 10.0*poisson_noise(pk_ph));
-			ADD_PX(fs+1, ss, 10.0*poisson_noise(pk_ph));
-			ADD_PX(fs, ss-1, 10.0*poisson_noise(pk_ph));
-			ADD_PX(fs, ss+1, 10.0*poisson_noise(pk_ph));
+			ADD_PX(fs, ss, 10.0*poisson_noise(rng, pk_ph));
+			ADD_PX(fs-1, ss, 10.0*poisson_noise(rng, pk_ph));
+			ADD_PX(fs+1, ss, 10.0*poisson_noise(rng, pk_ph));
+			ADD_PX(fs, ss-1, 10.0*poisson_noise(rng, pk_ph));
+			ADD_PX(fs, ss+1, 10.0*poisson_noise(rng, pk_ph));
 			n_strong++;
 		} else {
 			/* Absent peak */
@@ -226,6 +233,7 @@ int main(int argc, char *argv[])
 	free(image.det);
 	free(image.dp[0]);
 	free(image.dp);
+	gsl_rng_free(rng);
 
 	if ( fail ) return 1;
 
diff --git a/tests/ring_check.c b/tests/ring_check.c
index ac871581..f7e8c473 100644
--- a/tests/ring_check.c
+++ b/tests/ring_check.c
@@ -1,10 +1,14 @@
 /*
- * integration_check.c
+ * ring_check.c
  *
  * Check peak integration
  *
- * Copyright © 2012 Thomas White <taw@physics.org>
- * Copyright © 2012 Andrew Martin <andrew.martin@desy.de>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2011-2014 Thomas White <taw@physics.org>
+ *   2012      Andrew Martin <andrew.martin@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -41,7 +45,7 @@
 /* The third integration check draws a Poisson background and checks that, on
  * average, it gets subtracted by the background subtraction. */
 static void third_integration_check(struct image *image, int n_trials,
-                                    int *fail)
+                                    int *fail, gsl_rng *rng)
 {
 	double mean_intensity = 0.0;
 	double mean_bg = 0.0;
@@ -59,7 +63,8 @@ static void third_integration_check(struct image *image, int n_trials,
 
 		for ( fs=0; fs<image->width; fs++ ) {
 		for ( ss=0; ss<image->height; ss++ ) {
-			image->data[fs+image->width*ss] = poisson_noise(1000.0);
+			image->data[fs+image->width*ss]
+			                           = poisson_noise(rng, 1000.0);
 		}
 		}
 
@@ -100,7 +105,7 @@ static void third_integration_check(struct image *image, int n_trials,
  * top of it, then checks that the intensity of the peak is correctly recovered
  * accounting for the background. */
 static void fourth_integration_check(struct image *image, int n_trials,
-                                     int *fail)
+                                     int *fail, gsl_rng *rng)
 {
 	double mean_intensity = 0.0;
 	double mean_sigma = 0.0;
@@ -118,7 +123,7 @@ static void fourth_integration_check(struct image *image, int n_trials,
 		for ( fs=0; fs<image->width; fs++ ) {
 		for ( ss=0; ss<image->height; ss++ ) {
 			int idx = fs+image->width*ss;
-			image->data[idx] = poisson_noise(1000.0);
+			image->data[idx] = poisson_noise(rng, 1000.0);
 			if ( (fs-64)*(fs-64) + (ss-64)*(ss-64) > 9*9 ) continue;
 			image->data[idx] += 1000.0;
 			pcount++;
@@ -162,16 +167,19 @@ int main(int argc, char *argv[])
 	double fsp, ssp, intensity, sigma;
 	int fs, ss;
 	FILE *fh;
-	unsigned int seed;
+	unsigned long int seed;
 	int fail = 0;
 	const int n_trials = 100;
 	int r, npx;
 	double ex;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	fh = fopen("/dev/urandom", "r");
 	fread(&seed, sizeof(seed), 1, fh);
 	fclose(fh);
-	srand(seed);
+	gsl_rng_set(rng, seed);
 
 	image.data = malloc(128*128*sizeof(float));
 	image.flags = NULL;
@@ -261,10 +269,10 @@ int main(int argc, char *argv[])
 	}
 
 	/* Third check: Poisson background should get mostly subtracted */
-	third_integration_check(&image, n_trials, &fail);
+	third_integration_check(&image, n_trials, &fail, rng);
 
 	/* Fourth check: peak on Poisson background */
-	fourth_integration_check(&image, n_trials, &fail);
+	fourth_integration_check(&image, n_trials, &fail, rng);
 
 	/* Fifth check: uniform peak on uniform background */
 	npx = 0;
@@ -307,6 +315,7 @@ int main(int argc, char *argv[])
 	free(image.det->panels);
 	free(image.det);
 	free(image.data);
+	gsl_rng_free(rng);
 
 	if ( fail ) return 1;
 
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
index 7d25aa04..3affb652 100644
--- a/tests/transformation_check.c
+++ b/tests/transformation_check.c
@@ -3,7 +3,11 @@
  *
  * Check that unit cell transformations work
  *
- * Copyright © 2012 Thomas White <taw@physics.org>
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2012-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -111,6 +115,9 @@ int main(int argc, char *argv[])
 	int fail = 0;
 	UnitCell *cell, *cref;
 	UnitCellTransformation *tfn;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	cref = cell_new_from_parameters(50e-10, 55e-10, 70e-10,
 	                                deg2rad(67.0),
@@ -118,7 +125,7 @@ int main(int argc, char *argv[])
 	                                deg2rad(77.0));
 	if ( cref == NULL ) return 1;
 
-	cell = cell_rotate(cref, random_quaternion());
+	cell = cell_rotate(cref, random_quaternion(rng));
 	if ( cell == NULL ) return 1;
 	cell_free(cref);
 
@@ -171,6 +178,7 @@ int main(int argc, char *argv[])
 	tfn_free(tfn);
 
 	cell_free(cell);
+	gsl_rng_free(rng);
 
 	return fail;
 }