I am trying to integrate a Lorentzian offset by a smooth function which is defined on a 3D polar domain. Here is the code:

```
using Cubature
lorentzian_re(freq, ampl, freq0, fwhm) = ampl ./ (1 + ((freq - freq0) / fwhm).^2)
singlet_0_1Hz(freq) = lorentzian_re(freq, 100.0, 0.0, 0.1)
freq_range = linspace(-2.5, +2.5, 101)
function lineshape(reltol, abstol, f=freq_range)
ρ0, z0 = 2e-1, 8e-1
field = p -> 5 * p[1] * p[3] * sin(p[2]) # hcubature chokes on this function
#field = p -> 5 * p[1] * p[3] * cos(p[2]) # hcubature integrates this without any problem
#field = p -> 5 * p[1] * sin(p[2]) # no problem with this one either
inv_vol = 1.0 / (2π * ρ0^2 * z0)
integr(polar_point, v) = begin
# integral in polar coordinates is rho * f(rho, phi, z), normalized by the integration volume
ρ_norm = polar_point[1] * inv_vol
v[:] = ρ_norm * singlet_0_1Hz(f - field(polar_point))
end
hcubature(length(f), integr, [0.0, 0.0, -z0], [ρ0, 2π, +z0], reltol=reltol, abstol=abstol)
#pcubature(length(f), integr, [0.0, 0.0, -z0], [ρ0, 2π, +z0], reltol=reltol, abstol=abstol)
end
```

The integrand is a smooth function and the convergence should not be a problem. It is not a problem when the function field is `p[1]*p[3]*cos(2p[2])`

, but hcubature runs forever on the very similar function `p[1]*p[3]*sin(2p[2])`

, with the cosine replaced by the sine. `pcubature`

has no problem with either functions.

Does anybody have any idea why?

michele