Remove solid angle correction

It's never correct when using "bucket" integration, and always correct
when using "pixel" integration, so don't give the option.

8 files changed, 33 insertions, 65 deletions

diff --git a/src/geometry.c b/src/geometry.c
index 821aff2f..9bb41d12 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -283,7 +283,7 @@ double integrate_all(struct image *image, struct cpeak *cpeaks, int n)
 		float x, y, intensity;
 
 		if ( integrate_peak(image, cpeaks[i].x, cpeaks[i].y, &x, &y,
-                                    &intensity, NULL, NULL, 0, 0, 0) ) continue;
+                                    &intensity, NULL, NULL, 0, 0) ) continue;
 
 		itot += intensity;
 	}
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 83eee4d1..774cfde4 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -62,7 +62,6 @@ struct static_index_args
 	int config_polar;
 	int config_sanity;
 	int config_satcorr;
-	int config_sa;
 	int config_closer;
 	float threshold;
 	float min_gradient;
@@ -182,8 +181,6 @@ static void show_help(const char *s)
 "     --unpolarized        Don't correct for the polarisation of the X-rays.\n"
 "     --no-sat-corr        Don't correct values of saturated peaks using a\n"
 "                           table included in the HDF5 file.\n"
-"     --no-sa              Don't correct for the differing solid angles of\n"
-"                           the pixels.\n"
 "     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
 "     --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
 "                           Default: 100,000.\n"
@@ -415,7 +412,7 @@ static void process_image(void *pp, int cookie)
 	if ( config_nearbragg ) {
 		output_intensities(&image, image.indexed_cell,
 		                   pargs->static_args.output_mutex,
-		                   config_polar, pargs->static_args.config_sa,
+		                   config_polar,
 		                   pargs->static_args.config_closer,
 		                   pargs->static_args.ofh, 0, 0.1);
 	}
@@ -525,7 +522,6 @@ int main(int argc, char *argv[])
 	int config_polar = 1;
 	int config_sanity = 0;
 	int config_satcorr = 1;
-	int config_sa = 1;
 	int config_checkprefix = 1;
 	int config_closer = 1;
 	float threshold = 800.0;
@@ -586,7 +582,6 @@ int main(int argc, char *argv[])
 		{"check-sanity",       0, &config_sanity,      1},
 		{"no-sat-corr",        0, &config_satcorr,     0},
 		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
-		{"no-sa",              0, &config_sa,          0},
 		{"threshold",          1, NULL,               't'},
 		{"min-gradient",       1, NULL,                4},
 		{"no-check-prefix",    0, &config_checkprefix, 0},
@@ -882,7 +877,6 @@ int main(int argc, char *argv[])
 	qargs.static_args.config_polar = config_polar;
 	qargs.static_args.config_sanity = config_sanity;
 	qargs.static_args.config_satcorr = config_satcorr;
-	qargs.static_args.config_sa = config_sa;
 	qargs.static_args.config_closer = config_closer;
 	qargs.static_args.cellr = cellr;
 	qargs.static_args.threshold = threshold;
diff --git a/src/partialator.c b/src/partialator.c
index 7829f933..beef9b6f 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -208,8 +208,7 @@ static void integrate_image(struct image *image, ReflItemList *obs)
 		 * pattern? */
 		/* FIXME: Coordinates aren't whole numbers */
 		if ( integrate_peak(image, spots[j].x, spots[j].y,
-		                    &xc, &yc, &i_partial, NULL, NULL,
-		                    1, 1, 0) ) {
+		                    &xc, &yc, &i_partial, NULL, NULL, 1, 0) ) {
 			spots[j].valid = 0;
 			continue;
 		}
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 6c48c6ec..4584bb3c 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -565,7 +565,7 @@ int main(int argc, char *argv[])
 
 		if ( config_nearbragg ) {
 			find_projected_peaks(&image, cell, 0, 0.1);
-			output_intensities(&image, cell, NULL, 0, 1, 0, stdout,
+			output_intensities(&image, cell, NULL, 0, 0, stdout,
 			                   0, 0.1);
 			free(image.cpeaks);
 		}
diff --git a/src/peaks.c b/src/peaks.c
index 1ad29a03..88ac50bf 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -182,7 +182,7 @@ static int cull_peaks(struct image *image)
 int integrate_peak(struct image *image, int xp, int yp,
                    float *xc, float *yc, float *intensity,
                    double *pbg, double *pmax,
-                   int do_polar, int do_sa, int centroid)
+                   int do_polar, int centroid)
 {
 	signed int x, y;
 	int lim, out_lim;
@@ -204,7 +204,7 @@ int integrate_peak(struct image *image, int xp, int yp,
 	for ( x=-out_lim; x<+out_lim; x++ ) {
 	for ( y=-out_lim; y<+out_lim; y++ ) {
 
-		double val, sa, pix_area, Lsq, dsq, proj_area;
+		double val;
 		float tt = 0.0;
 		double phi, pa, pb, pol;
 		uint16_t flags;
@@ -238,30 +238,9 @@ int integrate_peak(struct image *image, int xp, int yp,
 
 		if ( val > max ) max = val;
 
-		if ( do_sa ) {
-
-			/* Area of one pixel */
-			pix_area = pow(1.0/p->res, 2.0);
-			Lsq = pow(p->clen, 2.0);
-
-			/* Area of pixel as seen from crystal (approximate) */
-			tt = get_tt(image, x+xp, y+yp);
-			proj_area = pix_area * cos(tt);
-
-			/* Calculate distance from crystal to pixel */
-			dsq = pow(((double)(x+xp) - p->cx) / p->res, 2.0);
-			dsq += pow(((double)(y+yp) - p->cy) / p->res, 2.0);
-
-			/* Projected area of pixel divided by distance squared */
-			sa = 1.0e7 * proj_area / (dsq + Lsq);
-
-			val /= sa;
-
-		}
-
 		if ( do_polar ) {
 
-			if ( !do_sa ) tt = get_tt(image, x+xp, y+yp);
+			tt = get_tt(image, x+xp, y+yp);
 
 			phi = atan2(y+yp, x+xp);
 			pa = pow(sin(phi)*sin(tt), 2.0);
@@ -420,7 +399,7 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
 		 * Don't bother doing polarisation/SA correction, because the
 		 * intensity of this peak is only an estimate at this stage. */
 		r = integrate_peak(image, mask_x, mask_y,
-		                   &fx, &fy, &intensity, NULL, NULL, 0, 0, 1);
+		                   &fx, &fy, &intensity, NULL, NULL, 0, 1);
 		if ( r ) {
 			/* Bad region - don't detect peak */
 			nrej_bad++;
@@ -702,7 +681,7 @@ static void output_header(FILE *ofh, UnitCell *cell, struct image *image)
 
 
 void output_intensities(struct image *image, UnitCell *cell,
-                        pthread_mutex_t *mutex, int polar, int sa,
+                        pthread_mutex_t *mutex, int polar,
                         int use_closer, FILE *ofh,
                         int circular_domain, double domain_r)
 {
@@ -766,7 +745,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 				 * revised coordinates. */
 				r = integrate_peak(image, f->x, f->y, &x, &y,
 					           &intensity, &bg, &max,
-					           polar, sa, 1);
+					           polar, 1);
 				if ( r ) {
 					/* The original peak (which also went
 					 * through integrate_peak(), but with
@@ -788,7 +767,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 					           image->cpeaks[i].x,
 					           image->cpeaks[i].y,
 					           &x, &y, &intensity, &bg, &max,
-					           polar, sa, 1);
+					           polar, 1);
 				if ( r ) {
 					/* Plain old ordinary peak veto */
 					n_veto++;
@@ -809,7 +788,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 				           image->cpeaks[i].x,
 				           image->cpeaks[i].y,
 				           &x, &y, &intensity, &bg, &max, polar,
-				           sa, 0);
+				           0);
 			if ( r ) {
 				/* Plain old ordinary peak veto */
 				n_veto++;
@@ -867,7 +846,7 @@ void output_intensities(struct image *image, UnitCell *cell,
 
 
 void output_pixels(struct image *image, UnitCell *cell,
-                   pthread_mutex_t *mutex, int do_polar, int do_sa,
+                   pthread_mutex_t *mutex, int do_polar,
                    FILE *ofh, int circular_domain, double domain_r)
 {
 	int i;
@@ -962,30 +941,27 @@ void output_pixels(struct image *image, UnitCell *cell,
 
 			if ( p->no_index ) continue;
 
-			if ( do_sa ) {
-
-				/* Area of one pixel */
-				pix_area = pow(1.0/p->res, 2.0);
-				Lsq = pow(p->clen, 2.0);
-
-				/* Area of pixel as seen from crystal  */
-				tt = get_tt(image, x, y);
-				proj_area = pix_area * cos(tt);
+			/* Area of one pixel */
+			pix_area = pow(1.0/p->res, 2.0);
+			Lsq = pow(p->clen, 2.0);
 
-				/* Calculate distance from crystal to pixel */
-				dsq = pow(((double)x - p->cx) / p->res, 2.0);
-				dsq += pow(((double)y - p->cy) / p->res, 2.0);
+			/* Area of pixel as seen from crystal  */
+			tt = get_tt(image, x, y);
+			proj_area = pix_area * cos(tt);
 
-				/* Projected area of pixel / distance squared */
-				sa = 1.0e7 * proj_area / (dsq + Lsq);
+			/* Calculate distance from crystal to pixel */
+			dsq = pow(((double)x - p->cx) / p->res, 2.0);
+			dsq += pow(((double)y - p->cy) / p->res, 2.0);
 
-				val /= sa;
+			/* Projected area of pixel / distance squared */
+			sa = 1.0e7 * proj_area / (dsq + Lsq);
 
-			}
+			/* Solid angle correction is needed in this case */
+			val /= sa;
 
 			if ( do_polar ) {
 
-				if ( !do_sa ) tt = get_tt(image, x, y);
+				tt = get_tt(image, x, y);
 
 				phi = atan2(y, x);
 				pa = pow(sin(phi)*sin(tt), 2.0);
diff --git a/src/peaks.h b/src/peaks.h
index ba1eea30..8b506e42 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -24,12 +24,12 @@ extern void search_peaks(struct image *image, float threshold,
 extern void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex);
 
 extern void output_intensities(struct image *image, UnitCell *cell,
-                               pthread_mutex_t *mutex, int polar, int sa,
+                               pthread_mutex_t *mutex, int polar,
                                int use_closer, FILE *ofh, int circular_domain,
                                double domain_r);
 
 extern void output_pixels(struct image *image, UnitCell *cell,
-                          pthread_mutex_t *mutex, int do_polar, int do_sa,
+                          pthread_mutex_t *mutex, int do_polar,
                           FILE *ofh, int circular_domain, double domain_r);
 
 extern int peak_sanity_check(struct image *image, UnitCell *cell,
@@ -39,6 +39,6 @@ extern int find_projected_peaks(struct image *image, UnitCell *cell,
 extern int integrate_peak(struct image *image, int xp, int yp,
                           float *xc, float *yc, float *intensity,
                           double *pbg, double *pmax,
-                          int do_polar, int do_sa, int centroid);
+                          int do_polar, int centroid);
 
 #endif	/* PEAKS_H */
diff --git a/src/post-refinement.c b/src/post-refinement.c
index d835a2ff..26bef459 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -246,7 +246,7 @@ double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
 		/* FIXME: Coordinates aren't whole numbers */
 		if ( integrate_peak(image, spots[h].x, spots[h].y,
 		                    &xc, &yc, &I_partial, NULL, NULL,
-		                    1, 1, 0) ) {
+		                    1, 0) ) {
 			continue;
 		}
 
@@ -318,7 +318,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 		/* FIXME: Coordinates aren't whole numbers */
 		if ( integrate_peak(image, spots[h].x, spots[h].y,
 		                    &xc, &yc, &I_partial, NULL, NULL,
-		                    1, 1, 0) ) {
+		                    1, 0) ) {
 			continue;
 		}
 
diff --git a/src/reintegrate.c b/src/reintegrate.c
index b5a388a6..9d882d4b 100644
--- a/src/reintegrate.c
+++ b/src/reintegrate.c
@@ -140,8 +140,7 @@ static void process_image(void *pg, int cookie)
 
 		output_intensities(&image, apargs->cell,
 		                   pargs->output_mutex, pargs->config_polar,
-		                   pargs->config_sa, pargs->config_closer,
-		                   pargs->ofh, 0, 0.1);
+		                   pargs->config_closer, pargs->ofh, 0, 0.1);
 	}
 
 	free(image.data);