Remove "sphere", "thin" and "gaussian" partiality models, add "scgaussian"

5 files changed, 25 insertions, 70 deletions

diff --git a/src/partial_sim.c b/src/partial_sim.c
index 10f5e58d..77b5a73b 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -357,7 +357,7 @@ static void run_job(void *vwargs, int cookie)
 		snprintf(wargs->image.filename, 255, "dummy.h5");
 	}
 
-	reflections = find_intersections(&wargs->image, cr, PMODEL_SPHERE);
+	reflections = find_intersections(&wargs->image, cr, PMODEL_SCSPHERE);
 	crystal_set_reflections(cr, reflections);
 
 	for ( i=0; i<NBINS; i++ ) {
diff --git a/src/partialator.c b/src/partialator.c
index 0d532997..bdee29c5 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -219,7 +219,7 @@ int main(int argc, char *argv[])
 	Stream *st;
 	Crystal **crystals;
 	char *pmodel_str = NULL;
-	PartialityModel pmodel = PMODEL_SPHERE;
+	PartialityModel pmodel = PMODEL_SCSPHERE;
 	int min_measurements = 2;
 	char *rval;
 	struct srdata srdata;
@@ -342,14 +342,10 @@ int main(int argc, char *argv[])
 	free(sym_str);
 
 	if ( pmodel_str != NULL ) {
-		if ( strcmp(pmodel_str, "sphere") == 0 ) {
-			pmodel = PMODEL_SPHERE;
-		} else if ( strcmp(pmodel_str, "unity") == 0 ) {
+		if ( strcmp(pmodel_str, "unity") == 0 ) {
 			pmodel = PMODEL_UNITY;
-		} else if ( strcmp(pmodel_str, "gaussian") == 0 ) {
-			pmodel = PMODEL_GAUSSIAN;
-		} else if ( strcmp(pmodel_str, "thin") == 0 ) {
-			pmodel = PMODEL_THIN;
+		} else if ( strcmp(pmodel_str, "scgaussian") == 0 ) {
+			pmodel = PMODEL_SCGAUSSIAN;
 		} else if ( strcmp(pmodel_str, "scsphere") == 0 ) {
 			pmodel = PMODEL_SCSPHERE;
 		} else {
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 4c16f2d8..e971f96a 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -83,24 +83,17 @@ static double partiality_gradient(double r, double profile_radius,
 		case PMODEL_UNITY:
 		return 0.0;
 
-		case PMODEL_SPHERE:
-		dqdr = 1.0 / (2.0*profile_radius);
-		return dpdq(r, profile_radius) * dqdr;
-
-		case PMODEL_GAUSSIAN:
-		/* FIXME: Get a proper gradient */
-		dqdr = 1.0 / (2.0*profile_radius);
-		return dpdq(r, profile_radius) * dqdr;
-
-		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;
 
+		case PMODEL_SCGAUSSIAN:
+		/* FIXME: Get a proper gradient */
+		dqdr = 1.0 / (2.0*profile_radius);
+		return dpdq(r, profile_radius) * dqdr;
+
 	}
 }
 
@@ -117,18 +110,16 @@ static double partiality_rgradient(double r, double profile_radius,
 		case PMODEL_UNITY:
 		return 0.0;
 
-		case PMODEL_SPHERE:
+		case PMODEL_SCSPHERE:
+		/* FIXME: wrong (?) */
 		dqdrad = -0.5 * r / (profile_radius * profile_radius);
 		return dpdq(r, profile_radius) * dqdrad;
 
-		case PMODEL_GAUSSIAN:
+		case PMODEL_SCGAUSSIAN:
 		/* FIXME: Get a proper gradient */
 		dqdrad = -0.5 * r / (profile_radius * profile_radius);
 		return dpdq(r, profile_radius) * dqdrad;
 
-		case PMODEL_THIN:
-		return 2.0*r*r*pow(profile_radius, -3.0);
-
 	}
 }
 
@@ -206,18 +197,16 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 			case PMODEL_UNITY:
 			return 0.0;
 
-			case PMODEL_GAUSSIAN:
+			case PMODEL_SCSPHERE:
 			gr  = partiality_rgradient(rlow, r, pmodel);
 			gr -= partiality_rgradient(rhigh, r, pmodel);
 			return gr;
 
-			case PMODEL_SPHERE:
+			case PMODEL_SCGAUSSIAN:
 			gr  = partiality_rgradient(rlow, r, pmodel);
 			gr -= partiality_rgradient(rhigh, r, pmodel);
 			return gr;
 
-			case PMODEL_THIN:
-			return 2.0*rlow*rlow/(r*r*r);
 		}
 	}
 
@@ -226,15 +215,14 @@ 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_SCSPHERE:
 			return 0.0; /* FIXME */
+
+			case PMODEL_SCGAUSSIAN:
+			return (ds*glow + ds*ghigh) / 2.0; /* FIXME: Wrong */
+
 		}
 	}
 
@@ -277,28 +265,6 @@ double p_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;
-	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);
-
-	ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
-
-	L = get_lorentz(refl);
-	return -ds*L*L / LORENTZ_SCALE;
-}
-
-
 static void apply_cell_shift(UnitCell *cell, int k, double shift)
 {
 	double asx, asy, asz;
@@ -564,10 +530,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 
 		/* Calculate all gradients for this reflection */
 		for ( k=0; k<NUM_PARAMS; k++ ) {
-			double gr;
-			gr  = p_gradient(cr, k, refl, pmodel) * l;
-			gr += l_gradient(cr, k, refl, pmodel) * p;
-			gradients[k] = gr;
+			gradients[k] = p_gradient(cr, k, refl, pmodel) * l;
 		}
 
 		for ( k=0; k<NUM_PARAMS; k++ ) {
diff --git a/src/post-refinement.h b/src/post-refinement.h
index e419c51d..0c564690 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,11 +3,11 @@
  *
  * Post refinement
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -71,12 +71,8 @@ struct prdata
 extern struct prdata pr_refine(Crystal *cr, const RefList *full,
                                PartialityModel pmodel);
 
-/* Exported so it can be poked by tests/pr_gradient_check */
+/* Exported so it can be poked by tests/pr_p_gradient_check */
 extern double p_gradient(Crystal *cr, int k, Reflection *refl,
                          PartialityModel pmodel);
 
-
-extern double l_gradient(Crystal *cr, int k, Reflection *refl,
-                         PartialityModel pmodel);
-
 #endif	/* POST_REFINEMENT_H */
diff --git a/src/process_image.c b/src/process_image.c
index 83973e50..1d1dac82 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -178,7 +178,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 
 	/* Integrate all the crystals at once - need all the crystals so that
 	 * overlaps can be detected. */
-	integrate_all_4(&image, iargs->int_meth, PMODEL_SPHERE, iargs->push_res,
+	integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE, iargs->push_res,
 	                iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
 	                iargs->int_diag, iargs->int_diag_h,
 	                iargs->int_diag_k, iargs->int_diag_l, results_pipe);