MWE.jl

using PyPlot
using LinearAlgebra
using Cuba
using Cubature

const s = 2
const λTF = (4*(s+2)/s)^(s/(s+2))
const RTF = (4*(s+2)/s)^(1/(s+2))
const RTF_fourier = sqrt(λTF)
positive_part(x) = x > 0 ? x : zero(x)
ind(b) = b ? 1.0 : 0.0
ρ_TF(r) = positive_part(λTF - r^s) / (4π)

function ATF(x)
    r = norm(x)
    r == 0 && return @SVector[0.0, 0.0]
    @SVector[-x[2], x[1]] / (4*r^2) * (RTF^s * r^2 - 2/(s+2)*r^(s+2))
end

function m(x, p, b)
    r = norm(x)
    ind(norm(p + b*ATF(x))^2 ≤ 4π*ρ_TF(norm(x)))
end

function quad_cuba(f, RTF)
    atol = 1e-3
    res = cuhre(((x,out)->(out[1] = f(@SVector[-RTF + 2*RTF*x[1],-RTF + 2*RTF*x[2]]))),
                2,1, rtol=0.0, atol=atol,maxevals=100_000, key=7)
    res[1][1] * (2RTF)^2
end

function marg_cuba(rp, b)
    val = quad(x -> m(x, @SVector[rp, 0.0], b), RTF)
    val / (2pi)^2
end
function marg_cubature(rp, b; maxevals=500_000)
    (val, err) = hcubature(x -> m(x, @SVector[rp, 0.0], b), [-RTF, -RTF], [RTF, RTF], reltol=1e-5, abstol=0.0, maxevals=maxevals)
    val / (2pi)^2
end

figure()
ps = range(0, 2RTF_fourier, length=100)
b = 1.5
plot(ps, marg_cuba.(ps, b), label="cuba")
plot(ps, marg_cubature.(ps, b), label="cubature")