function setup_2D(ltype_val,a_val,Gmax)
  global ltype,a,dset,Gset,VG

  V0=5      # eV
  kappa=0.5 # 1/nm
  r0=0.5    # nm
  
  ltype=ltype_val
  a=a_val
  
  if ltype=="square" || ltype=="hexagonal" dset=[[0.0,0.0]] end
  
  if ltype=="honeycomb" dset=[[0.0,-a/2],[0.0,+a/2]] end
  
  Gset=generate_G_set_2D(ltype,Gmax,a)
  VG=eval_atomic_VG_2D(Gset,eval_VPC(ltype,a),V0,kappa,r0)

end


function setup_3D(ltype_val,a_val,c_val,Gmax)
  global ltype,a,c,dset,Gset,VG

  V0=5      # eV
  kappa=0.5 # 1/nm
  r0=0.5    # nm
  
  ltype=ltype_val
  a=a_val
  c=c_val

  dset=[[0.0,0.0,0.0]]

  Gset=generate_G_set_3D(ltype,Gmax,a,c)
  VG=eval_atomic_VG_3D(Gset,eval_VPC(ltype,a,c),V0,kappa,r0)

end


function plot_bands_BZpath(Nk,Nband)
  global kpath,klen,kticks,E

  if ltype!="tetragonal"
    kpath,klen,kticks=generate_BZ_path(ltype,Nk,a)
  else
    kpath,klen,kticks=generate_BZ_path(ltype,Nk,a,c)
  end
  
  @time E=eval_bands(kpath,Gset,VG,dset)
  
  plt=plot(klen,E[:,1],xticks=kticks)
  for i in 2:Nband
    plot!(plt,klen,E[:,i])
  end
  display(plt)

end


function plot_bands_surf3D(kmin,kmax,Nk,Nband)
  global kx,ky,E

  @time kx,ky,E=eval_bands_grid(kmin,kmax,Nk,Gset,VG,dset)
  
  plt=surface(kx,ky,E[:,:,1])
  for i in 2:Nband
  	surface!(plt,kx,ky,E[:,:,i])
  end
  display(plt)

end


function plot_wf_set(k,xmax,Nx,Nband)
  global eigval

  Hmtx=eval_Hmtx(k,Gset,VG,dset,1e-3)
  eigval,eigvec=eigen(Hmtx)

  x=LinRange(-xmax,+xmax,Nx)
  y=LinRange(-xmax,+xmax,Nx)

  Psi=zeros(Complex{Float64},Nx,Nx)

  for n in 1:Nband
    Psi=eval_Psi(x,y,Gset,eigvec[:,n])
    rng=maximum(abs.(Psi))
    display(plot(heatmap(x,y,real.(Psi)',clims=(-rng,+rng),color=:coolwarm,aspect_ratio=:equal), 
                 heatmap(x,y,imag.(Psi)',clims=(-rng,+rng),color=:coolwarm,aspect_ratio=:equal)))
  end

end


function plot_DOS(NBZ,NE,Nband,sigma)

  kset,wght=generate_iBZ_covering(ltype,NBZ,a)

  println("iBZ points:",size(kset,1))

  @time Eband=eval_bands(kset,Gset,VG,dset)

  Emin=minimum(Eband[:,1])-0.01
  Emax=maximum(Eband[:,Nband])+0.01

  E=LinRange(Emin,Emax,NE)

  DOS=zeros(NE)
  
  for b in 1:Nband
    for i in 1:NE
      DOS[i]+=sum(wght.*exp.(-(E[i].-Eband[:,b]).^2 ./(2*sigma^2)))/(sqrt(2*pi)*sigma)
    end
  end

  display(plot(E,DOS))

end