√(x)=sqrt(x)
Δ=[1/√(2) -1/√(2) 0;1/√(6) 1/√(6) -2/√(6)]

#Representations: mat (2x2 matrix M)
#                 svd σ₁,σ₂,θ
#                 tri (3 vector of squared edge lengths)
#                 hemi λ,ϕ  (latitude, longitude)
#                 disk r,ϕ   (polar coordinates)
#                 xyz  (3 vector of cartesian coordinates on the hemisphere of radius 1/2)

# CONVERSIONS:

#mat2...

function mat2svd(M)
    (U,(σ₁,σ₂),V)=svd(M)
    # Pick the right  θ in [0,π)    
    θ=acos(V[1,1]*sign(det(M))*sign(V[1,2])) 
    (σ₁,σ₂,θ)
end

mat2tri(M) = diag( (M*Δ)'*(M*Δ)) 
mat2hemi(M) = svd2hemi( mat2svd(M)...)
mat2mat(M) = svd2mat(mat2svd(M)...) # Canonical Choice (throw away U)

function mat2disk(M)
    M/=normfro(M)
    MtM=M'*M
    rsc = [ MtM[1,2] ; ((MtM)[1,1]-(MtM)[2,2])/2]
    r=norm(rsc)
    ϕ=mod(atan2((rsc/norm(rsc))...),2π)
    (r,ϕ)
end

function mat2xyz(M)
    [M[1,1]^2+M[2,1]^2-M[1,2]^2-M[2,2]^2
        2*(M[1,1]*M[1,2]+M[2,1]*M[2,2])
        2*abs(M[1,1]*M[2,2]-M[2,1]*M[1,2])]/2
end

function mat2xyz2(M)   # Alternative in terms of MtM
    MtM=M'*M
    [MtM[1,1]-MtM[2,2]
        2*MtM[1,2] 
        2*√(det(MtM))]/2
end

#svd2...

function svd2tri(σ₁,σ₂,θ)
    (1-(σ₁^2-σ₂^2)*cos(2θ+[ 2π; -2π; 0]/3))/3
end

function svd2hemi(σ₁,σ₂,θ)
     λ=asin(2*σ₁*σ₂)
     ϕ=2θ
     λ,ϕ    
end

function svd2disk(σ₁,σ₂,θ)    
     r=(σ₁^2-σ₂^2)/2
     ϕ=2θ
     r,ϕ
end

function svd2xyz(σ₁,σ₂,θ)   
    hemi2xyz(svd2hemi(σ₁,σ₂,θ)...)
end

function svd2mat(σ₁,σ₂,θ) 
    diagm([σ₁,σ₂])*[cos(θ) sin(θ);-sin(θ) cos(θ)]
end

#hemi2...

function hemi2svd(λ,ϕ)
    σ₁=cos(λ/2)
    σ₂=sin(λ/2)
    θ=ϕ/2
    σ₁,σ₂,θ
end

function hemi2tri(λ,ϕ)
    r=cos(λ)/2
    disk2tri(r,ϕ)
end

function hemi2disk(λ,ϕ)
    r=cos(λ)/2
    r,ϕ
end

function hemi2xyz(λ,ϕ)  # Standard Spherical Coordinates
    [cos(λ)*cos(ϕ),cos(λ)*sin(ϕ),sin(λ)]/2
end

hemi2mat(λ,ϕ) = xyz2mat(hemi2xyz(λ,ϕ))
#tri2..

tri2svd(abc²) = disk2svd(tri2disk(abc²)...)

function tri2hemi(abc²)
    cs = √(6)*Δ*abc²
    λ = acos(norm(cs))
    ϕ =mod( atan2((cs/norm(cs))...), 2π)
    λ,ϕ
end

function tri2disk(abc²)
    rsc=√(3/2)*(Δ*abc²)
    r=norm(rsc)
    rsc/=r
    ϕ=mod(atan2(rsc...),2π)
    r,ϕ
end

tri2xyz(abc²) = hemi2xyz(tri2hemi(abc²)...)

tri2mat(abc²) = svd2mat(tri2svd(abc²)...)

#disk2..

function disk2svd(r,ϕ)
     σ₁= √(.5+r)
     σ₂= √(.5-r)
     θ=ϕ/2
     σ₁,σ₂,θ
end

function disk2tri(r,ϕ)
    abc² = (1-2r*cos(ϕ+[ 2π -2π 0]'/3))/3
end

function disk2hemi(r,ϕ)
    λ=acos(2r)
    λ,ϕ
end

disk2xyz(r,ϕ)=hemi2xyz(disk2hemi(r,ϕ)...)

disk2mat(r,ϕ)=hemi2mat(disk2hemi(r,ϕ)...)

# xyz2 ...

function xyz2hemi(xyz)
    λ=asin(2*xyz[3])
    ϕ=mod(atan2(xyz[2],xyz[1]),2π)
    λ,ϕ
end

xyz2svd(xyz) =  hemi2svd(xyz2hemi(xyz)...)
xyz2tri(xyz) =  hemi2tri(xyz2hemi(xyz)...)
xyz2disk(xyz) = hemi2disk(xyz2hemi(xyz)...)
xyz2mat(xyz) =  svd2mat(xyz2svd(xyz)...)

M=randn(2,2); M/=normfro(M)
println(mat2xyz(M))
println(mat2xyz2(M))

# Mat test
M=randn(2,2); M/=normfro(M)
println(mat2mat(M) )
println(svd2mat(mat2svd(M)...))
println(tri2mat(mat2tri(M)))
println(hemi2mat(mat2hemi(M)...))
println(disk2mat(mat2disk(M)...))
println(xyz2mat(mat2xyz(M)))

# SVD TEST
    sigma=sort(rand(2));sigma/=sum(sigma);
(σ₂,σ₁)=√(sigma)
θ = rand()*π
    println("-----")
    
abc   =  svd2tri(σ₁,σ₂,θ)
(λ,ϕ) = svd2hemi(σ₁,σ₂,θ)
(r,ϕ) =  svd2disk(σ₁,σ₂,θ)
xyz = svd2xyz(σ₁,σ₂,θ)
M = svd2mat(σ₁,σ₂,θ)
    println(σ₁," ",σ₂," ",θ)
    println(tri2svd(abc))
    println(hemi2svd(λ,ϕ))
    println( disk2svd(r,ϕ))
    println( xyz2svd(xyz))
     println(mat2svd(M))

# TRI TEST
  abc=rand(3); 
  abc/=sum(abc) 
    sl=sort(sqrt(abc))
    if sl[1]+sl[2]>sl[3]  # Throw away non triangle inequalities 
     (λ,ϕ) = tri2hemi(abc)   
     (σ₁,σ₂,θ)=tri2svd(abc)  
     (r,ϕ)=tri2disk(abc)
      xyz=tri2xyz(abc)
    M = tri2mat(abc)
    [abc svd2tri(σ₁,σ₂,θ) hemi2tri(λ,ϕ) disk2tri(r,ϕ) xyz2tri(xyz) mat2tri(M)]
end

# HEMI TEST 
    λ = rand()*pi/2
    ϕ = rand() *2*pi
    
   (σ₁,σ₂,θ)= hemi2svd(λ,ϕ)
    abc = hemi2tri(λ,ϕ)
    (r,ϕ) =  hemi2disk(λ,ϕ)
xyz= hemi2xyz(λ,ϕ)
M = hemi2mat(λ,ϕ)

    println(λ," ",ϕ)
    println(svd2hemi(σ₁,σ₂,θ))
    println(tri2hemi(abc))
    println( disk2hemi(r,ϕ))
    println(xyz2hemi(xyz))
println(mat2hemi(M))

# DISK TEST
    r = rand()*1/2
    ϕ = rand() *2*pi
    
    (σ₁,σ₂,θ)= disk2svd(r,ϕ)
    abc = disk2tri(r,ϕ)
    (λ,ϕ) =  disk2hemi(r,ϕ)   
xyz = disk2xyz(r,ϕ)   
M = disk2mat(r,ϕ)  
    println(r," ",ϕ)
    println(svd2disk(σ₁,σ₂,θ))
    println(tri2disk(abc))
    println(hemi2disk(λ,ϕ))
    println(xyz2disk(xyz))
    println(mat2disk(M))

#xyz test
xyz=randn(3);xyz[3]=abs(xyz[3]);xyz/=norm(xyz)*2
[xyz svd2xyz(xyz2svd(xyz)...) tri2xyz(xyz2tri(xyz)) hemi2xyz(xyz2hemi(xyz)...)  disk2xyz(xyz2disk(xyz)...)  mat2xyz(xyz2mat(xyz))]