optimize intepolate_nu

src/FFCT/nugrid_interpolate.jl

@@ -1,21 +1,22 @@
 @inbounds function interpolate_nu_single!(q::T, pos::NTuple{3, T}, k_x::Vector{T}, k_y::Vector{T}, r_z::Vector{T}, us_mat::Array{T, 3}, H_r::Array{Complex{T}, 3}, uspara::USeriesPara{T}, M_mid::Int) where{T}
 
     x, y, z = pos
-    for i in 1:size(H_r, 1), j in 1:size(H_r, 2), k in 1:size(H_r, 3)
+    for i in 1:size(H_r, 1), j in 1:size(H_r, 2)
         k_xi = k_x[i]
         k_yj = k_y[j]
-        r_zk = r_z[k]
-
         phase = exp( - T(1)im * (k_xi * x + k_yj * y))
 
-        val = zero(T)
+        for k in 1:size(H_r, 3)
+            r_zk = r_z[k]
+            val = zero(T)
 
-        for l in M_mid + 1:length(uspara.sw)
-            sl, wl = uspara.sw[l]
-            val += wl * (T(2) - T(4) * (z - r_zk)^2 / sl^2 + (k_xi^2 + k_yj^2) *  sl^2) * exp(- (z - r_zk)^2 / sl^2) * us_mat[i, j, l - M_mid]
-        end
+            for l in M_mid + 1:length(uspara.sw)
+                sl, wl = uspara.sw[l]
+                val += wl * (T(2) - T(4) * (z - r_zk)^2 / sl^2 + (k_xi^2 + k_yj^2) * sl^2) * exp(- (z - r_zk)^2 / sl^2) * us_mat[i, j, l - M_mid]
+            end
 
-        H_r[i, j, k] += q * π * phase * val
+            H_r[i, j, k] += q * π * phase * val
+        end
     end 
 
     return H_r