Move polarisation_correction() to geometry.c

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,