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