partialator: Account for Lorentz factor

4 files changed, 58 insertions, 15 deletions

diff --git a/src/post-refinement.c b/src/post-refinement.c
index 9e4649a2..cbb4e5b5 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -87,8 +87,9 @@ static double partiality_rgradient(double r, double profile_radius)
 }
 
 
-/* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
+/* Return the gradient of partiality wrt parameter 'k' given the current status
+ * of 'image'. */
+double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 {
 	double ds, azi;
 	double glow, ghigh;
@@ -200,6 +201,46 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 }
 
 
+/* Return the gradient of Lorentz factor wrt parameter 'k' given the current
+ * status of 'image'. */
+double l_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
+{
+	double ds;
+	signed int hs, ks, ls;
+
+	switch ( k ) {
+
+		/* Cell parameters do not affect Lorentz factor */
+		case REF_ASX :
+		case REF_BSX :
+		case REF_CSX :
+		case REF_ASY :
+		case REF_BSY :
+		case REF_CSY :
+		case REF_ASZ :
+		case REF_BSZ :
+		case REF_CSZ :
+		return 0.0;
+
+		/* Nor does change of radius */
+		case REF_R :
+		return 0.0;
+
+		default:
+		break;
+
+	}
+
+	assert(k == REF_DIV);
+
+	get_symmetric_indices(refl, &hs, &ks, &ls);
+
+	ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
+
+	return 0.0;  /* FIXME! */
+}
+
+
 static void apply_cell_shift(UnitCell *cell, int k, double shift)
 {
 	double asx, asy, asz;
@@ -383,7 +424,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 		double w;
 		double I_partial;
 		int k;
-		double p;
+		double p, l;
 		Reflection *match;
 		double gradients[NUM_PARAMS];
 		const double osf = crystal_get_osf(cr);
@@ -398,21 +439,23 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 			      " is it marked as refinable?\n", ha, ka, la);
 			continue;
 		}
-		I_full = osf * get_intensity(match);
+		I_full = get_intensity(match);
 
 		/* Actual measurement of this reflection from this pattern? */
-		I_partial = get_intensity(refl);
+		I_partial = osf * get_intensity(refl);
 		p = get_partiality(refl);
+		l = get_lorentz(refl);
 
 		/* Calculate the weight for this reflection */
 		w =  pow(get_esd_intensity(refl), 2.0);
-		w += p * I_full * pow(get_esd_intensity(match), 2.0);
+		w += l * p * I_full * pow(get_esd_intensity(match), 2.0);
 		w = pow(w, -1.0);
 
 		/* Calculate all gradients for this reflection */
 		for ( k=0; k<NUM_PARAMS; k++ ) {
 			double gr;
-			gr = gradient(cr, k, refl, pmodel);
+			gr  = p_gradient(cr, k, refl, pmodel) * l;
+			gr += l_gradient(cr, k, refl, pmodel) * p;
 			gradients[k] = gr;
 		}
 
@@ -429,7 +472,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 				if ( g > k ) continue;
 
 				M_c = gradients[g] * gradients[k];
-				M_c *= w * pow(osf * I_full, 2.0);
+				M_c *= w * pow(I_full, 2.0);
 
 				M_curr = gsl_matrix_get(M, k, g);
 				gsl_matrix_set(M, k, g, M_curr + M_c);
@@ -437,7 +480,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 
 			}
 
-			delta_I = I_partial - (p * I_full);
+			delta_I = I_partial - (l * p * I_full);
 			v_c = w * delta_I * I_full * gradients[k];
 			v_curr = gsl_vector_get(v, k);
 			gsl_vector_set(v, k, v_curr + v_c);
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 7b822938..f565a880 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -63,8 +63,8 @@ enum {
 extern void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel);
 
 /* Exported so it can be poked by tests/pr_gradient_check */
-extern double gradient(Crystal *cr, int k, Reflection *refl,
-                       PartialityModel pmodel);
+extern double p_gradient(Crystal *cr, int k, Reflection *refl,
+                         PartialityModel pmodel);
 
 
 #endif	/* POST_REFINEMENT_H */
diff --git a/src/scaling-report.c b/src/scaling-report.c
index a02a1796..1d7904a6 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -251,13 +251,13 @@ static void partiality_graph(cairo_t *cr, Crystal **crystals, int n,
 			if ( f == NULL ) continue;
 			if ( get_redundancy(f) < 2 ) continue;
 
-			Ipart = get_intensity(refl);
-			Ifull = crystal_get_osf(cryst) * get_intensity(f);
+			Ipart = crystal_get_osf(cryst) * get_intensity(refl);
+			Ifull = get_intensity(f);
 
 			//if ( Ifull < 10 ) continue;  /* FIXME: Ugh */
 
 			pobs = Ipart/Ifull;
-			pcalc = get_partiality(refl);
+			pcalc = get_lorentz(refl) * get_partiality(refl);
 
 			//STATUS("%4i %4i %4i : %9.6f %9.6f %e %e %e\n",
 			//       h, k, l, pobs, pcalc,
diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c
index 68a625a7..4600db9b 100644
--- a/tests/pr_gradient_check.c
+++ b/tests/pr_gradient_check.c
@@ -309,7 +309,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 			grad = (grad1 + grad2) / 2.0;
 			i++;
 
-			cgrad = gradient(cr, refine, refl, pmodel);
+			cgrad = p_gradient(cr, refine, refl, pmodel);
 
 			get_partial(refl, &r1, &r2, &p, &cl, &ch);