From e4a0b2c4f00ad396abd3c5c7e8cb55e2f39bc2cc Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 25 May 2010 17:20:35 +0200 Subject: Polarisation correction for extracted intensities --- src/indexamajig.c | 8 +++++++- src/pattern_sim.c | 2 +- src/peaks.c | 27 ++++++++++++++++++++------- src/peaks.h | 2 +- 4 files changed, 29 insertions(+), 10 deletions(-) (limited to 'src') 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 */ -- cgit v1.2.3