From a59b214e1c08f4b988cc4fb2e6ef1c5d10ec363f Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Tue, 21 Jul 2015 13:28:01 +0200
Subject: Add random partiality model

---
 libcrystfel/src/geometry.c | 40 +++++++++++++++++++++++++++++++++++++++-
 libcrystfel/src/geometry.h |  2 ++
 2 files changed, 41 insertions(+), 1 deletion(-)

(limited to 'libcrystfel/src')

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index d8b40a16..cdd936e2 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -164,7 +164,42 @@ double gaussian_fraction(double rlow, double rhigh, double R)
 }
 
 
+static double random_partiality(signed int h, signed int k, signed int l,
+                                int serial)
+{
+	gsl_rng *rng = gsl_rng_alloc(gsl_rng_mt19937);
+	unsigned long int seed;
+	double p;
+	int i;
+
+	gsl_rng_set(rng, serial);
+	for ( i=0; i<abs(h)+1; i++ ) {
+		seed = gsl_rng_get(rng);
+	}
+	if ( h < 0 ) seed = gsl_rng_get(rng);
+	gsl_rng_set(rng, seed);
+
+	for ( i=0; i<abs(k)+1; i++ ) {
+		seed = gsl_rng_get(rng);
+	}
+	if ( k < 0 ) seed = gsl_rng_get(rng);
+	gsl_rng_set(rng, seed);
+
+	for ( i=0; i<abs(l)+1; i++ ) {
+		seed = gsl_rng_get(rng);
+	}
+	if ( l < 0 ) seed = gsl_rng_get(rng);
+	gsl_rng_set(rng, seed);
+
+	p = gsl_rng_uniform(rng);
+	gsl_rng_free(rng);
+	return p;
+}
+
+
 static double partiality(PartialityModel pmodel,
+                         signed int h, signed int k, signed int l,
+                         int serial,
                          double rlow, double rhigh, double pr)
 {
 	double D = rlow - rhigh;
@@ -182,6 +217,9 @@ static double partiality(PartialityModel pmodel,
 		case PMODEL_SCGAUSSIAN:
 		return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D);
 
+		case PMODEL_RANDOM:
+		return random_partiality(h, k, l, serial);
+
 	}
 }
 
@@ -234,7 +272,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	     && (fabs(rhigh) > pr) ) return NULL;
 
 	/* Calculate partiality */
-	part = partiality(pmodel, rlow, rhigh, pr);
+	part = partiality(pmodel, h, k, l, image->serial, rlow, rhigh, pr);
 
 	if ( isnan(part) ) {
 		ERROR("Assigning NAN partiality!\n");
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 474755b2..6b035175 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -49,6 +49,7 @@ extern "C" {
  * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
  * @PMODEL_SCSPHERE : Sphere model with source coverage factor included
  * @PMODEL_SCGAUSSIAN : Gaussian model with source coverage factor included
+ * @PMODEL_RANDOM : Randomly assigned partialities
  *
  * A %PartialityModel describes a geometrical model which can be used to
  * calculate spot partialities and Lorentz correction factors.
@@ -58,6 +59,7 @@ typedef enum {
 	PMODEL_UNITY,
 	PMODEL_SCSPHERE,
 	PMODEL_SCGAUSSIAN,
+	PMODEL_RANDOM,
 
 } PartialityModel;
 
-- 
cgit v1.2.3