WIP on gradients for scsphere

2 files changed, 33 insertions, 14 deletions

diff --git a/src/post-refinement.c b/src/post-refinement.c
index 99ddf6b1..4c16f2d8 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -68,9 +68,14 @@ static double dpdq(double r, double profile_radius)
 
 /* Returns dp/dr at "r" */
 static double partiality_gradient(double r, double profile_radius,
-                                  PartialityModel pmodel)
+                                  PartialityModel pmodel,
+                                  double rlow, double rhigh, double p)
 {
-	double dqdr;  /* dq/dr */
+	double dqdr;    /* dq/dr */
+	double dpsphdr; /* dpsph / dr */
+	double A, D;
+
+	D = rlow - rhigh;
 
 	switch ( pmodel ) {
 
@@ -90,6 +95,12 @@ static double partiality_gradient(double r, double profile_radius,
 		case PMODEL_THIN:
 		return -(2.0*r)/(profile_radius*profile_radius);
 
+		case PMODEL_SCSPHERE:
+		dqdr = 1.0 / (2.0*profile_radius);
+		dpsphdr = dpdq(r, profile_radius) * dqdr;
+		A = (dpsphdr/D) - p*pow(rlow-rhigh, -2.0);
+		return 4.0*profile_radius*A/3.0;
+
 	}
 }
 
@@ -178,12 +189,12 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 
 	/* Calculate the gradient of partiality wrt excitation error. */
 	if ( clamp_low == 0 ) {
-		glow = partiality_gradient(rlow, r, pmodel);
+		glow = partiality_gradient(rlow, r, pmodel, rlow, rhigh, p);
 	} else {
 		glow = 0.0;
 	}
 	if ( clamp_high == 0 ) {
-		ghigh = partiality_gradient(rhigh, r, pmodel);
+		ghigh = partiality_gradient(rhigh, r, pmodel, rlow, rhigh, p);
 	} else {
 		ghigh = 0.0;
 	}
@@ -215,14 +226,15 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 
 			default:
 			case PMODEL_UNITY:
+			case PMODEL_THIN:
 			return 0.0;
 
 			case PMODEL_GAUSSIAN:
 			case PMODEL_SPHERE:
 			return (ds*glow + ds*ghigh) / 2.0;
 
-			case PMODEL_THIN:
-			return 0.0;
+			case PMODEL_SCSPHERE:
+			return 0.0; /* FIXME */
 		}
 	}
 
@@ -273,6 +285,9 @@ double l_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 	signed int hs, ks, ls;
 	double L;
 
+	/* L has a non-zero gradient only for div in sphere or gaussian model */
+	if ( (pmodel != PMODEL_SPHERE)
+	  && (pmodel != PMODEL_GAUSSIAN) ) return 0.0;
 	if ( k != REF_DIV ) return 0.0;
 
 	get_symmetric_indices(refl, &hs, &ks, &ls);
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index b88380eb..82ed2b02 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -45,7 +45,8 @@
 
 
 static void scan_partialities(RefList *reflections, RefList *compare,
-                              int *valid, long double *vals[3], int idx)
+                              int *valid, long double *vals[3], int idx,
+                              PartialityModel pmodel)
 {
 	int i;
 	Reflection *refl;
@@ -70,7 +71,7 @@ static void scan_partialities(RefList *reflections, RefList *compare,
 		}
 
 		get_partial(refl2, &r1, &r2, &p, &clamp_low, &clamp_high);
-		if ( clamp_low && clamp_high ) {
+		if ( clamp_low && clamp_high && (pmodel != PMODEL_SCSPHERE) ) {
 			if ( !within_tolerance(p, 1.0, 0.001) ) {
 
 				signed int h, k, l;
@@ -190,7 +191,7 @@ static void calc_either_side(Crystal *cr, double incr_val,
 		cr_new = new_shifted_crystal(cr, refine, -incr_val);
 		compare = find_intersections(image, cr_new, pmodel);
 		scan_partialities(crystal_get_reflections(cr), compare, valid,
-		                  vals, 0);
+		                  vals, 0, pmodel);
 		cell_free(crystal_get_cell(cr_new));
 		crystal_free(cr_new);
 		reflist_free(compare);
@@ -198,7 +199,7 @@ static void calc_either_side(Crystal *cr, double incr_val,
 		cr_new = new_shifted_crystal(cr, refine, +incr_val);
 		compare = find_intersections(image, cr_new, pmodel);
 		scan_partialities(crystal_get_reflections(cr), compare, valid,
-		                  vals, 2);
+		                  vals, 2, pmodel);
 		cell_free(crystal_get_cell(cr_new));
 		crystal_free(cr_new);
 		reflist_free(compare);
@@ -212,14 +213,14 @@ static void calc_either_side(Crystal *cr, double incr_val,
 		shift_parameter(&im_moved, refine, -incr_val);
 		compare = find_intersections(&im_moved, cr, pmodel);
 		scan_partialities(crystal_get_reflections(cr), compare,
-		                  valid, vals, 0);
+		                  valid, vals, 0, pmodel);
 		reflist_free(compare);
 
 		im_moved = *image;
 		shift_parameter(&im_moved, refine, +incr_val);
 		compare = find_intersections(&im_moved, cr, pmodel);
 		scan_partialities(crystal_get_reflections(cr), compare,
-		                  valid, vals, 2);
+		                  valid, vals, 2, pmodel);
 		reflist_free(compare);
 
 	}
@@ -271,7 +272,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 	}
 	for ( i=0; i<nref; i++ ) valid[i] = 1;
 
-	scan_partialities(reflections, reflections, valid, vals, 1);
+	scan_partialities(reflections, reflections, valid, vals, 1, pmodel);
 
 	calc_either_side(cr, incr_val, valid, vals, refine, pmodel);
 
@@ -450,7 +451,7 @@ int main(int argc, char *argv[])
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
-	for ( i=0; i<3; i++ ) {
+	for ( i=0; i<4; i++ ) {
 
 		UnitCell *rot;
 		double val;
@@ -467,6 +468,9 @@ int main(int argc, char *argv[])
 		} else if ( i == 2 ) {
 			pmodel = PMODEL_THIN;
 			STATUS("Testing Thin Ewald Sphere model:\n");
+		} else if ( i == 3 ) {
+			pmodel = PMODEL_SCSPHERE;
+			STATUS("Testing SCSphere model:\n");
 		} else {
 			ERROR("WTF?\n");
 			return 1;