import random
import math


def randomPointOnSphere(r=1., center=(0.,0.,0.)):
	"""compute random point on sphere with uniform distribution on surface.

	   Usage can be as simple as that:
	   >>> p = randomPointOnSphere()

	   We expect a point with 3 coordinates and distance 1 from the origin:
	   >>> len(p) # three coordinates
	   3
	   >>> r2 = sum(c*c for c in p) # radius squared
	   >>> round(r2, 10)
	   1.0

	   Scaling (i.e. having a radius different from 1) ...
	   >>> r = 2
	   >>> p = randomPointOnSphere(r)
	   >>> r2 = sum(c*c for c in p) # radius squared
	   >>> round(r2, 10)
	   4.0

	   ... and shifting (i.e. specifying a center different from the origin)
	   is also possible.
	   >>> center = (1,-1,2)
	   >>> p = randomPointOnSphere(center=center)
	   >>> p0 = (c-center[i] for i,c in enumerate(p)) # subtract center
	   >>> r2 = sum(c*c for c in p0) # radius squared
	   >>> round(r2, 10)
	   1.0

	   Let's test for uniform distribution by computing a bunch of random points:
	   >>> n = 2000
	   >>> pl = [randomPointOnSphere() for _ in range(n)]

	   Counting points in a cap of r/2 = 0.5 should yield roughly a quarter (0.25):
	   >>> xTopFrac = sum(p[0]>0.5 for p in pl)/n
	   >>> 0.2 < xTopFrac < 0.3
	   True
	   >>> yTopFrac = sum(p[1]>0.5 for p in pl)/n
	   >>> 0.2 <= yTopFrac <= 0.3
	   True
	   >>> zTopFrac = sum(p[2]>0.5 for p in pl)/n
	   >>> 0.2 <= zTopFrac <= 0.3
	   True
	   >>> # xTopFrac, yTopFrac, zTopFrac
	"""
	# A sphere surface has two degrees of freedom - so we need two random values.
	# The surface area of a spherical segment of a unit sphere *only* dependes on its height
	# (see http://mathworld.wolfram.com/Zone.html),
	# so when projecting points uniformly distributed on the surface to an axis,
	# the projected points are uniformly distributed in [-1,1] (our first random value).
	z0 = random.uniform(-1.,1.)
	r0_xy = math.sqrt(1. - z0*z0)
	# rotation angle around z-axis (our second random value)
	theta = 2.*math.pi*random.uniform(0., 1.)
	x0 = r0_xy*math.cos(theta)
	y0 = r0_xy*math.sin(theta)
	return (center[0] + r*x0,
	        center[1] + r*y0,
	        center[2] + r*z0)