Polarisation correction for extracted intensities

4 files changed, 29 insertions, 10 deletions

diff --git a/src/indexamajig.c b/src/indexamajig.c
index bd3b53e2..0f7154e8 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -57,6 +57,7 @@ struct process_args
 	int config_gpu;
 	int config_simulate;
 	int config_nomatch;
+	int config_unpolar;
 	IndexingMethod indm;
 	const double *intensities;
 	const unsigned int *counts;
@@ -118,6 +119,7 @@ static void show_help(const char *s)
 "      --no-match          Don't attempt to match the indexed cell to the\n"
 "                           model, just proceed with the one generated by the\n"
 "                           auto-indexing procedure.\n"
+"      --unpolarized       Don't correct for the polarisation of the X-rays.\n"
 "\n\nOptions for greater performance or verbosity:\n\n"
 "      --verbose           Be verbose about indexing.\n"
 "      --gpu               Use the GPU to speed up the simulation.\n"
@@ -237,6 +239,7 @@ static void *process_image(void *pargsv)
 	int config_gpu = pargs->config_gpu;
 	int config_simulate = pargs->config_simulate;
 	int config_nomatch = pargs->config_nomatch;
+	int config_unpolar = pargs->config_unpolar;
 	IndexingMethod indm = pargs->indm;
 	const double *intensities = pargs->intensities;
 	const unsigned int *counts = pargs->counts;
@@ -321,7 +324,7 @@ static void *process_image(void *pargsv)
 		/* Use original data (temporarily) */
 		simage->data = data_for_measurement;
 		output_intensities(simage, image.indexed_cell,
-		                   pargs->output_mutex);
+		                   pargs->output_mutex, config_unpolar);
 		simage->data = NULL;
 	}
 
@@ -382,6 +385,7 @@ int main(int argc, char *argv[])
 	int config_gpu = 0;
 	int config_verbose = 0;
 	int config_alternate = 0;
+	int config_unpolar = 0;
 	IndexingMethod indm;
 	char *indm_str = NULL;
 	UnitCell *cell;
@@ -417,6 +421,7 @@ int main(int argc, char *argv[])
 		{"intensities",        1, NULL,               'q'},
 		{"pdb",                1, NULL,               'p'},
 		{"prefix",             1, NULL,               'x'},
+		{"unpolarized",        0, &config_unpolar,     1},
 		{0, 0, NULL, 0}
 	};
 
@@ -570,6 +575,7 @@ int main(int argc, char *argv[])
 		pargs->config_gpu = config_gpu;
 		pargs->config_simulate = config_simulate;
 		pargs->config_nomatch = config_nomatch;
+		pargs->config_unpolar = config_unpolar;
 		pargs->cell = cell;
 		pargs->indm = indm;
 		pargs->intensities = intensities;
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index a0d81a9d..e09ab72a 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -346,7 +346,7 @@ int main(int argc, char *argv[])
 		record_image(&image, !config_nonoise);
 
 		if ( config_nearbragg ) {
-			output_intensities(&image, cell, NULL);
+			output_intensities(&image, cell, NULL, 1);
 		}
 
 		if ( config_powder ) {
diff --git a/src/peaks.c b/src/peaks.c
index cf4eacf3..af17ba73 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -41,6 +41,9 @@
  * for the purposes of integration. */
 #define PEAK_REALLY_CLOSE (10.0)
 
+/* Degree of polarisation of X-ray beam */
+#define POL (1.0)
+
 struct reflhit {
 	signed int h;
 	signed int k;
@@ -140,7 +143,8 @@ static void cull_peaks(struct image *image)
 
 
 static void integrate_peak(struct image *image, int xp, int yp,
-                           float *xc, float *yc, float *intensity)
+                           float *xc, float *yc, float *intensity,
+                           int do_polar)
 {
 	signed int x, y;
 	const int lim = INTEGRATION_RADIUS * INTEGRATION_RADIUS;
@@ -154,6 +158,7 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		struct panel *p;
 		double val, sa, pix_area, Lsq, dsq, proj_area;
 		float tt;
+		double phi, pa, pb, pol;
 
 		/* Circular mask */
 		if ( x*x + y*y > lim ) continue;
@@ -178,7 +183,12 @@ static void integrate_peak(struct image *image, int xp, int yp,
 		/* Projected area of pixel divided by distance squared */
 		sa = 1.0e7 * proj_area / (dsq + Lsq);
 
-		val = image->data[(x+xp)+image->width*(y+yp)] / sa;
+		phi = atan2(y+yp, x+xp);
+		pa = pow(sin(phi)*sin(tt), 2.0);
+		pb = pow(cos(tt), 2.0);
+		pol = 1.0 - 2.0*POL*(1-pa) + POL*(1.0+pb);
+
+		val = image->data[(x+xp)+image->width*(y+yp)] / (sa*pol);
 
 		total += val;
 
@@ -305,8 +315,10 @@ void search_peaks(struct image *image)
 			continue;
 		}
 
-		/* Centroid peak and get better coordinates */
-		integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity);
+		/* Centroid peak and get better coordinates.
+		 * Don't bother doing polarisation correction, because the
+		 * intensity of this peak is only an estimate at this stage. */
+		integrate_peak(image, mask_x, mask_y, &fx, &fy, &intensity, 0);
 
 		/* It is possible for the centroid to fall outside the image */
 		if ( (fx < 0.0) || (fx > image->width)
@@ -370,7 +382,7 @@ void dump_peaks(struct image *image, pthread_mutex_t *mutex)
 
 
 void output_intensities(struct image *image, UnitCell *cell,
-                        pthread_mutex_t *mutex)
+                        pthread_mutex_t *mutex, int unpolar)
 {
 	int x, y;
 	double ax, ay, az;
@@ -485,13 +497,14 @@ void output_intensities(struct image *image, UnitCell *cell,
 			/* f->intensity was measured on the filtered pattern,
 			 * so instead re-integrate using old coordinates.
 			 * This will produce further revised coordinates. */
-			integrate_peak(image, f->x, f->y, &x, &y, &intensity);
+			integrate_peak(image, f->x, f->y, &x, &y, &intensity,
+			               !unpolar);
 			intensity = f->intensity;
 
 		} else {
 
 			integrate_peak(image, hits[i].x, hits[i].y,
-			               &x, &y, &intensity);
+			               &x, &y, &intensity, !unpolar);
 
 		}
 
diff --git a/src/peaks.h b/src/peaks.h
index a5c771b7..fa912455 100644
--- a/src/peaks.h
+++ b/src/peaks.h
@@ -22,6 +22,6 @@
 extern void search_peaks(struct image *image);
 extern void dump_peaks(struct image *image, pthread_mutex_t *mutex);
 extern void output_intensities(struct image *image, UnitCell *cell,
-                               pthread_mutex_t *mutex);
+                               pthread_mutex_t *mutex, int unpolar);
 
 #endif	/* PEAKS_H */