From 3c43a43685fac5a86e30a99f4dc1d4bf1baf8521 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Mon, 30 Oct 2023 21:52:28 -0600
Subject: [PATCH] Optimize phase averaging

---
 src/Intensities/Interpolation.jl          |  4 ++--
 src/SampledCorrelations/BasisReduction.jl | 11 +++++++----
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/Intensities/Interpolation.jl b/src/Intensities/Interpolation.jl
index 648d6cde8..37f9756de 100644
--- a/src/Intensities/Interpolation.jl
+++ b/src/Intensities/Interpolation.jl
@@ -148,13 +148,13 @@ function intensities_interpolated(sc::SampledCorrelations, qs, formula::Classica
     intensities = zeros(return_type, size(qs)..., nω)
     
     # Call type-stable version of the function
-    intensities_interpolated!(intensities, sc, qs, ωvals, interp, formula, stencil_info, return_type)
+    intensities_interpolated!(intensities, sc, qs, ωvals, interp, formula, stencil_info, Val(return_type))
 
     return intensities
 end
 
 
-function intensities_interpolated!(intensities, sc::SampledCorrelations, q_targets::Array, ωvals, interp::InterpolationScheme{NInterp}, formula, stencil_info, T) where {NInterp}
+function intensities_interpolated!(intensities, sc::SampledCorrelations, q_targets::Array, ωvals, interp::InterpolationScheme{NInterp}, formula, stencil_info, ::Val{T}) where {NInterp, T}
     li_intensities = LinearIndices(intensities)
     ci_targets = CartesianIndices(q_targets)
 
diff --git a/src/SampledCorrelations/BasisReduction.jl b/src/SampledCorrelations/BasisReduction.jl
index c31963259..7c57031d1 100644
--- a/src/SampledCorrelations/BasisReduction.jl
+++ b/src/SampledCorrelations/BasisReduction.jl
@@ -1,13 +1,16 @@
 function phase_averaged_elements(data, q_absolute::Vec3, crystal::Crystal, ff_atoms, ::Val{NCorr}, ::Val{NAtoms}) where {NCorr, NAtoms}
     elems = zero(MVector{NCorr,ComplexF64})
+
+    # Form factor
     ffs = ntuple(i -> compute_form_factor(ff_atoms[i], q_absolute⋅q_absolute), NAtoms)
 
-    # Real space position of each atom within the unit cell
-    rs = ntuple(i -> crystal.latvecs * crystal.positions[i], NAtoms)
+    # Overall phase factor for each site
+    q = crystal.recipvecs \ q_absolute
+    r = crystal.positions
+    prefactor = ntuple(i -> ffs[i] * exp(- 2π*im * (q ⋅ r[i])), NAtoms)
 
     for j in 1:NAtoms, i in 1:NAtoms
-        phase = exp(im*(q_absolute ⋅ (rs[j] - rs[i])))
-        elems .+= phase .* ffs[i] .* ffs[j] .* view(data, :, i, j)
+        elems .+= (prefactor[i] * conj(prefactor[j])) .* view(data, :, i, j)
     end
 
     return SVector{NCorr,ComplexF64}(elems)