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()

Reply via email to