optimise for 3D

example_figure/test_2Dtrap_vortex.jl

@@ -16,9 +16,9 @@ function new_plot(;size=(400,200))
     return p
 end
 
-## square domain
+## square domain in oscillator units
 L = 22
-N = 256
+N = 512
 μ = 30.0  
 g = 0.1
 w = 1
@@ -151,3 +151,6 @@ icd = ic_density(k, Ψ)
 
 plot(k, [ikd,ckd,nkd,qkd,cqd,iqd,icd], label=[L"n^{i}" L"n^c" L"n^k" L"n^q" L"n^{cq}" L"n^{iq}" L"n^{ic}"], xlabel=L"k", ylabel=L"n(k)", title="U=0.2", legend=:topleft)
 
+
+plot(k*ξ, [nkd,ikd+ckd+qkd+cqd+iqd+icd], label=[L"n^k" "sum"], xlabel=L"k", ylabel=L"n(k)", title="U=0.2", legend=:topleft)
+xlims!(0,1)

src/analysis.jl

@@ -250,8 +250,16 @@ function auto_correlate(ψ,X,K)
     n = length(X)
     DX,DK = fft_differentials(X,K)
     ϕ = zeropad(ψ)
-	χ = fft(ϕ)*prod(DX)
-	return ifft(abs2.(χ))*prod(DK)*(2*pi)^(n/2) |> fftshift
+	fft!(ϕ)
+    @. ϕ *= $prod(DX)
+
+    Threads.@threads for i in eachindex(ϕ)
+        ϕ[i] = abs2(ϕ[i])
+    end
+
+    ifft!(ϕ)
+    @. ϕ *= $prod(DK)*(2*π)^(n/2) 
+	return ϕ |> fftshift
 end
 
 auto_correlate(psi::Psi{D}) where D = auto_correlate(psi.ψ,psi.X,psi.K)
@@ -649,7 +657,7 @@ function iq_density(k,psi::Psi{2})
     Wi, Wc = helmholtz(ux,uy,K...)
     wix,wiy = Wi 
 
-    psia = Psi(abs.(ψ) |> complex,X,K )
+    psia = Psi(a |> complex,X,K )
     wqx,wqy = gradient(psia)
 
     U = @. exp(im*angle(ψ))