const a1x=1/2
const a1y=-sqrt(3)/2
const a2x=1/2
const a2y=+sqrt(3)/2


function eval_expsum(kx,ky,Rmax)

  N=round(Int,ceil(2*Rmax/sqrt(3)))
 
  s=0
 
  for n1 in -N:+N
    for n2 in -N:+N
      Rx=n1*a1x+n2*a2x
      Ry=n1*a1y+n2*a2y
      if Rx^2+Ry^2 <= Rmax^2 
        s+=exp(-1im*(kx*Rx+ky*Ry))
      end
    end
  end
 
  return s
end


function eval_map(kx,ky,Rmax)

  Nx=size(kx,1)
  Ny=size(ky,1)

  s=zeros(Complex{Float64},Nx,Ny)

  for i in 1:Nx
    for j in 1:Ny
      s[i,j]=eval_expsum(kx[i],ky[j],Rmax)
    end
  end

  return s
end


kx=collect(-20:0.1:+20)
ky=collect(-20:0.1:+20)

s1=real(eval_map(kx,ky, 3.0))
s2=real(eval_map(kx,ky, 8.0))
s3=real(eval_map(kx,ky,12.0))

using Plots
pythonplot()

p1=heatmap(kx,ky,s1',aspect_ratio=:equal,colormap=:jet)
p2=heatmap(kx,ky,s2',aspect_ratio=:equal,colormap=:inferno,clims=(0,20))
p3=heatmap(kx,ky,s3',aspect_ratio=:equal,colormap=:viridis)
display(plot(p1,p2,p3,layout=(1,3)))

#plot(kx,s1[:,201],color=:green)
#plot!(kx,s2[:,201],color=:red)
#plot!(kx,s3[:,201],color=:blue)