Open ParadaCarleton opened 1 year ago
+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.
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
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?
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.
+1 Also need Hyperspherical coordinates, and the associated Jacobian and determinant.
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.)
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:
Hyperspherical
structshow
, isapprox
functions/dispatch::HypersphericalFromCartesian
and ::CartesianFromHyperspherical
struct transformscompose
) and conversions (convert
)Also... Hyperspherical
or HyperSpherical
?
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.
Spherical coordinates can be generalized to hyperspherical coordinates, as discussed here.