diff options
| -rw-r--r-- | libcrystfel/src/geometry.c | 42 | ||||
| -rw-r--r-- | libcrystfel/src/geometry.h | 2 | ||||
| -rw-r--r-- | src/process_hkl.c | 43 |
3 files changed, 45 insertions, 42 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index cdc62ab5..efb3f2ce 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -442,3 +442,45 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel) reflist_free(predicted); } + + +void polarisation_correction(RefList *list, UnitCell *cell, struct image *image) +{ + Reflection *refl; + RefListIterator *iter; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double pol, pa, pb, phi, tt, ool; + double intensity; + double xl, yl, zl; + signed int h, k, l; + + get_indices(refl, &h, &k, &l); + + /* Polarisation correction assuming 100% polarisation + * along the x direction */ + xl = h*asx + k*bsx + l*csx; + yl = h*asy + k*bsy + l*csy; + zl = h*asz + k*bsz + l*csz; + + ool = 1.0 / image->lambda; + tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool); + phi = atan2(yl, xl); + pa = pow(sin(phi)*sin(tt), 2.0); + pb = pow(cos(tt), 2.0); + pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb); + + intensity = get_intensity(refl); + set_intensity(refl, intensity / pol); + } +} diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index 2ac5982b..a6708b2d 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -63,6 +63,8 @@ extern RefList *find_intersections(struct image *image, Crystal *cryst); extern RefList *select_intersections(struct image *image, Crystal *cryst); extern void update_partialities(Crystal *cryst, PartialityModel pmodel); +extern void polarisation_correction(RefList *list, UnitCell *cell, + struct image *image); #ifdef __cplusplus } diff --git a/src/process_hkl.c b/src/process_hkl.c index a4f23021..c2589115 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -50,6 +50,7 @@ #include "image.h" #include "crystal.h" #include "thread-pool.h" +#include "geometry.h" static void show_help(const char *s) @@ -172,48 +173,6 @@ static double scale_intensities(RefList *reference, RefList *new, } -void polarisation_correction(RefList *list, UnitCell *cell, struct image *image) -{ - Reflection *refl; - RefListIterator *iter; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double pol, pa, pb, phi, tt, ool; - double intensity; - double xl, yl, zl; - signed int h, k, l; - - get_indices(refl, &h, &k, &l); - - /* Polarisation correction assuming 100% polarisation - * along the x direction */ - xl = h*asx + k*bsx + l*csx; - yl = h*asy + k*bsy + l*csy; - zl = h*asz + k*bsz + l*csz; - - ool = 1.0 / image->lambda; - tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool); - phi = atan2(yl, xl); - pa = pow(sin(phi)*sin(tt), 2.0); - pb = pow(cos(tt), 2.0); - pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb); - - intensity = get_intensity(refl); - set_intensity(refl, intensity / pol); - } -} - - static int merge_crystal(RefList *model, struct image *image, Crystal *cr, RefList *reference, const SymOpList *sym, double *hist_vals, signed int hist_h, |