Mathematica/Uniform Spherical Distribution

In this tutorial, we will plot a set of random, uniformly distributed, points on the surface of a sphere. This seems like a trivial task, but we will see that the "obvious" solution is actually incorrect. We will start off with this incorrect method, and then improve it to be correct.

Incorrect Method

The Spherical Coordinate system we are using.

In order to plot a point on a sphere's surface, we need two numbers: the colatitude, φ and the longitude, θ. We could use the latitude, which is given by 90°−φ, but colatitude is the conventional quantity to deal with.

Now, we will make the simplistic assumption that we can use a uniform distribution for both the colatitude and longitude. We will plot 6000 points over a half sphere only, as this will allow us to see the results more clearly.

Let's define these distributions:

theta := RandomReal[ {0, Pi} ]; (* Uniform on [0,Pi], longitude, only goes halfway round*)
phi   := RandomReal[ {0, Pi} ]; (* Uniform on [0,Pi], colatitude, from pole to pole *)

You can easily alter the angle covered by changing the {0, Pi} to your desired boundaries.

Notice that we have used the "set delayed" operator (:=). We need this because when we build a table of the values in the next step, we want the Random[] function to be called for each entry. If we do not do this, the values of theta and phi will be set above and never changed, resulting in 6000 identical points. We now make the table of data:

dataPoints = Table[
   {phi, theta},
   {6000}];

Let's have a look at the distribution of the points:

rectPlot = ListPlot[
  dataPoints,
  AxesLabel -> {"\[Phi]", "\[Theta]"},
  Ticks -> {
    {0, Pi/4, Pi/2, 3 Pi/4, Pi},
    {0, Pi/4, Pi/2, 3 Pi/4, Pi}}
  ]

This produces the following output. We see that the points indeed appear to be uniformly distributed.

We will now map these spherical coordinates into the Cartesian space. This transformation is standard:

$x=r\sin \phi \cos \theta \,$

$y=r\sin \phi \sin \theta \,$

$z=r\cos \phi \,$

For our purposes, we will set r=1. We make a function called sphericalToCartesian to transform each two-element coordinate {φ,θ} in the spherical system to a three-element coordinate,{x,y,z}, in the Cartesian system.

sphericalToCartesian[{phi_, theta_}] :=
  {Sin[phi]*Cos[theta],
   Sin[phi]*Sin[theta],
   Cos[phi]};

We have a mapping, so all that remains is to apply it using the "map" operator, /@, which will perform this on each point in dataPoints:

mappedPoints = sphericalToCartesian /@ dataPoints;

Now that we have a set of points in the Cartesian system, we would like to see them. We make a set of points located at each point in mappedPoints and show it. Here, we have the default view, and orthographic views from the top and front:

 points = Graphics3D[Point /@ mappedPoints];

 Show[points]
 Show[points, ViewPoint->{0,Infinity,0}]
 Show[points, ViewPoint->{0,0,Infinity}]

We see that the points are certainly not distributed evenly - they are much more dense at the poles. This is because the mapping from spherical to Cartesian coordinates does not preserve area - the spherical space is pinched and compressed at the poles by the mapping. In order to fix this, we must alter our original distributions.

The Correct Method

In order to find the correct distribution for theta and phi, we must consider the Jacobian determinant of the mapping sphericalToCartesian. This will tell us how the spherical space is pinched by the transformation at each point.

We must find the Jacobian matrix first. We will leave r in for the moment, as Mathematica will be crunching the numbers. The evaluated Jacobian matrix is shown below, but we will let Mathematica work it out from the beginning.

$J_{F}(r,\phi ,\theta )={\begin{bmatrix}{\dfrac {\partial x}{\partial r}}&{\dfrac {\partial x}{\partial \phi }}&{\dfrac {\partial x}{\partial \theta }}\\[3pt]{\dfrac {\partial y}{\partial r}}&{\dfrac {\partial y}{\partial \phi }}&{\dfrac {\partial y}{\partial \theta }}\\[3pt]{\dfrac {\partial z}{\partial r}}&{\dfrac {\partial z}{\partial \phi }}&{\dfrac {\partial z}{\partial \theta }}\\\end{bmatrix}}={\begin{bmatrix}\sin \phi \cos \theta &r\cos \phi \cos \theta &-r\sin \phi \sin \theta \\\sin \phi \sin \theta &r\cos \phi \sin \theta &r\sin \phi \cos \theta \\\cos \phi &-r\sin \phi &0\end{bmatrix}}.$

We will construct this matrix very easily using the Outer[] function, which takes every possible combinations of each element except the first, and gives them as arguments to the function in the first element. Again, we must not evaluate this yet, as if we do, Mathematica will take the inputs f_ and x_ as "atomic," meaning that they cannot be broken down and therefore are unsuitable for inputs to Outer[].

jacMat[f_, x_] := Outer[D, f, x]

This statement will take a list of directions (i.e. r, φ, θ), and a list of the same length of functions, which are the mappings taking these to the target space, in this case, Cartesian. This simple statement doesn't check that you have given it lists as inputs, and doesn't ensure they are same length. However, we know the lengths are right, so we won't bother with this.

We now define the list of directions. First, we must remove our previous definitions of phi and theta, (and r for good measure) as at the moment, they will inject a single random number each into the above statement, rather than a symbolic quantity. We also define the list of functions, which can be obtained by giving our now-symbolic directions as arguments to sphericalToCartesian.

sphericalToCartesian[{r_, phi_, theta_}] := {r*Sin[phi]*Cos[theta],  r*Sin[phi]*Sin[theta], r*Cos[phi]};

Remove[r, phi, theta]
dirs = {r, phi, theta}
map = sphericalToCartesian[{r, phi, theta}]

In order to find the Jacobian matrix and determinant, we now put:

jacMat[map, dirs]
Simplify[Det[%]]

The output is:

{{Cos[theta] Sin[phi], -r Sin[phi] Sin[theta],  r Cos[phi] Cos[theta]}, 
 {Sin[phi] Sin[theta],  r Cos[theta] Sin[phi],  r Cos[phi] Sin[theta]}, 
 {Cos[phi],             0,                     -r Sin[phi]}}

r^2 Sin[phi]

We see that the Jacobian matrix is rearranged, but in Jacobians, it is exactly equivalent, as it just depends on the order of specification of the directions and functions. The Jacobian determinant is independent of the longitude, theta, so our uniform distribution in spherical coordinates will be uniform in Cartesian space. However, phi will not be, as as it changes, the Jacobian also changes. r also does this, but as we can just deal with a unit sphere, we can cut it out of the equation now.

There exists a method called the "Inverse CDF Method" which allows us to generate samples of a given probability density function from a set of uniformly distributed numbers between 0 and 1 (just like Mathematica's Random[]).

To use this method, we must first find the CDF of the PDF given by the Jacobian determinant. We normalise (to give a total probability of 1) by dividing by the area under the curve. <>t is a dummy variable.

area = Integrate[Sin[t], {t, 0, Pi}];(*Area under PDF*)
cdf[y_] = Integrate[Sin[t]/area, {t, 0, y}](*CDF*)

We now find the inverse of the CDF:

Solve[x == cdf[y], y](*Get distribution transform*)

{{y -> -2 ArcSin[Sqrt[x]]}, {y -> 2 ArcSin[Sqrt[x]]}}

We take the positive solution, and we can now set our theta and phi distributions:

 theta := RandomReal[ {0, Pi} ]; (* Uniform on [0,Pi], longitude, only goes halfway round*)
 phi   := 2 ArcSin[Sqrt[Random[]]]; (* Colatitude, from pole to pole *)

We now plot the points exactly as before:

dataPoints = Table[{phi, theta}, {6000}];
rectPlot = ListPlot[
  dataPoints,
  AxesLabel -> {"\[Phi]", "\[Theta]"},
  Ticks -> {
    {0, Pi/4, Pi/2, 3 Pi/4, Pi},
    {0, Pi/4, Pi/2, 3 Pi/4, Pi}}
  ]
sphericalToCartesian[{phi_, theta_}] := 
  {Sin[phi]*Cos[theta], 
   Sin[phi]*Sin[theta], 
   Cos[phi]};

mappedPoints = sphericalToCartesian /@ dataPoints;
points = Graphics3D[Point /@ mappedPoints];
Show[points]
Show[points, ViewPoint -> {0, 0, Infinity}]
Show[points, ViewPoint -> {0, Infinity, 0}]

We have reduced the distribution density around the poles.

We have now made it so that the points are evenly distributed over the spherical area.