Separate partiality calculation from prediction

4 files changed, 95 insertions, 60 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index b968dc06..576655de 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -266,15 +266,12 @@ static double partiality(PartialityModel pmodel,
 
 
 static Reflection *check_reflection(struct image *image, Crystal *cryst,
-                                    PartialityModel pmodel,
                                     signed int h, signed int k, signed int l,
                                     double xl, double yl, double zl,
                                     Reflection *updateme)
 {
-	const int output = 0;
 	double tl;
-	double rlow, rhigh;     /* "Excitation error" */
-	double part;            /* Partiality */
+	double rlow, rhigh;    /* "Excitation error" */
 	double klow, khigh;    /* Wavenumber */
 	Reflection *refl;
 	double cet, cez;  /* Centre of Ewald sphere */
@@ -312,20 +309,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	     && (fabs(rlow) > pr)
 	     && (fabs(rhigh) > pr) ) return NULL;
 
-	/* Calculate partiality */
-	part = partiality(pmodel, h, k, l, image->serial, rlow, rhigh, pr);
-
-	if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) {
-		ERROR("Assigning NAN partiality!\n");
-		ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
-		      h, k, l, rlow, rhigh);
-		ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw);
-		ERROR("D = %e\n", rlow - rhigh);
-		ERROR("klow = %e  khigh = %e\n", klow, khigh);
-		ERROR("distance(..) = %e\n", distance(cet, cez, tl, zl));
-		abort();
-	}
-
 	if ( updateme == NULL ) {
 		refl = reflection_new(h, k, l);
 	} else {
@@ -373,15 +356,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		}
 	}
 
-	set_partial(refl, rlow, rhigh, part);
+	set_partial(refl, rlow, rhigh, 1.0); /* Actual partiality set by
+	                                      * calculate_partialities() */
 	set_lorentz(refl, 1.0);
 	set_symmetric_indices(refl, h, k, l);
 	set_redundancy(refl, 1);
 
-	if ( output ) {
-		printf("%3i %3i %3i %6f %5.2f\n", h, k, l, 0.0, part);
-	}
-
 	return refl;
 }
 
@@ -472,15 +452,7 @@ double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
 }
 
 
-RefList *find_intersections(struct image *image, Crystal *cryst,
-                            PartialityModel pmodel)
-{
-	return find_intersections_to_res(image, cryst, pmodel, INFINITY);
-}
-
-
-RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
-                                   PartialityModel pmodel, double max_res)
+RefList *predict_to_res(Crystal *cryst, double max_res)
 {
 	double ax, ay, az;
 	double bx, by, bz;
@@ -509,7 +481,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
-	mres = largest_q(image);
+	mres = largest_q(crystal_get_image(cryst));
 	if ( mres > max_res ) mres = max_res;
 
 	hmax = mres * modulus(ax, ay, az);
@@ -544,7 +516,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
 		yl = h*asy + k*bsy + l*csy;
 		zl = h*asz + k*bsz + l*csz;
 
-		refl = check_reflection(image, cryst, pmodel,
+		refl = check_reflection(crystal_get_image(cryst), cryst,
 		                        h, k, l, xl, yl, zl, NULL);
 
 		if ( refl != NULL ) {
@@ -559,12 +531,12 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
 }
 
 
-static void set_unity_partialities(Crystal *cryst)
+static void set_unity_partialities(RefList *list)
 {
 	Reflection *refl;
 	RefListIterator *iter;
 
-	for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
+	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
@@ -575,14 +547,85 @@ static void set_unity_partialities(Crystal *cryst)
 
 
 /**
- * update_partialities:
+ * calculate_partialities:
  * @cryst: A %Crystal
  * @pmodel: A %PartialityModel
  *
- * Updates the partialities of @cryst's reflections according to @pmodel and
+ * Calculates the partialities for the reflections in @cryst, given the current
+ * crystal and image parameters.  The crystal's image and reflection lists
+ * must be set.  The specified %PartialityModel will be used.
+ *
+ * You must not have changed the crystal or image parameters since you last
+ * called predict_to_res() or update_predictions(), because this function
+ * relies on the limiting wavelength values calculated by those functions.
+ */
+void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
+{
+	RefList *list;
+	Reflection *refl;
+	RefListIterator *iter;
+	double pr;
+	struct image *image;
+
+	list = crystal_get_reflections(cryst);
+	if ( list == NULL ) {
+		ERROR("No reflections for partiality calculation!\n");
+		return;
+	}
+
+	if ( pmodel == PMODEL_UNITY ) {
+		set_unity_partialities(list);
+		return;
+	}
+
+	image = crystal_get_image(cryst);
+	if ( image == NULL ) {
+		ERROR("No image for partiality calculation!\n");
+		return;
+	}
+
+	pr = crystal_get_profile_radius(cryst);
+
+	for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		double rlow, rhigh, part;
+		signed int h, k, l;
+
+		get_symmetric_indices(refl, &h, &k, &l);
+		get_partial(refl, &rlow, &rhigh, &part);
+
+		part = partiality(pmodel, h, k, l, image->serial,
+		                  rlow, rhigh, pr);
+
+		if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) {
+			ERROR("Assigning NAN partiality!\n");
+			ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
+			      h, k, l, rlow, rhigh);
+			ERROR("R = %e, bw = %e\n", pr, image->bw);
+			ERROR("D = %e\n", rlow - rhigh);
+			abort();
+	        }
+
+	        set_partiality(refl, part);
+
+	}
+}
+
+
+/**
+ * update_predictions:
+ * @cryst: A %Crystal
+ *
+ * Updates the predicted reflections (positions and excitation errors, but not
+ * the actual partialities) of @cryst's reflections according to
  * the current state of the crystal (e.g. its unit cell parameters).
+ *
+ * If you need to update the partialities as well, call calculate_partialities()
+ * afterwards.
  */
-void update_partialities(Crystal *cryst, PartialityModel pmodel)
+void update_predictions(Crystal *cryst)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -591,11 +634,6 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel)
 	double csx, csy, csz;
 	struct image *image = crystal_get_image(cryst);
 
-	if ( pmodel == PMODEL_UNITY ) {
-		set_unity_partialities(cryst);
-		return;
-	}
-
 	cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz,
 	                    &bsx, &bsy, &bsz, &csx, &csy, &csz);
 
@@ -613,8 +651,8 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel)
 		yl = h*asy + k*bsy + l*csy;
 		zl = h*asz + k*bsz + l*csz;
 
-		check_reflection(image, cryst, pmodel,
-		                 h, k, l, xl, yl, zl, refl);
+		check_reflection(image, cryst, h, k, l,
+		                 xl, yl, zl, refl);
 
 	}
 }
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 6b035175..4fbd0bca 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,12 +3,12 @@
  *
  * Geometry of diffraction
  *
- * Copyright © 2013-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2016 Thomas White <taw@physics.org>
  *   2012      Richard Kirian
  *
  * This file is part of CrystFEL.
@@ -85,15 +85,13 @@ enum gparam {
 };
 
 
-extern RefList *find_intersections(struct image *image, Crystal *cryst,
-                                   PartialityModel pmodel);
-extern RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
-                                          PartialityModel pmodel,
-					  double max_res);
+extern RefList *predict_to_res(Crystal *cryst, double max_res);
+
+extern void calculate_partialities(Crystal *cryst, PartialityModel pmodel);
 
 extern double r_gradient(UnitCell *cell, int k, Reflection *refl,
                          struct image *image);
-extern void update_partialities(Crystal *cryst, PartialityModel pmodel);
+extern void update_predictions(Crystal *cryst);
 extern void polarisation_correction(RefList *list, UnitCell *cell,
                                     struct image *image);
 
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index d50f8d6b..96eec92a 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1673,8 +1673,7 @@ void integrate_all_4(struct image *image, IntegrationMethod meth,
 		                          image->features);
 		crystal_set_resolution_limit(image->crystals[i], res);
 
-		list = find_intersections_to_res(image, image->crystals[i],
-		                                 pmodel, res+push_res);
+		list = predict_to_res(image->crystals[i], res+push_res);
 		crystal_set_reflections(image->crystals[i], list);
 
 	}
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index e702cf33..3246dbec 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -262,7 +262,7 @@ static int pair_peaks(struct image *image, Crystal *cr,
 		 * doesn't fall on this panel.  We're only interested
 		 * in how far away it is from the peak location.
 		 * The predicted position and excitation errors will be
-		 * filled in by update_partialities(). */
+		 * filled in by update_predictions(). */
 		set_panel(refl, f->p);
 
 		rps[n].refl = refl;
@@ -275,7 +275,7 @@ static int pair_peaks(struct image *image, Crystal *cr,
 	/* Get the excitation errors and detector positions for the candidate
 	 * reflections */
 	crystal_set_reflections(cr, all_reflist);
-	update_partialities(cr, PMODEL_SCSPHERE);
+	update_predictions(cr);
 
 	/* Pass over the peaks again, keeping only the ones which look like
 	 * good pairings */
@@ -339,7 +339,7 @@ void refine_radius(Crystal *cr, struct image *image)
 		return;
 	}
 	crystal_set_reflections(cr, reflist);
-	update_partialities(cr, PMODEL_SCSPHERE);
+	update_predictions(cr);
 	crystal_set_reflections(cr, NULL);
 
 	qsort(rps, n_acc, sizeof(struct reflpeak), cmpd2);
@@ -635,7 +635,7 @@ int refine_prediction(struct image *image, Crystal *cr)
 
 	/* Refine */
 	for ( i=0; i<MAX_CYCLES; i++ ) {
-		update_partialities(cr, PMODEL_SCSPHERE);
+		update_predictions(cr);
 		if ( iterate(rps, n, crystal_get_cell(cr), image,
 		             &total_x, &total_y, &total_z) ) return 1;
 		//STATUS("Residual after %i = %e\n", i,