function generate_G_set(b1,b2,Nmax,Gmax)

  Gset=[]

  for m1 in -Nmax:+Nmax
    for m2 in -Nmax:+Nmax
      G=m1*b1+m2*b2
      if sum(G.^2)<=Gmax^2
        push!(Gset,G)
      end
    end
  end

  return Gset
end


a=1.0

b1=4*pi/(sqrt(3)*a) * [ sqrt(3)/2, -1/2 ]
b2=4*pi/(sqrt(3)*a) * [ sqrt(3)/2, +1/2 ]

G=generate_G_set(b1,b2,50,50.0);

G2=zeros(2,size(G,1))

VG=Dict()

using QuadGK
using SpecialFunctions

V0=5 # eV
kappa=0.5 # 1/nm
r0=0.5 # nm

Vat(r)=-V0*exp(-kappa*r)*r0/(r+r0)

for i in eachindex(G)
  G2[:,i]=G[i]

  Glen=sqrt(sum(G[i].^2))
  Glen_rnd=round(Glen,digits=6)
  if !haskey(VG,Glen_rnd)
    val,err=quadgk(r -> r*Vat(r*10)*besselj0(Glen*r), 0, Inf, rtol=1e-12)
    VG[Glen_rnd]=2*pi*val
  end
end


function eval_Fourier_series(x,y,Gset,VG)

  val=0
  for G in Gset
    Glen=sqrt(sum(G.^2))
    Glen_rnd=round(Glen,digits=6)
    val+=VG[Glen_rnd]*cos(G[1]*x+G[2]*y)
  end

  return val
end

function eval_map(x,y,Gset,VG)

  Nx=size(x,1)
  Ny=size(y,1)

  V=zeros(Nx,Ny)

  for i in 1:Nx
    for j in 1:Ny
      V[i,j]=eval_Fourier_series(x[i],y[j],Gset,VG)
    end
  end

  return V
end

x=LinRange(-2,+2,101)
y=LinRange(-2,+2,101)

Vmap=eval_map(x,y,G,VG)

using Plots
pythonplot()

heatmap(x,y,Vmap')

#scatter(G2[1,:],G2[2,:])
#scatter(collect(keys(VG)),collect(values(VG)))

