k2: wavespace over 2D domain
function jacobian(nx,ny,dx,dy,wf,k2)
eps = 1.0e-6
kx = Array{Float64}(undef,nx)
ky = Array{Float64}(undef,ny)
#wave number indexing
hx = 2.0*pi/(nx*dx)
for i = 1:Int64(nx/2)
kx[i] = hx*(i-1.0)
kx[i+Int64(nx/2)] = hx*(i-Int64(nx/2)-1)
end
kx[1] = eps
ky = transpose(kx)
j1f = zeros(ComplexF64,nx,ny)
j2f = zeros(ComplexF64,nx,ny)
j3f = zeros(ComplexF64,nx,ny)
j4f = zeros(ComplexF64,nx,ny)