test.jl

using StaticArrays
using PyPlot
using QuadGK
const latt = @SMatrix[1 0;
                      0 1]
const α = .1
f(r) = exp(-α*r^2)
fp(r) = -α * 2 * r * exp(-α*r^2)
Rmax = 100

# sum_{r in latt, |r| ≤ R} f(r)
function S(R, Rmax=100)
    s = 0.0
    for i = -Rmax:Rmax
        for j = -Rmax:Rmax
            r = latt * @SVector[i, j]
            nr = norm(r)
            if nr ≥ R
                s += f(nr)
            end
        end
    end
    s
end

function Sbound(R)
    diam = norm(latt * ones(2)) # /!\ assumes the lattice is positively oriented
    area = det(latt)

    # S = ∫ N'(r) V(r) ≤ N(R) V(R) + ∫ N(r) |V'(r)|
    # N(r) = number of points inside a ball of radius r
    # N(r) ≤ 1/area area of ball of radius r+diam/2
    N(r) = 1/area * π * (r+diam/2)^2
    # N(r) = 1/area * π * (r)^2

    N(R) * abs(f(R)) + quadgk(r -> N(r) * abs(fp(r)), R, Inf)[1]
end

function Sbound2(R)
    diam = norm(latt * ones(2)) # /!\ assumes the lattice is positively oriented
    area = det(latt)

    1/area * 2π * quadgk(r -> r * min(f(R), f(r)), R-diam, Inf)[1]
end

Rs = 0:20
semilogy(Rs, S.(Rs))
semilogy(Rs, Sbound.(Rs))
semilogy(Rs, Sbound2.(Rs))