aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-05-15 14:45:44 +0200
committerThomas White <taw@physics.org>2013-05-15 14:45:44 +0200
commitb431c940907e5b9891b8eed7261a4a8bb922146d (patch)
tree47e224afc8d77adaade410671e59f305e246470e
parente8c8ed31440d008b34ed5cbfaa9646a6c1eac671 (diff)
process_hkl: Do polarisation correction before determining scale factors
Also, do the polarisation correction correctly..
-rw-r--r--src/process_hkl.c76
1 files changed, 48 insertions, 28 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c
index ffbc83f0..a4f23021 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -172,6 +172,48 @@ 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,
@@ -181,9 +223,12 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
Reflection *refl;
RefListIterator *iter;
double scale;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
+
+ /* First, correct for polarisation */
+ if ( !config_nopolar ) {
+ polarisation_correction(crystal_get_reflections(cr),
+ crystal_get_cell(cr), image);
+ }
if ( reference != NULL ) {
scale = scale_intensities(reference,
@@ -193,17 +238,11 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
}
if ( isnan(scale) ) return 1;
- cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
-
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
double refl_intensity, refl_sigma;
- double xl, yl, zl;
- double pol, pa, pb, phi, tt, ool;
signed int h, k, l;
int model_redundancy;
Reflection *model_version;
@@ -224,25 +263,6 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
refl_sigma = scale * get_esd_intensity(refl);
w = pow(refl_sigma, -2.0);
- /* FIXME: Should go before determining the scaling factor */
- if ( !config_nopolar ) {
-
- /* 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);
- refl_intensity /= pol;
-
- }
-
mean = get_intensity(model_version);
sumweight = get_temp1(model_version);
M2 = get_temp2(model_version);