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 ++
 src/partialator.c          |  2 ++
 src/post-refinement.c      | 14 ++++++++------
 4 files changed, 51 insertions(+), 7 deletions(-)

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;
 
diff --git a/src/partialator.c b/src/partialator.c
index 70c82960..d539f784 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -530,6 +530,8 @@ int main(int argc, char *argv[])
 			pmodel = PMODEL_SCGAUSSIAN;
 		} else if ( strcmp(pmodel_str, "scsphere") == 0 ) {
 			pmodel = PMODEL_SCSPHERE;
+		} else if ( strcmp(pmodel_str, "random") == 0 ) {
+			pmodel = PMODEL_RANDOM;
 		} else {
 			ERROR("Unknown partiality model '%s'.\n", pmodel_str);
 			return 1;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 9054ed89..12a29cb2 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -147,10 +147,11 @@ static double volume_fraction_rgradient(double r, double pr,
 
 		case PMODEL_SCGAUSSIAN :
 		return gaussian_fraction_rgradient(r, pr);
-	}
 
-	ERROR("No pmodel in volume_fraction_rgradient!\n");
-	return 1.0;
+		default :
+		ERROR("No pmodel in volume_fraction_rgradient!\n");
+		return 1.0;
+	}
 }
 
 
@@ -167,10 +168,11 @@ static double volume_fraction(double rlow, double rhigh, double pr,
 
 		case PMODEL_SCGAUSSIAN :
 		return gaussian_fraction(rlow, rhigh, pr);
-	}
 
-	ERROR("No pmodel in volume_fraction!\n");
-	return 1.0;
+		default :
+		ERROR("No pmodel in volume_fraction!\n");
+		return 1.0;
+	}
 }
 
 
-- 
cgit v1.2.3