process_hkl.jl: Add polarisation correction

1 files changed, 25 insertions, 3 deletions

diff --git a/julia/process_hkl.jl b/julia/process_hkl.jl
index 1d0bbdea..d2cebc86 100644
--- a/julia/process_hkl.jl
+++ b/julia/process_hkl.jl
@@ -1,4 +1,5 @@
 using CrystFEL
+using LinearAlgebra
 
 
 macro addmeasurement(measurement, weight,
@@ -26,7 +27,7 @@ function mergereflections(correction, crystalrefls, sym)
 
             indices = asymmetricindices(sym, refl.indices)
             model_version = get!(merged, indices)
-            @addmeasurement(correction(refl.intensity, cr), 1.0,
+            @addmeasurement(correction(refl, cr), 1.0,
                             model_version.intensity,
                             model_version.temp1,
                             model_version.temp2)
@@ -51,7 +52,28 @@ function mergereflections(correction, crystalrefls, sym)
 end
 
 
-let st = Stream("input.stream", "r")
-    merged = mergereflections((i,cr)->i, allcrystals(st), SymOpList("2/m_uab"))
+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
     savereflist!(merged, "merged.hkl")
 end