Factorise correction of intensity for G, B, p and L

There were no fewer than 8 different places in the code where these
factors needed to be applied.  They all need to agree on conventions
such as whether the intensities in the pattern should be multiplied or
divided by G to "correct" them.  They were not.

This commit reduces the number of places to three: one function
(actually two functions, so that the value before partiality can also be
calculated consistently), plus log_residual and scale_one_crystal, where
slightly different logarithmic versions are used instead.

5 files changed, 39 insertions, 33 deletions

diff --git a/src/merge.c b/src/merge.c
index 787afca3..8017182b 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -187,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie)
 		Reflection *f;
 		signed int h, k, l;
 		double mean, sumweight, M2, temp, delta, R;
-		double corr, res, w;
+		double res, w;
 		struct reflection_contributions *c;
 
 		if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
@@ -217,25 +217,16 @@ static void run_merge_job(void *vwargs, int cookie)
 			continue;
 		}
 
-		/* Total (multiplicative) correction factor */
-		corr = 1.0/G * exp(B*res*res) * get_lorentz(refl)
-		        / get_partiality(refl);
-		if ( isnan(corr) ) {
-			ERROR("NaN corr:\n");
-			ERROR("G = %f, B = %e\n", G, B);
-			ERROR("res = %e\n", res);
-			ERROR("p = %f\n", get_partiality(refl));
-			unlock_reflection(f);
-			continue;
-		}
-
 		/* Reflections count less the more they have to be scaled up */
 		w = get_partiality(refl);
 
 		/* Running mean and variance calculation */
 		temp = w + sumweight;
-		if (ln_merge) delta = log(get_intensity(refl)*corr) - mean;
-		else delta = get_intensity(refl)*corr - mean;
+		if ( ln_merge ) {
+			delta = log(correct_reflection(refl, G, B, res)) - mean;
+		} else {
+			delta = correct_reflection(refl, G,  B, res) - mean;
+		}
 		R = delta * w / temp;
 		set_intensity(f, mean + R);
 		set_temp2(f, M2 + sumweight * delta * R);
@@ -352,6 +343,25 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
 }
 
 
+/* Correct intensity in pattern for scaling and Lorentz factors,
+ *  but not partiality nor polarisation */
+double correct_reflection_nopart(Reflection *refl, double osf, double Bfac,
+                                 double res)
+{
+	double corr = osf * exp(-Bfac*res*res);
+	return (get_intensity(refl) / corr) / get_lorentz(refl);
+}
+
+
+/* Correct intensity in pattern for scaling, partiality and Lorentz factors,
+ *  but not polarisation */
+double correct_reflection(Reflection *refl, double osf, double Bfac, double res)
+{
+	double Ipart = correct_reflection_nopart(refl, osf, Bfac, res);
+	return Ipart / get_partiality(refl);
+}
+
+
 double residual(Crystal *cr, const RefList *full, int free,
                 int *pn_used, const char *filename)
 {
@@ -368,7 +378,7 @@ double residual(Crystal *cr, const RefList *full, int free,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double w, corr, res;
+		double w, res;
 		signed int h, k, l;
 		Reflection *match;
 		double I_full;
@@ -384,8 +394,7 @@ double residual(Crystal *cr, const RefList *full, int free,
 
 		if ( get_redundancy(match) < 2 ) continue;
 
-		corr = get_lorentz(refl) / (G * exp(-B*res*res));
-		int1 = get_intensity(refl) * corr;
+		int1 = correct_reflection_nopart(refl, G, B, res);
 		int2 = get_partiality(refl)*I_full;
 		w = 1.0;
 
diff --git a/src/merge.h b/src/merge.h
index 476f8755..2c2cfeda 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -47,6 +47,12 @@ extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
                                   int min_meas, double push_res, int use_weak,
                                   int ln_merge);
 
+extern double correct_reflection_nopart(Reflection *refl, double osf,
+                                        double Bfac, double res);
+
+extern double correct_reflection(Reflection *refl, double osf, double Bfac,
+                                 double res);
+
 extern double residual(Crystal *cr, const RefList *full, int free,
                        int *pn_used, const char *filename);
 
diff --git a/src/partialator.c b/src/partialator.c
index 10befba8..8cf1dad9 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -638,7 +638,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 	{
 		signed int h, k, l;
 		double pobs, pcalc;
-		double res, corr, Ipart;
+		double res, Ipart;
 		Reflection *match;
 
 		if ( !get_flag(refl) ) continue;  /* Not free-flagged */
@@ -660,8 +660,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 		pcalc = get_partiality(refl);
 
 		/* Observed partiality */
-		corr = G * exp(B*res*res) * get_lorentz(refl);
-		Ipart = get_intensity(refl) * corr;
+		Ipart = correct_reflection_nopart(refl, G, B, res);
 		pobs = Ipart / get_intensity(match);
 
 		fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n",
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 0a799b9e..b1af92db 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -559,7 +559,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double corr, Ipart, Ifull, pobs, pcalc;
+		double Ipart, Ifull, pobs, pcalc;
 		double res;
 		signed int h, k, l;
 		Reflection *match;
@@ -573,8 +573,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
 		/* Don't calculate pobs if reference reflection is weak */
 		if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
 
-		corr = G * exp(B*res*res) * get_lorentz(refl);
-		Ipart = get_intensity(refl) * corr;
+		Ipart = correct_reflection_nopart(refl, G, B, res);
 		Ifull = get_intensity(match);
 		pobs = Ipart / Ifull;
 		pcalc = get_partiality(refl);
diff --git a/src/rejection.c b/src/rejection.c
index 6ec59151..5d1041fa 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -143,11 +143,7 @@ static int calculate_refl_mean_var(RefList *full)
 
 			G = crystal_get_osf(c->contrib_crystals[j]);
 			B = crystal_get_Bfac(c->contrib_crystals[j]);
-			Ii = get_intensity(c->contribs[j]);
-
-			/* Total (multiplicative) correction factor */
-			Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j])
-			          / get_partiality(c->contribs[j]);
+			Ii = correct_reflection(c->contribs[j], G, B, res);
 
 			Ex += Ii - K;
 			Ex2 += (Ii - K) * (Ii - K);
@@ -232,16 +228,13 @@ static double calculate_cchalf(RefList *template, RefList *full,
 		while ( exrefl != NULL ) {
 
 			double G, B;
-			double Ii = get_intensity(exrefl);
 
 			G = crystal_get_osf(exclude);
 			B = crystal_get_Bfac(exclude);
 
 			if ( get_partiality(exrefl) > MIN_PART_MERGE ) {
 
-				/* Total (multiplicative) correction factor */
-				Ii *= 1.0/G * exp(B*res*res) * get_lorentz(exrefl)
-				          / get_partiality(exrefl);
+				double Ii = correct_reflection(exrefl, G, B, res);
 
 				/* Remove contribution of this reflection */
 				Ex -= Ii - K;