Julia: Do polarisation correction via CrystFEL function

The Julia-native correction was 50% slower.

2 files changed, 22 insertions, 21 deletions

diff --git a/julia/CrystFEL/src/mergeutils.jl b/julia/CrystFEL/src/mergeutils.jl
index b75bd2e8..701992db 100644
--- a/julia/CrystFEL/src/mergeutils.jl
+++ b/julia/CrystFEL/src/mergeutils.jl
@@ -1,9 +1,12 @@
 module MergeUtils
 
+import ..CrystFEL: libcrystfel
 using ..CrystFEL.RefLists
 using ..CrystFEL.Symmetry
+import ..CrystFEL.UnitCells: InternalUnitCell
 export @addmeasurement, cstddev, mergereflections
 
+
 macro addmeasurement(measurement, weight,
                      mean, sumweight, wksp)
     return quote
@@ -19,12 +22,28 @@ end
 cstddev(nmeas, work1, work2) = sqrt(work2/work1)/sqrt(nmeas)
 
 
+struct Polarisation
+    fraction::Cdouble
+    angle::Cdouble
+    disable::Cint
+end
+
+function polarisation_correction!(reflist, cell, polfrac, polangle)
+    pol = Polarisation(polfrac, rad2deg(polangle), 0)
+    @ccall libcrystfel.polarisation_correction(reflist.internalptr::Ptr{InternalRefList},
+                                               cell.internalptr::Ptr{InternalUnitCell},
+                                               pol::Ref{Polarisation})::Cvoid
+end
+
+
 function mergereflections(correction, crystalrefls, sym)
 
     merged = RefList{MergedReflection}(sym)
 
     for (cr,reflections) in crystalrefls
 
+        polarisation_correction!(reflections, cr.cell, 1.0, 0.0)
+
         for refl in reflections
 
             indices = asymmetricindices(sym, refl.indices)
@@ -53,5 +72,7 @@ function mergereflections(correction, crystalrefls, sym)
 
 end
 
+mergereflections(crystalrefls, sym) = mergereflections((refl,cr)->refl.intensity, crystalrefls, sym)
+
 
 end  # of module
diff --git a/julia/process_hkl.jl b/julia/process_hkl.jl
index 80b78e52..34e6104e 100644
--- a/julia/process_hkl.jl
+++ b/julia/process_hkl.jl
@@ -1,27 +1,7 @@
 using CrystFEL
 using LinearAlgebra
 
-function anglebetween(v1, v2)
-    let v1n = norm(v1), v2n = norm(v2)
-        return 2*atan(norm(v1*v2n - v2*v1n),
-                      norm(v1*v2n + v2*v1n))
-    end
-end
-
 let st = Stream("/home/twhite/experiments/cxidb-193/short.stream", "r"),
-    merged = mergereflections(allcrystals(st), SymOpList("mmm")) do refl,crystal
-
-        polfrac = 1.0
-        polangle = 0.0
-
-        lp = transpose(crystal.cell.reciprocalcartesian) * refl.symmetricindices
-        tt = anglebetween([0,0,1], lp+[0,0,refl.kpred])
-        phi = atan(lp[2], lp[1]) - polangle
-		pol =         polfrac*(1.0 - cos(phi)*cos(phi)*sin(tt)*sin(tt)) +
-		        (1.0-polfrac)*(1.0 - sin(phi)*sin(phi)*sin(tt)*sin(tt))
-
-        return refl.intensity / pol
-
-    end
+    merged = mergereflections(allcrystals(st), SymOpList("mmm"))
     savereflist!(merged, "merged.hkl")
 end