spherical.py

# encoding: utf-8
from math import cos, sin, acos, asin, atan2, sqrt, pi

"""
A simple sketch to make sure spherical projection and unprojections work
properly. There are many ways to encode them, so this one assumes
THETA as rotating around Y (horizontal plane positionning) and PHI as the
rotation around Z (vertical positioning), so that THETA can be considered
the spherical X and PHI the spherical Y.

It's important to note that PHI is in [-90,+90], not the full range as otherwise
this will conflict with THETA.
"""

def rad(degrees):
    return degrees * pi / 180.0

def deg(radians):
    r = (radians * 180.0 / pi) % 360
    return 360 + r if r < 0 else r

def spherical( theta, phi, radius=1.0 ):
    theta, phi = rad(theta), rad(phi)
    # This assumes Theta rotating aroud Y and Phi rotating around Z
    # phi should be restriced to [-90,+90] or [-PI/2,+PI/2]
    rs = cos(phi) * radius
    x  = cos(theta) * rs
    y  = sin(phi)
    z  = sin(theta) * rs
    return (x,y,z)

def unspherical( x, y, z):
    r     = sqrt(x*x+y*y+z*z)
    theta = atan2(z,x)
    phi   = asin(y)
    return (deg(theta),deg(phi),r)

def format( p ):
    return "({0:03.2f},{1:03.2f},{2:03.2f})".format(*p)

def transform( *p ):
    sp  = spherical(*p)
    pp  = unspherical(*sp)
    spp = spherical(*pp)
    return "sph:{0} → crt:{1} → sph:{2} → crt:{3}".format(format(p), format(sp), format(pp), format(spp))


print ("phi=0")
print (transform(  0,0,1))
print (transform( 90,0,1))
print (transform(180,0,1))
print (transform(270,0,1))

print ("phi=90")
print (transform(  0, 90,1))
print (transform( 90, 90,1))
print (transform(180, 90,1))
print (transform(270, 90,1))

print ("phi=-90")
print (transform(  0, -90,1))
print (transform( 90, -90,1))
print (transform(180, -90,1))
print (transform(270, -90,1))

for theta in range(0,360):
	for phi in range(-90,90):
		phi = phi if phi > 0 else 360 + phi
		p   = (theta, phi, 1.0)
		sp  = spherical(*p)
		pp  = unspherical(*sp)
		spp = spherical(*pp)
		try:
			assert format(p) == format(pp  ), "theta={0} phi={1} {2} != {3}".format(theta, phi, format(p), format(pp))
			assert format(sp) == format(spp), "theta={0} phi={1} {2} != {3}".format(theta, phi, format(sp), format(spp))
		except Exception as e:
			print (e)
# EOF