Placing a satellite dish, aiming a random spotlight inside a cone-shaped beam, or scattering simulated meteorite impacts evenly across a planet's surface all reduce to the same underlying question: how do you draw a point on a sphere so that every direction is exactly as likely as every other? The obvious parametrization - latitude and longitude - turns out to be the wrong tool for the job, because sampling those two angles uniformly does not produce a uniform result. This chapter derives the correct construction from the sphere's own area element, then extends it to spherical caps aimed in any direction.
1,500
Live preview of the construction derived in this chapter: uniformly random points on a sphere. Drag to rotate, use the slider to change the sample count.
Why the Obvious Parametrization Fails
Every point on the unit sphere \(S^2\) can be written with the usual spherical coordinates, polar angle \(\theta\in[0,\pi]\) measured from the north pole and azimuth \(\varphi\in[0,2\pi)\):
\[ \mathbf{p}(\theta,\varphi)=(\sin\theta\cos\varphi,\ \sin\theta\sin\varphi,\ \cos\theta). \]Drawing \(\theta\) and \(\varphi\) independently and uniformly looks harmless, but it packs far too many points near the poles. The reason is geometric: the circles of constant \(\theta\) ("circles of latitude") shrink down to a single point as \(\theta\to0\) or \(\theta\to\pi\), so a band of fixed angular width \(d\theta\) wraps around a rapidly shrinking ring of actual surface near the poles, while the same angular width near the equator wraps around a much larger ring. Equal steps in \(\theta\) therefore do not correspond to equal patches of area on the sphere, and quantifying exactly how much they differ is the first thing to derive.
The Area Element
The area swept out by small changes \(d\theta\) and \(d\varphi\) is the magnitude of the cross product of the two tangent vectors of \(\mathbf{p}\):
\[ \frac{\partial\mathbf{p}}{\partial\theta}=(\cos\theta\cos\varphi,\ \cos\theta\sin\varphi,\ -\sin\theta), \qquad \frac{\partial\mathbf{p}}{\partial\varphi}=(-\sin\theta\sin\varphi,\ \sin\theta\cos\varphi,\ 0). \]Expanding the cross product component by component,
\[ \frac{\partial\mathbf{p}}{\partial\theta}\times\frac{\partial\mathbf{p}}{\partial\varphi} =\left(\sin^2\theta\cos\varphi,\ \sin^2\theta\sin\varphi,\ \sin\theta\cos\theta\right) =\sin\theta\cdot\mathbf{p}(\theta,\varphi). \]and since \(\mathbf{p}\) has unit length, the area element is simply
\[ \boxed{dA=\sin\theta\,d\theta\,d\varphi}. \]This is exactly the culprit: near the poles (\(\theta\) close to \(0\) or \(\pi\)), \(\sin\theta\) is small, so a fixed-size step \(d\theta\) covers very little actual area there, while the same step near the equator (\(\theta=\tfrac{\pi}{2}\)) covers much more. Sampling \(\theta\) uniformly hands out equal numbers of samples to equal-\(\theta\) bands regardless of how little area they actually contain, which is precisely what over-populates the poles.
Uniform Points on a Sphere
The fix is to stop parametrizing by \(\theta\) and switch to \(z=\cos\theta\) instead. Since \(dz=-\sin\theta\,d\theta\), and flipping the direction of integration cancels the minus sign,
\[ dA=\sin\theta\,d\theta\,d\varphi=dz\,d\varphi. \]There is no leftover factor of \(\theta\) anywhere in this expression: in the coordinates \((z,\varphi)\), area is exactly uniform. Integrating both sides over the whole sphere confirms it: \(\int_{-1}^{1}\int_0^{2\pi}dz\,d\varphi=2\cdot 2\pi=4\pi\), the surface area of the unit sphere. Geometrically, this is the classical fact that projecting a sphere radially outward onto the surface of its circumscribing cylinder preserves area exactly, known since antiquity as Archimedes' hat-box theorem.
Left: six bands at equal steps of \(30^\circ\) in \(\theta\) - thin near the poles, thick near the equator. Right: six bands at equal steps of \(z\) - identical thickness everywhere, hence identical area everywhere.
So instead of sampling angles, sample the pair \((Z,\Phi)\) uniformly on the rectangle \([-1,1]\times[0,2\pi)\) directly, with \(\Phi\) playing the role of \(\varphi\) and \(Z\) the role of \(\cos\theta\). Recovering Cartesian coordinates from this pair means inverting \(Z=\cos\theta\): since \(\theta\in[0,\pi]\), \(\sin\theta\) is never negative there, so the Pythagorean identity gives \(\sin\theta=\sqrt{1-\cos^2\theta}=\sqrt{1-Z^2}\) unambiguously, with no sign to track. Substituting \(\cos\theta=Z\), \(\sin\theta=\sqrt{1-Z^2}\) and \(\varphi=\Phi\) into the original parametrization \(\mathbf{p}(\theta,\varphi)\) from the top of this chapter recovers Cartesian coordinates without ever computing \(\theta\) itself:
\[ \boxed{ Z\sim\mathrm{unif}(-1,1),\quad \Phi\sim\mathrm{unif}(0,2\pi),\quad \mathbf{p}=\left(\sqrt{1-Z^2}\cos\Phi,\ \sqrt{1-Z^2}\sin\Phi,\ Z\right) }. \]Restricting to a Spherical Cap
A spherical cap of half-angle \(\psi\) centered on the north pole is the set of points with polar angle \(\theta\in[0,\psi]\), equivalently \(z\in[\cos\psi,\,1]\). Because \((Z,\Phi)\) are exactly uniform on the full rectangle \([-1,1]\times[0,2\pi)\), conditioning them on landing inside a sub-rectangle simply restricts each coordinate to that sub-range and renormalizes it - nothing about the shape of the density changes, since it was constant to begin with. That gives a uniform sample restricted to the cap directly, with no rejection needed:
\[ \boxed{ Z=\cos\psi+U(1-\cos\psi),\qquad \Phi=2\pi V,\qquad U,V\sim\mathrm{unif}(0,1) }. \]
The spherical cap of half-angle \(\psi\) around the north pole: exactly the band \(z\in[\cos\psi,1]\).
Aiming the Cap in an Arbitrary Direction
So far the cap sits on the north pole \(\mathbf{n}_0=(0,0,1)\). To point it along an arbitrary unit direction \(\mathbf{n}\), build a right-handed orthonormal frame \(\{\mathbf{u},\mathbf{v},\mathbf{n}\}\) and treat \(\mathbf{n}\) as the local "north pole" for that frame. Pick any reference vector \(\mathbf{e}\) not parallel to \(\mathbf{n}\) - the standard basis vector least aligned with \(\mathbf{n}\) always works and keeps the construction numerically stable - and remove its component along \(\mathbf{n}\):
\[ \mathbf{u}=\frac{\mathbf{e}-(\mathbf{e}\cdot\mathbf{n})\,\mathbf{n}}{\lVert\mathbf{e}-(\mathbf{e}\cdot\mathbf{n})\,\mathbf{n}\rVert}, \qquad \mathbf{v}=\mathbf{n}\times\mathbf{u}. \]By construction \(\mathbf{u}\perp\mathbf{n}\) and \(\lVert\mathbf{u}\rVert=1\); since \(\mathbf{v}\) is the cross product of two orthonormal vectors it is automatically a unit vector orthogonal to both. Right-handedness follows from the vector triple product identity \(\mathbf{a}\times(\mathbf{b}\times\mathbf{a})=\mathbf{b}\) for orthonormal \(\mathbf{a},\mathbf{b}\):
\[ \mathbf{u}\times\mathbf{v}=\mathbf{u}\times(\mathbf{n}\times\mathbf{u})=\mathbf{n}(\mathbf{u}\cdot\mathbf{u})-\mathbf{u}(\mathbf{u}\cdot\mathbf{n})=\mathbf{n}. \]A sample \((x,y,z)\) from the north-pole cap above then maps into world space by treating \(x\) and \(y\) as coordinates along \(\mathbf{u}\) and \(\mathbf{v}\), and \(z\) as the coordinate along \(\mathbf{n}\):
\[ \boxed{\mathbf{p}_{\mathbf{n}}=x\,\mathbf{u}+y\,\mathbf{v}+z\,\mathbf{n}}. \]Seeing the Distributions in 3D
The difference between the constructions above is easiest to judge by eye. Each sphere below is a scatter plot of 2,200 samples; drag any sphere to spin it around its own center - the three respond independently. The left sphere uses naive \((\theta,\varphi)\)-uniform sampling and visibly piles up at its poles, the middle sphere uses the \((z,\varphi)\)-uniform construction and covers the surface evenly, and the right sphere uses that same construction restricted to a cap of half-angle \(30^\circ\) aimed along an arbitrary tilted direction \(\mathbf{n}\) (marked with a black arrow), built with the orthonormal frame derived above.
naive (θ,φ)-uniform (z,φ)-uniform (correct) cap around tilted n direction n
Where This Goes Next
The same core idea - find the one substitution whose derivative reproduces the density that naive angular sampling gets wrong, so its cumulative distribution function becomes the identity - extends one dimension up to sampling uniformly random 3D orientations. See Uniform Random Rotations for that construction.
References
- Marsaglia72G. Marsaglia (1972) Choosing a Point from the Surface of a Sphere, Annals of Mathematical Statistics, 43(2), 645-646