RNG overhaul

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

21 files changed, 277 insertions, 158 deletions

diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index 5cd0ce0f..174abefe 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -311,7 +311,7 @@ double get_tt(struct image *image, double fs, double ss, int *err)
 }
 
 
-void record_image(struct image *image, int do_poisson)
+void record_image(struct image *image, int do_poisson, gsl_rng *rng)
 {
 	int x, y;
 	double total_energy, energy_density;
@@ -371,7 +371,7 @@ void record_image(struct image *image, int do_poisson)
 		sa = proj_area / (dsq + Lsq);
 
 		if ( do_poisson ) {
-			counts = poisson_noise(intensity * ph_per_e * sa);
+			counts = poisson_noise(rng, intensity * ph_per_e * sa);
 		} else {
 			cf = intensity * ph_per_e * sa;
 			counts = cf;
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index cfe0808a..eaa5618e 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -153,7 +153,7 @@ extern double get_tt(struct image *image, double xs, double ys, int *err);
 
 extern int in_bad_region(struct detector *det, double fs, double ss);
 
-extern void record_image(struct image *image, int do_poisson);
+extern void record_image(struct image *image, int do_poisson, gsl_rng *rng);
 
 extern struct panel *find_panel(struct detector *det, double fs, double ss);
 
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 2ab4397a..110d8701 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -793,6 +793,13 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
 	double bx, by, bz;
 	double cx, cy, cz;
 
+	gsl_rng *rng;
+
+	/* This is a rubbish RNG, but it serves for this purpose: nothing more
+	 * than "I couldn't be bothered to think of n_trials sets of random
+	 * indices". */
+	rng = gsl_rng_alloc(gsl_rng_taus2);
+
 	cell_get_reciprocal(c1, &asx, &asy, &asz,
 	                        &bsx, &bsy, &bsz,
 	                        &csx, &csy, &csz);
@@ -807,9 +814,9 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
 		double dev;
 		signed int h2i, k2i, l2i;
 
-		h = flat_noise(0, 10);
-		k = flat_noise(0, 10);
-		l = flat_noise(0, 10);
+		h = gsl_rng_uniform_int(rng, 10);
+		k = gsl_rng_uniform_int(rng, 10);
+		l = gsl_rng_uniform_int(rng, 10);
 
 		/* Position of this (randomly selected)
 		 * reciprocal lattice point */
@@ -844,6 +851,8 @@ static int check_twinning(UnitCell *c1, UnitCell *c2, int verbose)
 		STATUS("%i duplicates.\n", n_dup);
 	}
 
+	gsl_rng_free(rng);
+
 	if ( n_dup > 10 ) return 1;
 	return 0;
 }
diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
index 2f5107e3..4c85bf9d 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -3,11 +3,11 @@
  *
  * Utility stuff
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2012 Thomas White <taw@physics.org>
+ *   2009-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -165,26 +165,26 @@ void progress_bar(int val, int total, const char *text)
 }
 
 
-double random_flat(double max)
+double random_flat(gsl_rng *rng, double max)
 {
-	return max * (double)random()/RAND_MAX;
+	return max * gsl_rng_uniform(rng);
 }
 
 
-double flat_noise(double expected, double width)
+double flat_noise(gsl_rng *rng, double expected, double width)
 {
-	double noise = random_flat(2.0*width);
+	double noise = random_flat(rng, 2.0*width);
 	return expected+noise-width;
 }
 
 
-double gaussian_noise(double expected, double stddev)
+double gaussian_noise(gsl_rng *rng, double expected, double stddev)
 {
 	double x1, x2, noise;
 
 	/* A uniformly distributed random number between 0 and 1 */
-	x1 = ((double)random()/RAND_MAX);
-	x2 = ((double)random()/RAND_MAX);
+	x1 = gsl_rng_uniform(rng);
+	x2 = gsl_rng_uniform(rng);
 
 	noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2);
 
@@ -192,14 +192,14 @@ double gaussian_noise(double expected, double stddev)
 }
 
 
-static int fake_poisson_noise(double expected)
+static int fake_poisson_noise(gsl_rng *rng, double expected)
 {
-	double rf = gaussian_noise(expected, sqrt(expected));
+	double rf = gaussian_noise(rng, expected, sqrt(expected));
 	return (int)rf;
 }
 
 
-int poisson_noise(double expected)
+int poisson_noise(gsl_rng *rng, double expected)
 {
 	double L;
 	int k = 0;
@@ -207,7 +207,7 @@ int poisson_noise(double expected)
 
 	/* For large values of the mean, we get big problems with arithmetic.
 	 * In such cases, fall back on a Gaussian with the right variance. */
-	if ( expected > 100.0 ) return fake_poisson_noise(expected);
+	if ( expected > 100.0 ) return fake_poisson_noise(rng, expected);
 
 	L = exp(-expected);
 
@@ -216,7 +216,7 @@ int poisson_noise(double expected)
 		double r;
 
 		k++;
-		r = (double)random()/(double)RAND_MAX;
+		r = gsl_rng_uniform(rng);
 		p *= r;
 
 	} while ( p > L );
@@ -282,14 +282,14 @@ struct quaternion normalise_quaternion(struct quaternion q)
  *
  * Returns: a randomly generated, normalised, quaternion.
  **/
-struct quaternion random_quaternion()
+struct quaternion random_quaternion(gsl_rng *rng)
 {
 	struct quaternion q;
 
-	q.w = 2.0*(double)random()/RAND_MAX - 1.0;
-	q.x = 2.0*(double)random()/RAND_MAX - 1.0;
-	q.y = 2.0*(double)random()/RAND_MAX - 1.0;
-	q.z = 2.0*(double)random()/RAND_MAX - 1.0;
+	q.w = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.x = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.y = 2.0*gsl_rng_uniform(rng) - 1.0;
+	q.z = 2.0*gsl_rng_uniform(rng) - 1.0;
 	q = normalise_quaternion(q);
 
 	return q;
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index f4477717..81270916 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -3,11 +3,11 @@
  *
  * Utility stuff
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2012 Thomas White <taw@physics.org>
+ *   2009-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -40,6 +40,7 @@
 #include <pthread.h>
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_vector.h>
+#include <gsl/gsl_rng.h>
 
 
 #include "thread-pool.h"
@@ -89,7 +90,7 @@ struct quaternion {
 
 extern struct quaternion normalise_quaternion(struct quaternion q);
 extern double quaternion_modulus(struct quaternion q);
-extern struct quaternion random_quaternion(void);
+extern struct quaternion random_quaternion(gsl_rng *rng);
 extern int quaternion_valid(struct quaternion q);
 extern struct rvec quat_rot(struct rvec q, struct quaternion z);
 
@@ -114,10 +115,10 @@ 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 random_flat(double max);
-extern double flat_noise(double expected, double width);
-extern double gaussian_noise(double expected, double stddev);
-extern int poisson_noise(double expected);
+extern double random_flat(gsl_rng *rng, double max);
+extern double flat_noise(gsl_rng *rng, double expected, double width);
+extern double gaussian_noise(gsl_rng *rng, double expected, double stddev);
+extern int poisson_noise(gsl_rng *rng, double expected);
 
 /* Keep these ones inline, to avoid function call overhead */
 static inline struct quaternion invalid_quaternion(void)
diff --git a/src/diffraction.c b/src/diffraction.c
index 826aac77..92e18521 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -456,7 +456,7 @@ static struct sample *get_gaussian_spectrum(double eV_cen, double eV_step,
 }
 
 
-static int add_sase_noise(struct sample *spectrum, int nsteps)
+static int add_sase_noise(struct sample *spectrum, int nsteps, gsl_rng *rng)
 {
 	struct sample *noise;
 	int i, j;
@@ -483,7 +483,7 @@ static int add_sase_noise(struct sample *spectrum, int nsteps)
 		noise[i].weight = 0.0;
 
 		/* Gaussian noise with mean = 0, std = 1 */
-		gaussianNoise[i] = gaussian_noise(noise_mean, noise_sigma);
+		gaussianNoise[i] = gaussian_noise(rng, noise_mean, noise_sigma);
 		gaussianNoise[i+nsteps] = gaussianNoise[i];
 		gaussianNoise[i+2*nsteps] = gaussianNoise[i];
 	}
@@ -546,7 +546,7 @@ struct sample *generate_tophat(struct image *image)
 }
 
 
-struct sample *generate_SASE(struct image *image)
+struct sample *generate_SASE(struct image *image, gsl_rng *rng)
 {
 	struct sample *spectrum;
 	int i;
@@ -555,7 +555,7 @@ struct sample *generate_SASE(struct image *image)
 	const double jitter_sigma_eV = 8.0;
 
 	/* Central wavelength jitters with Gaussian distribution */
-	eV_cen = gaussian_noise(ph_lambda_to_eV(image->lambda),
+	eV_cen = gaussian_noise(rng, ph_lambda_to_eV(image->lambda),
 	                        jitter_sigma_eV);
 
 	/* Convert FWHM to standard deviation.  Note that bandwidth is taken to
@@ -571,7 +571,7 @@ struct sample *generate_SASE(struct image *image)
 	spectrum = get_gaussian_spectrum(eV_cen, eV_step, sigma, spec_size);
 
 	/* Add SASE-type noise to Gaussian spectrum */
-	add_sase_noise(spectrum, spec_size);
+	add_sase_noise(spectrum, spec_size, rng);
 
 	/* Normalise intensity (before taking restricted number of samples) */
 	double total_weight = 0.0;
diff --git a/src/diffraction.h b/src/diffraction.h
index a903346e..5b67c198 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -35,10 +35,13 @@
 #ifndef DIFFRACTION_H
 #define DIFFRACTION_H
 
+#include <gsl/gsl_rng.h>
+
 #include "image.h"
 #include "cell.h"
 #include "symmetry.h"
 
+
 typedef enum {
 	GRADIENT_MOSAIC,
 	GRADIENT_INTERPOLATE,
@@ -52,6 +55,6 @@ extern void get_diffraction(struct image *image, int na, int nb, int nc,
 
 extern struct sample *generate_tophat(struct image *image);
 
-extern struct sample *generate_SASE(struct image *image);
+extern struct sample *generate_SASE(struct image *image, gsl_rng *rng);
 
 #endif	/* DIFFRACTION_H */
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 454737b0..8ba77ebf 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -97,6 +97,9 @@ static void poisson_reflections(RefList *list, double adu_per_photon)
 {
 	Reflection *refl;
 	RefListIterator *iter;
+	gsl_rng *rng;
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
@@ -106,10 +109,12 @@ static void poisson_reflections(RefList *list, double adu_per_photon)
 
 		val = get_intensity(refl);
 
-		c = adu_per_photon * poisson_noise(val/adu_per_photon);
+		c = adu_per_photon * poisson_noise(rng, val/adu_per_photon);
 		set_intensity(refl, c);
 
 	}
+
+	gsl_rng_free(rng);
 }
 
 
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 7a3517a3..d9c9ae9f 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -3,11 +3,11 @@
  *
  * Generate partials for testing scaling
  *
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2013 Thomas White <taw@physics.org>
+ *   2011-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -39,6 +39,7 @@
 #include <getopt.h>
 #include <assert.h>
 #include <pthread.h>
+#include <gsl/gsl_rng.h>
 
 #include "image.h"
 #include "utils.h"
@@ -54,7 +55,7 @@
 #define NBINS 50
 
 
-static void mess_up_cell(Crystal *cr, double cnoise)
+static void mess_up_cell(Crystal *cr, double cnoise, gsl_rng *rng)
 {
 	double ax, ay, az;
 	double bx, by, bz;
@@ -65,15 +66,15 @@ static void mess_up_cell(Crystal *cr, double cnoise)
 	//cell_print(cell);
 
 	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	ax = flat_noise(ax, cnoise*fabs(ax)/100.0);
-	ay = flat_noise(ay, cnoise*fabs(ay)/100.0);
-	az = flat_noise(az, cnoise*fabs(az)/100.0);
-	bx = flat_noise(bx, cnoise*fabs(bx)/100.0);
-	by = flat_noise(by, cnoise*fabs(by)/100.0);
-	bz = flat_noise(bz, cnoise*fabs(bz)/100.0);
-	cx = flat_noise(cx, cnoise*fabs(cx)/100.0);
-	cy = flat_noise(cy, cnoise*fabs(cy)/100.0);
-	cz = flat_noise(cz, cnoise*fabs(cz)/100.0);
+	ax = flat_noise(rng, ax, cnoise*fabs(ax)/100.0);
+	ay = flat_noise(rng, ay, cnoise*fabs(ay)/100.0);
+	az = flat_noise(rng, az, cnoise*fabs(az)/100.0);
+	bx = flat_noise(rng, bx, cnoise*fabs(bx)/100.0);
+	by = flat_noise(rng, by, cnoise*fabs(by)/100.0);
+	bz = flat_noise(rng, bz, cnoise*fabs(bz)/100.0);
+	cx = flat_noise(rng, cx, cnoise*fabs(cx)/100.0);
+	cy = flat_noise(rng, cy, cnoise*fabs(cy)/100.0);
+	cz = flat_noise(rng, cz, cnoise*fabs(cz)/100.0);
 	cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
 
 	//STATUS("Changed:\n");
@@ -89,7 +90,7 @@ static void calculate_partials(Crystal *cr,
                                pthread_rwlock_t *full_lock,
                                unsigned long int *n_ref, double *p_hist,
                                double *p_max, double max_q, double full_stddev,
-                               double noise_stddev)
+                               double noise_stddev, gsl_rng *rng)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -124,7 +125,7 @@ static void calculate_partials(Crystal *cr,
 				rfull = find_refl(full, h, k, l);
 				if ( rfull == NULL ) {
 					rfull = add_refl(full, h, k, l);
-					If = fabs(gaussian_noise(0.0,
+					If = fabs(gaussian_noise(rng, 0.0,
 					                         full_stddev));
 					set_intensity(rfull, If);
 					set_redundancy(rfull, 1);
@@ -160,7 +161,7 @@ static void calculate_partials(Crystal *cr,
 			       res, bin,  p);
 		}
 
-		Ip = gaussian_noise(Ip, noise_stddev);
+		Ip = gaussian_noise(rng, Ip, noise_stddev);
 
 		set_intensity(refl, Ip);
 		set_esd_intensity(refl, noise_stddev);
@@ -226,6 +227,7 @@ struct queue_args
 	double p_max[NBINS];
 
 	Stream *stream;
+	gsl_rng **rngs;
 };
 
 
@@ -280,14 +282,15 @@ static void run_job(void *vwargs, int cookie)
 	crystal_set_image(cr, &wargs->image);
 
 	do {
-		osf = gaussian_noise(1.0, qargs->osf_stddev);
+		osf = gaussian_noise(qargs->rngs[cookie], 1.0,
+		                     qargs->osf_stddev);
 	} while ( osf <= 0.0 );
 	crystal_set_osf(cr, osf);
 	crystal_set_mosaicity(cr, 0.0);
 	crystal_set_profile_radius(cr, wargs->image.beam->profile_radius);
 
 	/* Set up a random orientation */
-	orientation = random_quaternion();
+	orientation = random_quaternion(qargs->rngs[cookie]);
 	crystal_set_cell(cr, cell_rotate(qargs->cell, orientation));
 
 	snprintf(wargs->image.filename, 255, "dummy.h5");
@@ -305,10 +308,10 @@ static void run_job(void *vwargs, int cookie)
 	                   &qargs->full_lock,
 	                   wargs->n_ref, wargs->p_hist, wargs->p_max,
 	                   qargs->max_q, qargs->full_stddev,
-	                   qargs->noise_stddev);
+	                   qargs->noise_stddev, qargs->rngs[cookie]);
 
 	/* Give a slightly incorrect cell in the stream */
-	mess_up_cell(cr, qargs->cnoise);
+	mess_up_cell(cr, qargs->cnoise, qargs->rngs[cookie]);
 
 	image_add_crystal(&wargs->image, cr);
 }
@@ -367,6 +370,7 @@ int main(int argc, char *argv[])
 	double osf_stddev = 2.0;
 	double full_stddev = 1000.0;
 	double noise_stddev = 20.0;
+	gsl_rng *rng_for_seeds;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -623,6 +627,18 @@ int main(int argc, char *argv[])
 	qargs.noise_stddev = noise_stddev;
 	qargs.max_q = largest_q(&image);
 
+	qargs.rngs = malloc(n_threads * sizeof(gsl_rng *));
+	if ( qargs.rngs == NULL ) {
+		ERROR("Failed to allocate RNGs\n");
+		return 1;
+	}
+	rng_for_seeds = gsl_rng_alloc(gsl_rng_mt19937);
+	for ( i=0; i<n_threads; i++ ) {
+		qargs.rngs[i] = gsl_rng_alloc(gsl_rng_mt19937);
+		gsl_rng_set(qargs.rngs[i], gsl_rng_get(rng_for_seeds));
+	}
+	gsl_rng_free(rng_for_seeds);
+
 	for ( i=0; i<NBINS; i++ ) {
 		qargs.n_ref[i] = 0;
 		qargs.p_hist[i] = 0.0;
@@ -683,6 +699,9 @@ int main(int argc, char *argv[])
 
 	}
 
+	for ( i=0; i<n_threads; i++ ) {
+		gsl_rng_free(qargs.rngs[i]);
+	}
 	pthread_rwlock_destroy(&qargs.full_lock);
 	close_stream(stream);
 	cell_free(cell);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 81a94c37..5ea328bf 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -3,11 +3,13 @@
  *
  * Simulate diffraction patterns from small crystals
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2009-2012 Thomas White <taw@physics.org>
+ *   2009-2014 Thomas White <taw@physics.org>
+ *   2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de>
+ *   2013      Alexandra Tolstikova
  *
  * This file is part of CrystFEL.
  *
@@ -278,6 +280,7 @@ int main(int argc, char *argv[])
 	char *sym_str = NULL;
 	SymOpList *sym;
 	int nsamples = 3;
+	gsl_rng *rng;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -407,15 +410,6 @@ int main(int argc, char *argv[])
 
 	}
 
-	if ( config_random ) {
-		FILE *fh;
-		unsigned int seed;
-		fh = fopen("/dev/urandom", "r");
-		fread(&seed, sizeof(seed), 1, fh);
-		fclose(fh);
-		srand(seed);
-	}
-
 	if ( random_size == 1 ) {
 		ERROR("You must specify both --min-size and --max-size.\n");
 		return 1;
@@ -565,6 +559,16 @@ int main(int argc, char *argv[])
 	image.features = NULL;
 	image.flags = NULL;
 
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+	if ( config_random ) {
+		FILE *fh;
+		unsigned long int seed;
+		fh = fopen("/dev/urandom", "r");
+		fread(&seed, sizeof(seed), 1, fh);
+		fclose(fh);
+		gsl_rng_set(rng, seed);
+	}
+
 	powder = calloc(image.width*image.height, sizeof(*powder));
 
 	/* Splurge a few useful numbers */
@@ -597,7 +601,7 @@ int main(int argc, char *argv[])
 
 		/* Read quaternion from stdin */
 		if ( config_randomquat ) {
-			orientation = random_quaternion();
+			orientation = random_quaternion(rng);
 		} else {
 			orientation = read_quaternion();
 		}
@@ -620,7 +624,7 @@ int main(int argc, char *argv[])
 			break;
 
 			case SPECTRUM_SASE :
-			image.spectrum = generate_SASE(&image);
+			image.spectrum = generate_SASE(&image, rng);
 			break;
 
 		}
@@ -650,7 +654,7 @@ int main(int argc, char *argv[])
 			goto skip;
 		}
 
-		record_image(&image, !config_nonoise);
+		record_image(&image, !config_nonoise, rng);
 
 		if ( powder_fn != NULL ) {
 
@@ -723,5 +727,7 @@ skip:
 	free(sym_str);
 	free_symoplist(sym);
 
+	gsl_rng_free(rng);
+
 	return 0;
 }
diff --git a/src/scaling-report.c b/src/scaling-report.c
index 89b2f3f2..85ed0ad2 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -196,6 +196,7 @@ static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
 	double pcalcmin[nbins];
 	double pcalcmax[nbins];
 	int num_nondud;
+	gsl_rng *rng;
 
 	show_text_simple(cr, "Observed partiality", -20.0, g_height/2.0,
 	                      NULL, -M_PI_2, J_CENTER);
@@ -223,6 +224,10 @@ static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
 		num_nondud++;
 	}
 
+	/* The reflections chosen for the graph will be the same every time
+	 * (given the same sequence of input reflections, scalabilities etc) */
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+
 	cairo_set_source_rgb(cr, 0.0, 0.7, 0.0);
 	prob = 1.0 / num_nondud;
 	for ( i=0; i<n; i++ ) {
@@ -279,13 +284,15 @@ static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
 
 			bin = nbins * pcalc;
 
-			if ( random_flat(1.0) < prob ) {
+			if ( random_flat(rng, 1.0) < prob ) {
 				plot_point(cr, g_width, g_height, pcalc, pobs);
 			}
 		}
 
 	}
 
+	gsl_rng_free(rng);
+
 	cairo_new_path(cr);
 	cairo_rectangle(cr, 0.0, 0.0, g_width, g_height);
 	cairo_clip(cr);
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;
 }