function ground()
#using PyPlot
print("start")
const R = 320.
const Ks = 2^11
const x = linspace(-R, R, Ks+1)
const dx = 2R/Ks
const y = x
const dt = 0.1
const pm = 2π/2dx
px = circshift( linspace(-pm, pm, Ks+1), round(Int64, Ks/2) )
py = px
FFTW.set_num_threads(8)
#Potential(x::Float64, y::Float64) = -2/sqrt((x-2.)^2+1)
-2/sqrt((y-2)^2+1) -1/sqrt((x+2.)^2+1) -1/sqrt((y+2)^2+1) +
1/sqrt((x-y)^2+1)
Potential(x::Float64, y::Float64) = -2./sqrt(x^2+1.) -2./sqrt(y^2+1.)
+ 1./sqrt((x-y)^2+1.)
Ko(x::Float64,y::Float64) = x^2+y^2
V = Float64[Potential(ix, iy) for ix in x, iy in y]
ϕ = V.*(1.+0im)
ϕ /= maximum(abs(ϕ))
ko = Float64[e^(-dt/4* (Ko(ipx,ipy))) for ipx in px, ipy in py]
ko2 = ko.^2
#ϕ = ground(x,V)
#function ground(x, V)
ϕ = fft(ϕ)
ϕ .*= ko
ϕ = ifft(ϕ)
Δϕ = 1.
nstep = 0
while Δϕ > 1.e-15
ϕ′ = ϕ
ϕ .*= e.^(-dt*V)
ϕ /= maximum(abs(ϕ))
ϕ = fft(ϕ)
ϕ .*= ko2
ϕ = ifft(ϕ)
Δϕ = maximum(abs(ϕ-ϕ′))
nstep += 1
if mod(nstep,200)==0
print(nstep)
print("Δϕ=",Δϕ)
end
end
ϕ .*= exp(-dt*V)
ϕ = fft(ϕ)
ϕ .*= ko
ϕ = ifft(ϕ)
ϕ /= maximum(abs(ϕ))
#imshow(abs(ϕ)^2, interpolation="none", extent=[-R, R ,-R ,R])
#show()
end
@time ground()