Add generalized kernel.

src/SpinWaveTheory/Test_scripts/simple_AFM.jl

@@ -81,13 +81,35 @@ display(p)
     Plots.heatmap(qvals, energies, is'; clims=(0.0, 10), xlabel="[H,0,0]",ylabel="E (mev)")
 end
 =#
+
+#####################################
+# Old code without general Kernel
+#=
 @time begin
     qvals = 0.:0.005:1.0
     qs = [(q, 0, 0) for q in qvals]
     energies = 0.0:0.01:4.25
     is = Sunny.KPMintensities(swt_SUN, qs, energies, 200)
     p=Plots.heatmap(qvals, energies, is'; clims=(0.0, 10), xlabel="[H,0,0]",ylabel="E (mev)")
 end
+=#
+
+M = 200
+
+kT = 10 * Sunny.meV_per_K
+σ=0.1
+function lorentzian(ω, x, σ)
+    return (1/π) * (σ / ((x - ω)^2 + σ^2))
+end
+
+@time begin
+    qvals = 0.:0.005:1.0
+    qs = [(q, 0, 0) for q in qvals]
+    energies = 0.0:0.01:4.25
+    is = Sunny.KPMintensities(swt_SUN, qs, energies, M,kT,σ,lorentzian)
+    p=Plots.heatmap(qvals, energies, is'; clims=(0.0, 10), xlabel="[H,0,0]",ylabel="E (mev)")
+end
+
 
 using PlotUtils
 #savefig("20x20x1_LSWT_0p1.png")

src/SpinWaveTheory/kpm.jl

@@ -24,8 +24,20 @@ function Cheby(x,m)
     return T
 end
 
+function bose_function(kT, x) 
+    return 1 / (exp(x/kT) - 1)
+end
 
-function KPMdssf(swt::SpinWaveTheory, qs,ωlist,P::Int64)
+function get_all_coefficients(M, ωs, broadening, σ, kT,γ)
+    f(ω, x) = tanh((x/σ)^2) * broadening(ω, x*γ, σ) * (1 + bose_function(kT, x))
+    output = OffsetArray(zeros(M, length(ωs)), 0:M-1, 1:length(ωs))
+    for i in eachindex(ωs)
+        output[:, i] = cheb_coefs(M, 2M, x -> f(ωs[i], x), (-1, 1))
+    end
+    return output
+end
+
+function KPMdssf(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening)
     # This stuff all has no dependence on m:
     # P is the max Chebyshyev coefficient - relabel because 1/M expansion (Hao's code)
     (; sys, s̃_mat, T̃_mat, Q̃_mat,positions_chem) = swt
@@ -42,7 +54,7 @@ function KPMdssf(swt::SpinWaveTheory, qs,ωlist,P::Int64)
     # preallocate
     Hmat = zeros(ComplexF64, 2*nmodes, 2*nmodes) # preallocate Hmat (maybe combine with next line later)
     Avec_pref = zeros(ComplexF64, Nm) # initialize array of some prefactors   
-    chebyshev_moments = zeros(ComplexF64,3,3,length(qs),P)
+    chebyshev_moments = OffsetArray(zeros(ComplexF64,3,3,length(qs),P),1:3,1:3,1:length(qs),0:P-1)
     Sαβs = zeros(ComplexF64,3,3,length(qs),length(ωlist))
     ωdep = zeros(length(ωlist),P)
     # Calculate DSSF 
@@ -79,31 +91,39 @@ function KPMdssf(swt::SpinWaveTheory, qs,ωlist,P::Int64)
             mul!(α0,Ĩ,u[β,:]) # calculate α0
             mul!(α1,A,α0) # calculate α1   
             for α=1:3
-                chebyshev_moments[α,β,qidx,1] =  real(dot(u[α,:],α0))
-                chebyshev_moments[α,β,qidx,2] =  real(dot(u[α,:],α1)) 
+                chebyshev_moments[α,β,qidx,0] =  real(dot(u[α,:],α0))
+                chebyshev_moments[α,β,qidx,1] =  real(dot(u[α,:],α1)) 
             end
 
             for m=2:P-1
                 αnew = zeros(ComplexF64,2*nmodes) # preallocate array for αnew
                 mul!(αnew,A,α1)
                 @. αnew = 2*αnew - α0
                 for α=1:3
-                    chebyshev_moments[α,β,qidx,m+1] = real(dot(u[α,:],αnew))
+                    chebyshev_moments[α,β,qidx,m] = real(dot(u[α,:],αnew))
                 end
                 (α1, α0) = (αnew, α1)
             end
         end
+        ####################################
+        # Start of delta specific code
+        #=
         for m=0:P-1
             for w=1:length(ωlist)
                 ωdep[w,m+1]=Cheby(ωlist[w]/γ,m)*JacksonKernelDamping(m,P-1)
             end
         end
         ωdep[:,1] *=(1/π)
         ωdep[:,2:end] *=(2/π)
+        ############################
+        =#
+        ωdep =  get_all_coefficients(P,ωlist,broadening,σ,kT,γ)
+
         for w=1:length(ωlist)
             for α=1:3
                 for β=1:3
-                    Sαβs[α,β,qidx,w] = sum(chebyshev_moments[α,β,qidx,:] .*  ωdep[w,:]) * wlist[w]/γ
+                    #Sαβs[α,β,qidx,w] = sum(chebyshev_moments[α,β,qidx,:] .*  ωdep[w,:]) * wlist[w]/γ # old version for explicit δ function kernel
+                    Sαβs[α,β,qidx,w] = sum(chebyshev_moments[α,β,qidx,:] .*  ωdep[:,w])
                 end
             end
         end 
@@ -113,10 +133,10 @@ end
 
 
 
-function KPMintensities(swt::SpinWaveTheory, qs, ωvals,P::Int64)
+function KPMintensities(swt::SpinWaveTheory, qs, ωvals,P::Int64,kT,σ,broadening)
     (; sys) = swt
     qs = Sunny.Vec3.(qs)
-    Sαβs = KPMdssf(swt, qs,ωvals,P)
+    Sαβs = KPMdssf(swt, qs,ωvals,P,kT,σ,broadening)
     num_ω = length(ωvals)
     is = zeros(Float64, size(qs)..., num_ω)
     for qidx in CartesianIndices(qs)