New partiality model from Ginn et al.

4 files changed, 15 insertions, 52 deletions

diff --git a/src/partial_sim.c b/src/partial_sim.c
index 4a348742..b2581096 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -386,7 +386,7 @@ static void run_job(void *vwargs, int cookie)
 
 	reflections = predict_to_res(cr, largest_q(&wargs->image));
 	crystal_set_reflections(cr, reflections);
-	calculate_partialities(cr, PMODEL_SCSPHERE);
+	calculate_partialities(cr, PMODEL_XSPHERE);
 
 	for ( i=0; i<NBINS; i++ ) {
 		wargs->n_ref[i] = 0;
diff --git a/src/partialator.c b/src/partialator.c
index 428b0d68..5ad106c7 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -770,7 +770,7 @@ int main(int argc, char *argv[])
 	Stream *st;
 	Crystal **crystals;
 	char *pmodel_str = NULL;
-	PartialityModel pmodel = PMODEL_SCSPHERE;
+	PartialityModel pmodel = PMODEL_XSPHERE;
 	int min_measurements = 2;
 	char *rval;
 	int polarisation = 1;
@@ -975,10 +975,8 @@ int main(int argc, char *argv[])
 	if ( pmodel_str != NULL ) {
 		if ( strcmp(pmodel_str, "unity") == 0 ) {
 			pmodel = PMODEL_UNITY;
-		} else if ( strcmp(pmodel_str, "scgaussian") == 0 ) {
-			pmodel = PMODEL_SCGAUSSIAN;
-		} else if ( strcmp(pmodel_str, "scsphere") == 0 ) {
-			pmodel = PMODEL_SCSPHERE;
+		} else if ( strcmp(pmodel_str, "xsphere") == 0 ) {
+			pmodel = PMODEL_XSPHERE;
 		} else if ( strcmp(pmodel_str, "random") == 0 ) {
 			pmodel = PMODEL_RANDOM;
 		} else {
diff --git a/src/post-refinement.c b/src/post-refinement.c
index bf30d299..74075503 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -85,20 +85,6 @@ const char *str_prflag(enum prflag flag)
 }
 
 
-/* Returns dp(gauss)/dr at "r" */
-static double gaussian_fraction_gradient(double r, double R)
-{
-	const double ng = 2.6;
-	const double sig = R/ng;
-
-	/* If the Ewald sphere isn't within the profile, the gradient is zero */
-	if ( r < -R ) return 0.0;
-	if ( r > +R ) return 0.0;
-
-	return exp(-pow(r/sig, 2.0)/2.0) / (sig*sqrt(2.0*M_PI));
-}
-
-
 /* Returns dp(sph)/dr at "r" */
 static double sphere_fraction_gradient(double r, double pr)
 {
@@ -130,14 +116,10 @@ static double partiality_gradient(double r, double pr,
 		case PMODEL_UNITY:
 		return 0.0;
 
-		case PMODEL_SCSPHERE:
+		case PMODEL_XSPHERE:
 		A = sphere_fraction_gradient(r, pr)/D;
 		return 4.0*pr*A/3.0;
 
-		case PMODEL_SCGAUSSIAN:
-		A = gaussian_fraction_gradient(r, pr)/D;
-		return 4.0*pr*A/3.0;
-
 	}
 }
 
@@ -152,19 +134,6 @@ static double sphere_fraction_rgradient(double r, double R)
 }
 
 
-static double gaussian_fraction_rgradient(double r, double R)
-{
-	const double ng = 2.6;
-	const double sig = R/ng;
-
-	/* If the Ewald sphere isn't within the profile, the gradient is zero */
-	if ( r < -R ) return 0.0;
-	if ( r > +R ) return 0.0;
-
-	return -(ng*r/(sqrt(2.0*M_PI)*R*R))*exp(-r*r/(2.0*sig*sig));
-}
-
-
 static double volume_fraction_rgradient(double r, double pr,
                                        PartialityModel pmodel)
 {
@@ -173,12 +142,9 @@ static double volume_fraction_rgradient(double r, double pr,
 		case PMODEL_UNITY :
 		return 1.0;
 
-		case PMODEL_SCSPHERE :
+		case PMODEL_XSPHERE :
 		return sphere_fraction_rgradient(r, pr);
 
-		case PMODEL_SCGAUSSIAN :
-		return gaussian_fraction_rgradient(r, pr);
-
 		default :
 		ERROR("No pmodel in volume_fraction_rgradient!\n");
 		return 1.0;
@@ -194,12 +160,9 @@ static double volume_fraction(double rlow, double rhigh, double pr,
 		case PMODEL_UNITY :
 		return 1.0;
 
-		case PMODEL_SCSPHERE :
+		case PMODEL_XSPHERE :
 		return sphere_fraction(rlow, rhigh, pr);
 
-		case PMODEL_SCGAUSSIAN :
-		return gaussian_fraction(rlow, rhigh, pr);
-
 		default :
 		ERROR("No pmodel in volume_fraction!\n");
 		return 1.0;
@@ -212,12 +175,14 @@ static double volume_fraction(double rlow, double rhigh, double pr,
 double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 {
 	double glow, ghigh;
-	double rlow, rhigh, p;
+	double rlow, rhigh;
 	struct image *image = crystal_get_image(cr);
 	double R = crystal_get_profile_radius(cr);
 	double gr;
 
-	get_partial(refl, &rlow, &rhigh, &p);
+	/* FIXME ! */
+	rlow = 0.0;
+	rhigh = 0.0;
 
 	if ( k == GPARAM_R ) {
 
@@ -607,7 +572,7 @@ static void write_residual_graph(Crystal *cr, const RefList *full)
 		                          bsx, bsy, bsz,
 		                          csx, csy, csz);
 		update_predictions(cr);
-		calculate_partialities(cr, PMODEL_SCSPHERE);
+		calculate_partialities(cr, PMODEL_XSPHERE);
 		res = residual(cr, full, 0, &n, NULL);
 		fprintf(fh, "%i %e %e %i\n", i, asx, res, n);
 	}
@@ -616,7 +581,7 @@ static void write_residual_graph(Crystal *cr, const RefList *full)
 	                          bsx, bsy, bsz,
 	                          csx, csy, csz);
 	update_predictions(cr);
-	calculate_partialities(cr, PMODEL_SCSPHERE);
+	calculate_partialities(cr, PMODEL_XSPHERE);
 	fclose(fh);
 }
 
diff --git a/src/process_image.c b/src/process_image.c
index 87089289..dbd4f427 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -246,7 +246,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	if ( iargs->fix_bandwidth >= 0.0 ) {
 		image.bw = iargs->fix_bandwidth;
 	} else {
-		image.bw = 0.00000001;
+		image.bw = 0.0013;
 	}
 
 	if ( image_feature_count(image.features) < iargs->min_peaks ) {
@@ -301,7 +301,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	/* Integrate! */
 	time_accounts_set(taccs, TACC_INTEGRATION);
 	sb_shared->pings[cookie]++;
-	integrate_all_5(&image, iargs->int_meth, PMODEL_SCSPHERE,
+	integrate_all_5(&image, iargs->int_meth, PMODEL_XSPHERE,
 	                iargs->push_res,
 	                iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
 	                iargs->int_diag, iargs->int_diag_h,