Hcubature + spherical from CoordinateTransformations

Specific Domains
1

I am trying integration over a spherical shell by combining hcubature and spherical from CoordinateTransformations. I do get a method error. Code and error message below. Advice is much appreciated.

function foo(r)
    return CartesianFromSpherical()(r)
end 

hcubature(r -> foo(r)*r.r^2*sin(r.θ), (1,0,-pi/2), (2,2*π,π/2))

MethodError: no method matching (::CartesianFromSpherical)(::SVector{3, Float64})
Closest candidates are:
(::CartesianFromSpherical)(::Spherical) at ~/.julia/packages/CoordinateTransformations/jYKYg/src/coordinatesystems.jl:182

Stacktrace:
[1] foo(r::SVector{3, Float64})
@ Main ./In[61]:2
[2] (::var"#33#34")(r::SVector{3, Float64})
@ Main ./In[62]:1
[3] (::HCubature.GenzMalik{3, Float64})(f::var"#33#34", a::SVector{3, Float64}, b::SVector{3, Float64}, norm::typeof(norm))
@ HCubature ~/.julia/packages/HCubature/gOo1d/src/genz-malik.jl:121
[4] hcubature_(f::var"#33#34", a::SVector{3, Float64}, b::SVector{3, Float64}, norm::typeof(norm), rtol_::Int64, atol::Int64, maxevals::Int64, initdiv::Int64, buf::Nothing)
@ HCubature ~/.julia/packages/HCubature/gOo1d/src/HCubature.jl:110
[5] hcubature_(f::Function, a::Tuple{Int64, Int64, Float64}, b::Tuple{Int64, Float64, Float64}, norm::Function, rtol::Int64, atol::Int64, maxevals::Int64, initdiv::Int64, buf::Nothing)
@ HCubature ~/.julia/packages/HCubature/gOo1d/src/HCubature.jl:182
[6] hcubature#4
@ ~/.julia/packages/HCubature/gOo1d/src/HCubature.jl:235 [inlined]
[7] hcubature(f::Function, a::Tuple{Int64, Int64, Float64}, b::Tuple{Int64, Float64, Float64})
@ HCubature ~/.julia/packages/HCubature/gOo1d/src/HCubature.jl:235
[8] top-level scope
@ In[62]:1

2

Is this a good solution?

function foo(r)
    r2 = Spherical(r[1],r[2],r[3])
    return CartesianFromSpherical()(r2)
end 

hcubature(r -> foo(r)*r[1]^2*sin(r[2]), (1,0,-pi/2), (2,2*π,π/2))
1 Like
3

Seems reasonable to me, except shouldn’t your third coordinate go from [0,\pi]?

(In this particular example, I believe your integral is zero; you might want to specify an absolute tolerance atol if this is a possibility in your real problem. I would also generally specify a relative tolerance rtol, as the default tolerance 1e-8 is fairly aggressive in multidimensional problems.)

4

Many thanks again!

1/ The way I understand the definition of the spherical coordinate system in CoordinateTransform.jl is that the latitude coordinate \phi runs from -pi/2 (negative z-axis) to +pi/2 (positive z-axis) (see docstring of Spherical). But it is indeed good to cross check.

2/ Thanks for the guideline on the tolerance settings. Will do.

3/ sin(r[2]) in the Jacobian should be cos(r[3]) instead

1 Like