Hyperspherical coordinates

[ParadaCarleton]

Spherical coordinates can be generalized to hyperspherical coordinates, as discussed here.

[Ju-Rien]

+1

I can't find a Julia package that implements hyperspherical coordinates yet.

Knowing that these coordinates don't have a single definition for $n$ dimensions (you have to choose which hyperspherical coordinate goes from $[0, 2\pi)$ while the others go from $[0, \pi]$), it might require a more complicated approach than what is coded in CoordinateTransformations.jl.

[Ju-Rien]

So I made myself functions to move from $\mathbb{R}^n$ cartesian coordinates to $n$-spherical coordinates, if it can help out!

Cartesian to $n$-spherical:

function Ξ(x::Vector{Float64})
    n = length(x)

    if n == 1
        return x
    end
    
    only_nulls = false
    
    ϕ = zeros(n)
    ϕ[1] = norm(x) # r
    
    for i in n:-1:1
        # Cas r
        if i == 1
            continue
        # Cas n
        elseif i == n
            if (x[n] ≠ 0.0 || x[n-1] ≠ 0.0)
                ϕ[n] = 2*atan(x[n]/( x[n-1] + norm([ x[n-1], x[n] ]) ))
                only_nulls = false
            else
                ϕ[n] = 0
                only_nulls = true
            end
        else
            if only_nulls == true && x[i] == 0.0 && x[i-1] ≠ 0.0
                if x[i-1] < 0.0
                    ϕ[i] = Float64(pi)
                elseif x[i-1] > 0.0
                    ϕ[i] = 0.0
                else
                    throw("UNREACHABLE")
                end;
                only_nulls = false
            elseif only_nulls == true && x[i] == 0.0 && x[i-1] == 0.0
                ϕ[i] = 0
            else
                ϕ[i] = atan(norm([x[j] for j in i:n])/x[i-1])
            end
        end
    end
    return ϕ
end

And $n$-spherical to cartesian:

function Ξ⁻¹(rϕ)
    n = length(rϕ)
    if n == 1
        return rϕ
    elseif n == 2
        x = zeros(n)
        x[1] = rϕ[1] * cos(rϕ[2])
        x[2] = rϕ[1] * sin(rϕ[2])
        return x
    end
    r = rϕ[1]
    x = zeros(n)
    x[1] = r * cos(rϕ[2])
    for i in 2:(n-1)
        x[i] = r*prod([ sin(rϕ[j]) for j in 2:i ])*cos(rϕ[i+1])
    end
    x[n] = r*prod([ sin(rϕ[j]) for j in 2:n ])
    return x
end

[ParadaCarleton]

So I made myself functions to move from Rn cartesian coordinates to n-spherical coordinates, if it can help out!

Cartesian to n-spherical:

function Ξ(x::Vector{Float64})
    n = length(x)

    if n == 1
        return x
    end
    
    only_nulls = false
    
    ϕ = zeros(n)
    ϕ[1] = norm(x) # r
    
    for i in n:-1:1
        # Cas r
        if i == 1
            continue
        # Cas n
        elseif i == n
            if (x[n] ≠ 0.0 || x[n-1] ≠ 0.0)
                ϕ[n] = 2*atan(x[n]/( x[n-1] + norm([ x[n-1], x[n] ]) ))
                only_nulls = false
            else
                ϕ[n] = 0
                only_nulls = true
            end
        else
            if only_nulls == true && x[i] == 0.0 && x[i-1] ≠ 0.0
                if x[i-1] < 0.0
                    ϕ[i] = Float64(pi)
                elseif x[i-1] > 0.0
                    ϕ[i] = 0.0
                else
                    throw("UNREACHABLE")
                end;
                only_nulls = false
            elseif only_nulls == true && x[i] == 0.0 && x[i-1] == 0.0
                ϕ[i] = 0
            else
                ϕ[i] = atan(norm([x[j] for j in i:n])/x[i-1])
            end
        end
    end
    return ϕ
end

And n-spherical to cartesian:

function Ξ⁻¹(rϕ)
    n = length(rϕ)
    if n == 1
        return rϕ
    elseif n == 2
        x = zeros(n)
        x[1] = rϕ[1] * cos(rϕ[2])
        x[2] = rϕ[1] * sin(rϕ[2])
        return x
    end
    r = rϕ[1]
    x = zeros(n)
    x[1] = r * cos(rϕ[2])
    for i in 2:(n-1)
        x[i] = r*prod([ sin(rϕ[j]) for j in 2:i ])*cos(rϕ[i+1])
    end
    x[n] = r*prod([ sin(rϕ[j]) for j in 2:n ])
    return x
end

This looks amazing, thank you so much! Could you make a pull request adding this?

[Ju-Rien]

This looks amazing, thank you so much! Could you make a pull request adding this?

Surely! It's going to take a while (weeks?), since I'm quite busy these days. I'll have to understand the codebase and fit at the right place, while also documenting what is being done here.

[3f6a]

+1 Also need Hyperspherical coordinates, and the associated Jacobian and determinant.

[stevengj]

See this post for sketch of how to do this kind of transformation in a completely unrolled way for SVector (StaticArrays.jl) arguments.

x[i] = r*prod([ sin(rϕ[j]) for j in 2:i ])*cos(rϕ[i+1])

This is particularly inefficient: not only does it allocate a new temporary array for every loop iteration i, but it also does $O(d^2)$ multiplications and sin evaluations rather than $O(d)$.

(In my linked code, I compute sin once on each angle and then use cumprod to compute the partial products in linear time, since I don't have to worry about allocations with static-length tuples. With a dynamic length, you could just write out the cumprod loop without allocating any array other than the result.)

[Ju-Rien]

This was a good time to poke this thread, I just quit my job and have free time and energy to look at this again.

At a glance, adding hyperspherical coordinates is a "big" change for this package, since only 2D and 3D coordinate systems are implemented. I'll have to look deeper into adding an ND Hyperspherical coordinate struct/type and see if the other functionalities will adapt easily (e.g. affine maps).

Wikipedia's spherical coordinates definition seems like a good start point for the transforms.

Thank you @stevengj for pointing the inefficiency, I'll study your code snippet.

A non-exhaustive todo list:

• Add Hyperspherical struct
• Add constructors, show, isapprox functions/dispatch
• Add ::HypersphericalFromCartesian and ::CartesianFromHyperspherical struct transforms
• Add derivatives and inverse for the transforms
• Add basic compositions (compose) and conversions (convert)

Also... Hyperspherical or HyperSpherical?

[3f6a]

Also... Hyperspherical or HyperSpherical?

+1 for Hyperspherical. For example see this paper title: "Hyperspherical functions with arbitrary permutational symmetry ", https://journals.aps.org/pra/abstract/10.1103/PhysRevA.59.1135. Also I never see HyperSphere, I always see Hypersphere.

[Ju-Rien]

Small update. I'm still working on it, albeit slowly.

Current todo list:

• Add struct Hyperspherical{T, A, N}
• Add constructor Hyperspherical(r::T, φ::SVector{N_minus_1, A}) where {T, A, N_minus_1}
• Add constructor Hyperspherical(args...)
• Add constructor Hyperspherical(x::AbstractVector)
• Prepare tests for the hyperspherical coordinates transforms (TDD)
• Add functions show, isapprox functions/dispatch
• Add ::HypersphericalFromCartesian struct transforms
• Add::CartesianFromHyperspherical struct transforms
• Add derivatives and inverse for the transforms
• Add basic compositions (compose) and conversions (convert)