From bbb2b3aca4db8227918cde9b14c0fe160b98db53 Mon Sep 17 00:00:00 2001 From: Sam Quinn Date: Fri, 18 Aug 2023 17:47:07 -0400 Subject: [PATCH 1/2] Implement the Mourigal limit --- src/Intensities/ElementContraction.jl | 29 +++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/src/Intensities/ElementContraction.jl b/src/Intensities/ElementContraction.jl index 89a51ff62..f4a8e5d72 100644 --- a/src/Intensities/ElementContraction.jl +++ b/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 ################################################################################ From 5c331b23ded284015422a765038acc22d988a9e5 Mon Sep 17 00:00:00 2001 From: Sam Quinn Date: Mon, 21 Aug 2023 13:56:05 -0400 Subject: [PATCH 2/2] Update tests to reflect unphysicality of Q=0 --- test/test_binning.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/test_binning.jl b/test/test_binning.jl index 7e69ce5aa..b6adcd4db 100644 --- a/test/test_binning.jl +++ b/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