This is my first julia code, I am happy it did the right thing, but
compared with the Matlab code that did the same thing, it runs so slowly.
The Matlab code takes about 90s, but the julia below takes 130s.
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
function ground()
R = 320.
Ks = 2^11
x = linspace(-R, R, Ks+1)
dx = 2R/Ks
dt = 0.1
pm = 2π/2dx
px = circshift( linspace(-pm, pm, Ks+1), round(Int64, Ks/2) )
FFTW.set_num_threads(8)
V = Float64[Potential(ix, iy) for ix in x, iy in x]
ko = Float64[e^(-dt/4* (Ko(ipx,ipy))) for ipx in px, ipy in px]
ko2 = ko.^2
vo = Float64[e^(-dt* (Potential(ix,iy))) for ix in x, iy in x]
ϕ = Array(Complex128,(Ks+1,Ks+1))
Pl = plan_fft!(ϕ; flags=FFTW.MEASURE)
ϕ = V.*complex(1.)
normphi = sqrt(sumabs2(ϕ))*dx
ϕ /= normphi
ϕ = Pl*ϕ
ϕ .*= ko
ϕ = Pl\ϕ
Δϕ = 1.
nstep = 0
print("start loop")
while Δϕ > 1.e-15
ϕ′ = ϕ
ϕ .*= vo
normphi = sqrt(sumabs2(ϕ))*dx
ϕ /= normphi
ϕ = Pl*ϕ
ϕ .*= ko2
ϕ = Pl\ϕ
# if check Δϕ for every step, 35s is needed.
if nstep>500
Δϕ = maxabs(ϕ-ϕ′)
end
nstep += 1
if mod(nstep,200)==0
print(nstep," ")
print("Δϕ=",Δϕ," ")
end
end
ϕ .*= vo
ϕ = Pl*ϕ
ϕ .*= ko
ϕ = Pl\ϕ
normphi = sqrt(sumabs2(ϕ))*dx
ϕ /= normphi
end