Polish + revert :full to unmodified behavior

src/Integrators.jl

@@ -92,6 +92,7 @@ end
 function step!(sys::System{0}, integrator::ImplicitMidpoint)
     s = sys.dipoles
     (; Δt, atol) = integrator
+	isnan(Δt) && return
 
     (∇E, s̄, ŝ, s̄′) = get_dipole_buffers(sys, 4)
 

src/Intensities/ElementContraction.jl

@@ -20,6 +20,7 @@ end
 
 struct FullTensor{NCorr,NSquare,NObs,NObs2} <: Contraction{SMatrix{NObs, NObs, ComplexF64, NObs2}}
     indices :: SVector{NSquare, Int64}
+    unilateral_to_bilateral :: Bool
 end
 
 struct AllAvailable{NCorr} <: Contraction{SVector{NCorr, ComplexF64}}
@@ -77,14 +78,14 @@ function Element(obs::ObservableInfo, pair::Tuple{Symbol,Symbol})
     Element(only(lookup_correlations(obs,[pair]; err_msg = pair -> "Missing correlation $(pair), which was requested.")))
 end
 
-function FullTensor(obs::ObservableInfo)
+function FullTensor(obs::ObservableInfo; unilateral_to_bilateral = true)
     n_obs = num_observables(obs)
     tensor_elements = Matrix{Tuple{Symbol,Symbol}}(undef,n_obs,n_obs)
     for (ki,i) = obs.observable_ixs, (kj,j) = obs.observable_ixs
       tensor_elements[i,j] = (ki,kj) # Required to put matrix in correct order
     end
     indices = lookup_correlations(obs, collect(tensor_elements); err_msg = αβ -> "Missing correlation $(αβ). All correlations are required to return the full tensor.")
-    FullTensor{num_correlations(obs),length(indices),n_obs,n_obs*n_obs}(indices)
+    FullTensor{num_correlations(obs),length(indices),n_obs,n_obs*n_obs}(indices, unilateral_to_bilateral)
 end
 
 ################################################################################
@@ -109,19 +110,14 @@ end
 # usually on order 1e-17 and is due to roundoff in phase_averaged_elements.
 function contract(diagonal_elements, _, traceinfo::Trace)
     if traceinfo.unilateral_to_bilateral
+        # Factor of 2 due to Laplace -> Fourier transform glueing
         2 * sum(real(diagonal_elements))
     else
         sum(real(diagonal_elements))
     end
 end
 
 function contract(dipole_elements, k::Vec3, dipoleinfo::DipoleFactor)
-    dip_factor = polarization_matrix(k)
-
-    # Note, can just take the real part since:
-    #   (1) diagonal elements are real by construction, and 
-    #   (2) pairs of off diagonal contributions have the form x*conj(y) + conj(x)*y = 2real(x*conj(y)).
-
     Sab = reshape(dipole_elements,3,3)
 
     #display(Sab)
@@ -151,6 +147,7 @@ function contract(dipole_elements, k::Vec3, dipoleinfo::DipoleFactor)
     # Since dip_factor is symmteric, and Sab is (now) guaranteed
     # to be conjugate-symmetric, dipole_intensity is real up to
     # machine precision
+    dip_factor = polarization_matrix(k)
     dipole_intensity = sum(dip_factor .* Sab)
 
     # This assertation catches the case where the user set
@@ -167,7 +164,11 @@ end
 contract(specific_element, _, ::Element) = only(specific_element)
 
 function contract(all_elems, _, full::FullTensor{NCorr,NSquare,NObs,NObs2}) where {NCorr, NSquare,NObs,NObs2}
-    reshape(all_elems,NObs,NObs)
+    Sab = reshape(all_elems,NObs,NObs)
+    if full.unilateral_to_bilateral
+        Sab = Sab + Sab'
+    end
+    Sab
 end
 
 function contract(all_elems, _, ::AllAvailable{NCorr}) where NCorr
@@ -194,13 +195,13 @@ Base.zeros(::Contraction{T}, dims...) where T = zeros(T, dims...)
 function contractor_from_mode(source, mode::Symbol)
     if mode == :trace
         contractor = Trace(source.observables; unilateral_to_bilateral = true)
-        string_formula = "Tr S′\n\n with S′ = S + S′"
+        string_formula = "Tr S′\n\n with S′ = S + S†"
     elseif mode == :perp
         contractor = DipoleFactor(source.observables; unilateral_to_bilateral = true)
         string_formula = "∑_ij (I - Q⊗Q){i,j} S′{i,j}\n\n(i,j = Sx,Sy,Sz) and with S′ = S + S†"
     elseif mode == :full
-        contractor = FullTensor(source.observables)
-        string_formula = "S{α,β}"
+        contractor = FullTensor(source.observables; unilateral_to_bilateral = true)
+        string_formula = "S′{α,β}\n\n with S′ = S + S†"
     elseif mode == :all_available
         corrs = keys(source.observables.correlations)
         contractor = AllAvailable{length(corrs)}()

src/SpinWaveTheory/SWTCalculations.jl

@@ -750,10 +750,16 @@ function intensity_formula(f::Function,swt::SpinWaveTheory,corr_ix::AbstractVect
             end
 
             # The meaning of the (A,B) correlation at this point is
-            # that the fourier transform, normalized as
+            # that the (bilateral) fourier transform, normalized as
             #   ∫±∞ exp(-iωt) <A(t) B> dt
             # of the correlation <A(t) B> has a term:
             #   2π * corrs[corr_ix_AB] * δ(disp[band] - ω)
+            #
+            # The factor of 2 (which is not included in `corrs') is
+            # incorporated later by the `unilateral_to_bilateral = true'
+            # kwarg to DipoleFactor(...), FullTensor(...), etc. The
+            # factor of π is due to the width of the Lorentzian, and
+            # Sunny currently does not incorporate this factor.
             intensity[band] = f(q_absolute, disp[band], corrs[corr_ix])
         end
 

test/test_binning.jl

@@ -109,8 +109,7 @@
 
     is2_golden = ComplexF64[0.20692326628022634 + 0.0im -0.1729875452235063 - 0.08960830762607443im -0.1321381427681472 + 0.27533711824849455im; -0.1729875452235063 + 0.08960830762607443im 0.18342229117272907 + 0.0im -0.008767695763007954 - 0.28740396674625im; -0.1321381427681472 - 0.27533711824849455im -0.008767695763007954 + 0.28740396674625im 0.4507517165000102 + 0.0im]
 
-    # is_golden is bilateral
-    @test isapprox(is[2] .* 2,is2_golden;atol = 1e-12)
+    @test isapprox(is[2],is2_golden;atol = 1e-12)
     @test all(counts .== 1.)
 
     # TODO: Test AABB

test/test_contraction.jl

@@ -1,13 +1,13 @@
 @testitem "Dipole Factor Ordering" begin
     using LinearAlgebra
     obs = Sunny.parse_observables(3; observables=nothing, correlations=nothing)
-    dipoleinfo = Sunny.DipoleFactor(obs)
-    fake_intensities = [1. 3 5; 0 7 11; 0 0 13]
+    dipoleinfo = Sunny.DipoleFactor(obs;unilateral_to_bilateral = false)
+    fake_intensities = [1. 3 5; 17 7 11; 19 23 13]
     # This is the order expected by contract(...)
-    fake_dipole_elements = fake_intensities[[1,4,5,7,8,9]]
+    fake_dipole_elements = fake_intensities[:]
     k = Sunny.Vec3(1,2,3)
     mat = Sunny.polarization_matrix(k)
     out = Sunny.contract(fake_dipole_elements,k,dipoleinfo)
-    out_true = dot(Hermitian(fake_intensities), mat)
+    out_true = dot(fake_intensities, mat)
     @test isapprox(out,out_true)
 end