Use linear scale when scaling against reference

1 files changed, 8 insertions, 86 deletions

diff --git a/src/post-refinement.c b/src/post-refinement.c
index 165080b4..5624b15d 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -44,6 +44,7 @@
 #include "cell.h"
 #include "cell-utils.h"
 #include "reflist-utils.h"
+#include "scaling.h"
 
 
 struct prdata
@@ -79,104 +80,23 @@ const char *str_prflag(enum prflag flag)
 }
 
 
-static int linear_scale(RefList *list1, const RefList *list2, double *G)
-{
-	Reflection *refl1;
-	Reflection *refl2;
-	RefListIterator *iter;
-	int max_n = 256;
-	int n = 0;
-	double *x;
-	double *y;
-	double *w;
-	int r;
-	double cov11;
-	double sumsq;
-
-	x = malloc(max_n*sizeof(double));
-	w = malloc(max_n*sizeof(double));
-	y = malloc(max_n*sizeof(double));
-	if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
-		ERROR("Failed to allocate memory for scaling.\n");
-		return 1;
-	}
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		signed int h, k, l;
-		double Ih1, Ih2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;
-
-		Ih1 = get_intensity(refl1);
-		Ih2 = get_intensity(refl2);
-
-		if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue;
-		if ( isnan(Ih1) || isinf(Ih1) ) continue;
-		if ( isnan(Ih2) || isinf(Ih2) ) continue;
-
-		if ( n == max_n ) {
-			max_n *= 2;
-			x = realloc(x, max_n*sizeof(double));
-			y = realloc(y, max_n*sizeof(double));
-			w = realloc(w, max_n*sizeof(double));
-			if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
-				ERROR("Failed to allocate memory for scaling.\n");
-				return 1;
-			}
-		}
-
-		x[n] = Ih2;
-		y[n] = Ih1;
-		w[n] = get_partiality(refl1);
-		n++;
-
-	}
-
-	if ( n < 2 ) {
-		ERROR("Not enough reflections for scaling\n");
-		return 1;
-	}
-
-	r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq);
-
-	if ( r ) {
-		ERROR("Scaling failed.\n");
-		return 1;
-	}
-
-	free(x);
-	free(y);
-	free(w);
-
-	return 0;
-}
-
-
-
 double residual(Crystal *cr, const RefList *full, int free,
                 int *pn_used, const char *filename)
 {
-	double G;
 	Reflection *refl;
 	RefListIterator *iter;
 	int n_used = 0;
 	double num = 0.0;
 	double den = 0.0;
-
-	if ( linear_scale(crystal_get_reflections(cr), full, &G) ) {
-		return INFINITY;
-	}
+	double G = crystal_get_osf(cr);
+	double B = crystal_get_Bfac(cr);
+	UnitCell *cell = crystal_get_cell(cr);
 
 	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double p, w;
+		double p, w, corr, res;
 		signed int h, k, l;
 		Reflection *match;
 		double I_full;
@@ -185,6 +105,7 @@ double residual(Crystal *cr, const RefList *full, int free,
 		if ( free && !get_flag(refl) ) continue;
 
 		get_indices(refl, &h, &k, &l);
+		res = resolution(cell, h, k, l);
 		match = find_refl(full, h, k, l);
 		if ( match == NULL ) continue;
 		I_full = get_intensity(match);
@@ -194,7 +115,8 @@ double residual(Crystal *cr, const RefList *full, int free,
 		p = get_partiality(refl);
 		//if ( p < 0.2 ) continue;
 
-		int1 = get_intensity(refl) / G;
+		corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+		int1 = get_intensity(refl) * corr;
 		int2 = p*I_full;
 		w = 1.0;