sebastien · February 6, 2019 16:19 · Feb 6, 2019
diff --git a/spherical.py b/spherical.py
@@ -0,0 +1,78 @@
+# 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