Slow cubature

For example:

# map t ∈ (0,1)ⁿ first to hyper-spherical coordinates (ρ,θ₁,θ₂,…,θₙ₋₁)
# with a singular coordinate transform on t[1], a scaling by π of
# t[2:end-1], and a scaling by 2π of t[end], and then transforming
# to a Cartesian vector x, returning x and the Jacobian factor.
# (Throws an error if n < 2.)
function spherical_coords(t::SVector{n}) where {n}
    piT = float(eltype(t))(pi)
    ρ = t[1] / (1 - t[1]) # map [0,1) to [0,∞)
    jacobian = inv((1 - t[1])^2) * (piT^(n-1) * 2) # from coordinate scalings
    θ = (ntuple(i -> t[i+1] * piT, Val{n-2}())..., t[n] * 2piT)
    sinθ = sin.(θ)
    sinθprod = cumprod(sinθ)
    x′ = ntuple(i -> ρ * cos(θ[i+1]) * sinθprod[i], Val{n-2}())
    x = SVector(ρ * cos(θ[1]), x′..., ρ * sinθprod[n-1])
    jacobian *= ρ^(n-1) * prod(ntuple(i -> sinθ[i]^(n-1-i), Val{n-2}()))
    return x, jacobian
end

(This is just copying the formulas from the Wikipedia hyperspherical coordinates article, heavily using ntuple with Val to make sure that the compiler knows all of the lengths statically and the routine is non-allocating.) I haven’t checked it for correctness, so you should test this before using.

3 Likes