Implement the Mourigal limit

src/Intensities/ElementContraction.jl

@@ -83,8 +83,33 @@ end
 # Contraction helper functions
 ################################################################################
 @inline function polarization_matrix(k::Vec3)
-    k /= norm(k) + 1e-12
-    return SMatrix{3, 3, Float64, 9}(I(3) - k * k')
+    nk = norm(k)
+    if nk > 0 # Common case, k > 0
+        k /= norm(k)
+        return SMatrix{3, 3, Float64, 9}(I(3) - k * k')
+    else # Exceptional case, k = 0
+        # N.B.: Where does this 2/3 come from??
+        # =====================================
+        # When deriving the "dipole correction" (δ - q⊗q) [with q a unit vector],
+        # the following identity is used:
+        #
+        #     ∫ exp(i(Q+q)⋅R) d³R = (2π)³δ(Q+q)       [see 4.29 in Boothroyd]
+        #
+        # This integral is taken over an infinite volume (d³R over all space).
+        # In reality, the sample volume is finite, so instead of a sharp delta function
+        # on the right hand side, should really have a slightly blurred kernel whose shape
+        # depends reciprocally on the shape of the coherent volume within the sample.
+        # To zeroth order, this is just a sphere of some finite size.
+        #
+        # As a result, the (δ - q⊗q) should be averaged over nearby q in some way.
+        # For nonzero q, this changes essentially nothing. But for q comparable to
+        # 1/(sample length), the averaging is significant.
+        #
+        # Here, we assume a 3D spherical volume, so that the averaging is entirely
+        # isotropic. Averaging the 3x3 matrix diag(0,1,1) over all rotations
+        # (isotropically) gives the following matrix:
+        return SMatrix{3, 3, Float64, 9}((2/3) .* I(3)) # The Mourigal limit
+    end
 end
 
 ################################################################################

test/test_binning.jl

@@ -63,14 +63,18 @@
     is, counts = intensities_binned(sc, params, formula)
 
     is_golden = [2.452071781061995; 0.8649599530836397; 1.1585615432377976; 0.2999470844988036;;;; 0; 0; 0; 0;;;; 0; 0; 0; 0]
-    @test isapprox(is,is_golden;atol = 1e-12)
+    mask = ones(Float64,size(is_golden)) # Mask out irrelevant Q=0 intensity
+    mask[1,1,1,1] = 0
+    @test isapprox(mask .* is,mask .* is_golden;atol = 1e-12)
     @test all(counts .== 1.)
 
     is, counts = powder_average_binned(sc, (0,6π,6π/4), formula)
 
     is_golden = [4.475593277383433 0 0; 17.95271052224501 0 0; 51.13888001854976 0 0; 45.72331040682036 0 0]
+    mask = ones(Float64,size(is_golden))
+    mask[1,1] = 0
     counts_golden = [3.0 3.0 3.0; 15.0 15.0 15.0; 28.0 28.0 28.0; 39.0 39.0 39.0]
-    @test isapprox(is,is_golden;atol = 1e-12)
+    @test isapprox(mask .* is,mask .* is_golden;atol = 1e-8)
     @test isapprox(counts,counts_golden;atol = 1e-12)
 
     # TODO: Test AABB