# how to get $n$ equidistributed points on the unit sphere

We can get $n$ equidistributed points in the unit circle using CirclePoints. But how do you get $n$ equidistributed points on the unit sphere(surface of a ball)? The preliminary idea is to suppose that the points are some electric charges with same electric quantity on a sphere. And they have the same effective working area. I think their position is what I want. This is my current solution:

point = RandomPoint[Sphere[], 10^6];
Graphics3D[{Sphere[], Blue, PointSize[.02],
Point[Union[point,SameTest -> (EuclideanDistance[#1, #2] < 0.4 &)]]},
Boxed -> False]


But this solution cannot draw the specified number of points, and the space is not very equidistributed in some places.

• The vertices of the Platonic and Archimedean solids are in a way the most optimal way to place points on a sphere. If the number of points you need is not the same as the vertex count of any of those solids, you might consider looking up spherical codes or spherical designs, or look at this article… my point is, this is a math problem and not a Mathematica problem. Jun 25, 2016 at 10:53
• What do you mean by equally spaced? (Points on a circle naturally divide the circle into segments, and we can ask that all the segments to be equal; however, points do not divide a sphere.) Jun 25, 2016 at 12:45
• Perhaps the question and your criteria need to be updated? (Something like you want a stable configuration that minimizes the potential energy of an inverse-square law, such as that for electron on a sphere. The choice of a specific objective function to minimize, as opposed to elementary geometric criteria or another function, seems an important omission. Also requiring stability rules out interesting configurations, although I don't know if that means there is always a unique solution up to rotation/reflection.) Jun 25, 2016 at 14:54
• In any event: equidistributed would be the more appropriate term than equispaced in this context. Jun 25, 2016 at 15:45
• Ifaik the problem is not solved for arbitrary number of points. Electrostatic approach gives some answers, but I don't remember a reference. Mar 18, 2018 at 9:07

Aha~ I suppose this question is created while solving this. Am I correct @yode :P

So here's an easy solution, simple, elegant, and may I say even quite fast after some optimization?

pt = With[{p =
Table[{x[i], y[i], z[i]}, {i, 80(*number of charges*)}]},
p /. Last@
NMinimize[
Total[1/Norm[Normalize[#1] - Normalize[#2]] & @@@
Subsets[p, {2}]], Flatten[p, 1]]];
Graphics3D[{Opacity@.3, Darker@Green, Sphere[], Opacity@1,
PointSize@Large, Darker@Blue, Point@*Normalize /@ pt}]


The result is quite good:

the setting of the minimization variable is crucial, or the point will not be on surface. But fortunately, our 'kindergarten physics' taught us that when charges are freely scattering in a sphere, they'll always be on surface evenly! Thus this must be some sort of 'most even' form of scattering as it follows physical laws.

• You got the idea of God. :)
– yode
Jun 25, 2016 at 14:01
• but i still cannot get the idea of upvoters~ :P
– Wjx
Jun 25, 2016 at 14:02
• This is precisely the Thomson problem mentioned by kirma. Jun 27, 2016 at 3:57
• The generated set of vectors have no unit length, so the sphere is not a unit sphere. Jul 3, 2016 at 9:36
• @plasmacel er, please Normalize them
– Wjx
Jul 3, 2016 at 9:49

If more ad hoc, inexact approaches are welcome, one way to generate relatively uniform density of points on a sphere is to use Monte Carlo Lloyd's algorithm (modified for the spherical case):

With[{points = 200, samples = 40000, iterations = 20},
Nest[With[{randoms = Join[#, RandomPoint[Sphere[], samples]]},
Table[Normalize@Mean@Extract[randoms, Position[#, {i}]], {i, points}] &@
Nearest[# -> Automatic, randoms,
DistanceFunction -> (1 - Dot[#1, #2] &)]] &,
RandomPoint[Sphere[], points], iterations]] //
Graphics3D[{Sphere[{0, 0, 0}, 0.999], Red, Point@#}] &


EDIT:

The above can be written in more concise and much more efficient form as:

With[{points = 200, samples = 40000, iterations = 20},
Nest[
With[{randoms = Join[#, RandomPoint[Sphere[], samples]]},
Normalize@Mean@randoms[[#]] & /@
Values@PositionIndex@Nearest[#, randoms]] &,
RandomPoint[Sphere[], points], iterations]] //
Graphics3D[{Sphere[{0, 0, 0}, 0.999], Red, Point@#}] &

• P.S. Nearest[# -> Range@Length@#, randoms] is more compactly done as Nearest[# -> Automatic, randoms] Jun 27, 2016 at 7:38
• @J.M. I have a strange tendency to miss these Automatic variations... I'll improve the answer. Jun 27, 2016 at 7:52
• I like the idea, but when I tried to use it in high dimensions it got quite slow. I ended up with this approach: Repeatedly find the nearest neighbour for each point; orthogonalize by subtracting the shared component; Reapeat. This converges pretty quickly even for 1000 dimensions. Dec 18, 2020 at 20:01

For an approximately even distribution of points on any surface with cylindrical symmetry, we can use the Golden Angle, the same way that the sunflower uses it on the plane.

To place N points on the surface of a sphere, define an axis. Divide the surface into N equal area strips perpendicular to the axis. For k in 0 to N-1, on the kth strip, place a point at an angle of k*ga, in the centre of the its width. ga is the golden angle, 1/(phi+1) of a circle, about 137.5 degrees / 2.34 rads.

This construction can be generalised to the surface of any volume of revolution, for instance a vase or turned table leg, by keeping the area of each strip constant.

Obviously what is being done here is that as each strip is equal area, the construction automatically makes each point 'serve' the same amount of space. Use of the 'most irrational fraction' then does a reasonable job of spreading the points round the axis without any long range structure developing.

Edit by J. M.

As I noted in a comment to this answer, the phyllotactic arrangement of points on a sphere has been previously featured on the Wolfram Blog. The code there is more general than what is needed here, so I took the liberty to simplify the code a bit for the spherical case, and also used the fact that GoldenAngle is now a built-in constant:

With[{n = Floor[4 π 100]},
Graphics3D[{Sphere[Table[{2 Sqrt[(1 - i/n) i/n] Cos[i GoldenAngle],
2 Sqrt[(1 - i/n) i/n] Sin[i GoldenAngle],
1 - 2 i/n}, {i, n}], 100/n]}, Boxed -> False]]


The 100 in the expression for n controls the point density; increase or decrease as seen fit.

• This is a Wolfram Blog entry that discusses the phyllotaxis strategy. Of note: GoldenAngle is now built-in. Jun 25, 2016 at 22:46
• Thanks for the link, it was a report on the Starshine 3 satellite that first brought this method to my attention, needless to say I couldn't google up the name when I wanted it! Jun 26, 2016 at 6:24
• If you don't mind, I'd like to edit your post later to include some of the code from that link, just to make your post even more prominent (this is a site that expects some amount of code, after all). Jun 27, 2016 at 4:14
• Go for it, I'm all for cooperation. I was surprised to see no mention of phyllotaxis so far in this thread, hence the text-only post from a non-Mathematica user (I'm a python+matplotlib+numpy man). Even given your link I was reluctant to transplant code as I'm not in a position to verify it. I was surprised to see phyllotaxis used for Starshine, as I can't imagine a hard constraint on N. I built a LED sphere recently, and triangular tiled the faces of an icoshedron, a straightforward construction, mathematically, mechanically and electrically. Jun 27, 2016 at 8:43

Correct me if I'm wrong... but I suppose, pedantically speaking, there are only five solutions defined by Platonic solids, and trivial solutions for 0-3 points (extension of CirclePoints).

Thus:

ClearAll[spherePoints];

spherePoints[r_: 1, n_ /; 0 <= n < 4] := {##, 0} & @@@ CirclePoints[r, n];

(spherePoints[r_: 1, PolyhedronData[#, "VertexCount"]] :=
r (Normalize /@ PolyhedronData[#, "VertexCoordinates"])) & /@
PolyhedronData["Platonic"];


Now:

Graphics3D[{Red, PointSize[Large], Point[spherePoints[12]],
Opacity[1/2], White, Sphere[{0, 0, 0}, 0.99]}]


For non-Platonic configurations, function stays unevaluated:

spherePoints[5]


spherePoints[5]

Tammes and Thomson problems can be considered as an extension for non-Platonic cases. For instance, see Math Overflow: Distributing points evenly on a sphere. Sadly, general solutions to these problems are not trivial.

• If I put n electric charges in the sphere,we will get the stable layout?Then we can get the layout.
– yode
Jun 25, 2016 at 8:37
• @yode I think Thompson problem is exactly that. Solving this using numerical optimization on Mma is not that hard, but I don't know how close to optimal solutions achieved this way actually are. Jun 25, 2016 at 8:45
• @yode The following numerically solves the Thomson problem for a ten-point case: With[{points = Array[c, {10, 3}]}, points /. Last@NMinimize[Total[1 / EuclideanDistance @@ # & /@ Subsets[points, {2}]], # \[Element] Sphere[] & /@ points]]. It is not really very fast, but it's reasonably simple... Jun 25, 2016 at 9:00
• And the same for Tammes problem: With[{points = Array[c, {10, 3}]}, points /. Last@NMaximize[Min[EuclideanDistance @@@ Subsets[points, {2}]], # \[Element] Sphere[] & /@ points]]. Jun 25, 2016 at 9:10
• The evaluation can be eased a bit by fixing one of the points to be a "standard" point, such as any of the unit basis vectors, and constraining a second point to lie on a coordinate plane containing the fixed point set earlier. Jun 27, 2016 at 3:47
With[{n = {3, 4, 5, 6, 7, 20}},
Partition[show[points[#]] & /@ n, 3]] // GraphicsGrid


With the code below, with piecewise spring forces (exclusively repulsive) and Cartesian vectors. A different model for force and switch to spherical vectors would probably be wise.

show[points_] :=
Graphics3D[{
Opacity[.5], Sphere[],
Opacity[1], PointSize[Medium], Point[points]},
ImageSize -> Tiny,
Boxed -> False]

points[n_] :=
With[{
steps = 50,
start = RandomPoint[Sphere[], n]},
Nest[move[#, .1] &, start, steps]]

move[points_, dt_] :=
With[{n = Length[points]},
Module[{copy = points},
Do[
copy[[i]] += dt Sum[force[copy[[j]], copy[[i]]],
{j, DeleteCases[Range[n], i]}];
copy[[i]] = copy[[i]]/Norm[copy[[i]]],
{i, n}];
copy]]

force[p1_, p2_] :=
With[{d = VectorAngle[p1, p2], range = Pi},
If[d > range, 0, Normalize[p2 - p1]*(range - d)]]


## Update: Simulated Annealing

Start with n random points on the unit sphere. Then do the main loop: pick a point randomly and make a random move. If the move lowers the energy (electro-static potential), accept the move into a new configuration. Accept the move also, this is crucial, if the energy difference is less then Exp[-difference/T], where T is a temperature-like control parameter. Do K1 * Length[points] of these inner iterations (about 10, so each point is picked about ten times) and K2 outer passes of the inner iterating. Before each outer pass T is set, i.e. gradually lowered from some starting value (few times higher than a typical temperature change) to zero.

Results agree with literature on this, so-called Thomson problem.

start = points[18];
result = Reap[simulation[start, 10, 10]]; // AbsoluteTiming

(* {2.7106, Null} *)


Energy is recorded in the main loop:

ListPlot[result[[2, 1]],
Joined -> True,
ImageSize -> Medium]


Table[Graphics3D[{
PointSize[Medium], Point[xyz /@ p],
Opacity[.75], Sphere[{0, 0, 0}, .99]},
ImageSize -> Small,
Boxed -> False,
SphericalRegion -> True],
{p, {start, result[[1]]}}] // GraphicsRow


Code:

points[n_Integer] :=
With[{x := RandomReal[]},
Table[{2 Pi \[Xi], ArcCos[2 \[Xi] - 1]}, n]]

xyz[{u_, v_}] :=
{Sin[v] Cos[u], Sin[v] Sin[u], Cos[v]}

energy[points_] :=
With[{n = Length[points]},
Sum[1/EuclideanDistance @@ (xyz /@ points[[{i, j}]]),
{i, n - 1}, {j, i + 1, n}]]

(* O(n) instead of energy O(n^2) *)
dE[points_, i_, r_] :=
With[{n = Length[points]},
Module[{indices, before, after},
indices = DeleteCases[Range[n], i];
before = Sum[
1/Norm[xyz@points[[j]] - xyz@points[[i]]], {j, indices}];
after = Sum[
1/Norm[xyz@points[[j]] - xyz@r], {j, indices}];
after - before]]

simulation[start_, K1_, K2_] :=
With[{
n = Length[start],
du = .1 Pi,
dv = .1 Pi,
x := RandomReal[]},
Module[{T, temp = start},
Do[
T = 10 (1 - (k1 - 1)/K1)^5;
temp = Nest[
Module[{i, p, change},
(* monitor energy; must Reap *)
Sow[energy[#]];
i = RandomInteger[{1, n}];
p = #[[i]] + {du (2 x - 1), dv (2 x - 1)};
change = dE[#, i, p];
If[change < 0 || (x < Exp[-change/T]),
ReplacePart[#, i -> p], #]] &, temp, K2 n], {k1, K1}];
temp]]

• Upvote.But it seem the points not very equidistributed such as this angle of view in some place.If you have fixed this,I think I should accept this answer in consideration of the code efficiency
– yode
Jun 25, 2016 at 17:01
• @yode One can tweak the steps, force etc., and see what works best. I remembered another approach which lends to this problem, see update.
– BoLe
Jul 3, 2016 at 18:19

Mathematica 11.1 has a built-in function for this task: SpherePoints. This roughly corresponds CirclePoints for three-dimensional case; thus I wouldn't rely on it to get any sort of randomness on the point placement.

EDIT: An implementation with SpherePoints, with an added random rotation:

Outer[Dot, SpherePoints[200],
RandomVariate@CircularRealMatrixDistribution[3], 1] //
Graphics3D[{Sphere[{0, 0, 0}, 0.999], Red, Point@#}] &


This works with larger number of points (10000), too:

• I think the Wolfram Research guys have seen this question,then go home to build a SpherePoints. :)
– yode
Mar 16, 2017 at 5:18
• And If you want to get a random point,you just rotate a random angle around the centre of sphere
– yode
Mar 16, 2017 at 5:20
• @yode Ah, now I understand your point. Indeed, although set of more or less equidistant points on a sphere are not absolutely random, you can improve the situation by choosing a random direction (random point on a sphere once again), and perform an uniformly distributed random rotation around it. You can alternatively pick a variate from CircularRealMatrixDistribution[3] which provides uniform samples from random 3D rotation matrices - which you can use to rotate the point set. Mar 16, 2017 at 11:17
• Yes,that is what I mean. :)
– yode
Mar 16, 2017 at 15:06

For points uniformly spaced in angular variables, you can use CirclePoints.

spherepoint[m_, n_] := Union@Flatten[Table[Join[{Cos[q]},
Sin[q] #] & /@ CirclePoints[n], {q, 0, Pi, Pi/m}], 1]

ListPointPlot3D[spherepoint[20, 30], BoxRatios -> 1]


For uniformly spaced in Cartesian coordinates, things would be complicated. The best I can suggest is to go for scaled PolyhedronData and try to map them on the sphere. You can find your polyhedron for any given number of face (or vertices, there is a relation I forgot). For example

pts = PolyhedronData["Dodecahedron", "VertexCoordinates"];
PolyhedronData["Dodecahedron", "VertexCount"]
PolyhedronData["Dodecahedron", "FaceCount"]
r = Norm[pts[[1]]];
Show[PolyhedronData["Dodecahedron"],
SphericalPlot3D[r, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},
Mesh -> Full, PlotStyle -> Opacity[0.5]]]


20

12

You have to scale the data with r to bring them on unit circle.

You might find this interesting for any arbitrary polygon

Draw an arbitrary convex polyhedron without excess diagonals drawn

• It's seem the spaced isn't equal each other?
– yode
Jun 25, 2016 at 8:06
• They are uniformly spaced in angular variable. Are you looking for uniform spacing in cartesian coordinates? Jun 25, 2016 at 8:20
• Yes,uniform spacing in cartesian coordinates.Can you help me update your solution?
– yode
Jun 25, 2016 at 8:31
• "for any given number of face (or vertices, there is a relation I forgot)." - you might want to look up dual polyhedra. ;) Jun 25, 2016 at 10:59

I posted an implementation of the spherical Lloyd algorithm in Wolfram Community. Applied to this problem, we have the following:

n = 50; (* number of points *)
BlockRandom[SeedRandom[1337, Method -> "MersenneTwister"];
sp = Normalize /@ RandomVariate[NormalDistribution[], {n, 3}]];

With[{maxit = 45,(* maximum iterations *)
tol = 0.001 (* distance tolerance *)},
lp = FixedPoint[Function[pts,
Block[{ch, polys, verts, vor},
ch = ConvexHullMesh[pts];
verts = MeshCoordinates[ch];
polys = First /@ MeshCells[ch, 2];
vor = Normalize[Cross[verts[[#2]] - verts[[#1]],
verts[[#3]] - verts[[#1]]]] &
@@@ polys;
SphericalPolygonCentroid[vor[[#]]] & /@
ch["VertexFaceConnectivity"]]], sp, maxit,
SameTest -> (Max[MapThread[cosDistance, {#1, #2}]] < tol &)]];

Graphics3D[{{Opacity[.75], Sphere[]}, {Green, Sphere[lp, 0.02]}},
Boxed -> False]


To keep this answer self contained, here are the required auxiliary routines:

(* https://mathematica.stackexchange.com/a/97854 *)
vecang[v1_?VectorQ, v2_?VectorQ] := Module[{n1 = Norm[v1], n2 = Norm[v2]},
2 ArcTan[Norm[v1 n2 + n1 v2], Norm[v1 n2 - n1 v2]]]

(* https://mathematica.stackexchange.com/a/167114 *)
vectorRotate[vv1_?VectorQ, vv2_?VectorQ] :=
Module[{v1 = Normalize[vv1], v2 = Normalize[vv2], c, d, d1, d2, t1, t2},
d = v1.v2;
If[TrueQ[Chop[1 + d] == 0],
c = UnitVector[3, First[Ordering[Abs[v1], 1]]];
t1 = c - v1; t2 = c - v2; d1 = t1.t1; d2 = t2.t2;
IdentityMatrix[3] - 2 (Outer[Times, t2, t2]/d2 -
2 t2.t1 Outer[Times, t2, t1]/(d2 d1) + Outer[Times, t1, t1]/d1),

c = Cross[v1, v2];
d IdentityMatrix[3] + Outer[Times, c, c]/(1 + d) - LeviCivitaTensor[3].c]]

(* https://mathematica.stackexchange.com/a/154112 *)
sphereExp[q_?VectorQ, p_?VectorQ] /; Length[q] == Length[p] + 1 :=
With[{n = Norm[p]}, vectorRotate[{0, 0, 1}, q].Append[p Sinc[n], Cos[n]]]

sphereLog[q_?VectorQ, p_?VectorQ] /; Length[q] == Length[p] :=
Most[vectorRotate[q, {0, 0, 1}].p]/Sinc[vecang[p, q]]

SphericalPolygonCentroid[pts_?MatrixQ] := Module[{k = 0, n = Length[pts], cp, h, pr},
cp = Normalize[Total[pts]];
pr = InternalEffectivePrecision[pts];
While[cp = sphereExp[cp, h = Sum[sphereLog[cp, p], {p, pts}]/n];
k++; Norm[h] > 10^(-pr/2) && k <= 30];
cp]

(* https://mathematica.stackexchange.com/a/159465 *)
cosDistance[v1_?VectorQ, v2_?VectorQ] :=
Module[{n1 = Normalize[v1], n2 = Normalize[v2], y},
y = Norm[n1 - n2]^2; 2 y/(Norm[n1 + n2]^2 + y)]


Of course,we have SpherePoints in 11.1,but if we are in before version,we can use IGLayoutSphere to do it,which is in a package IGraphM writed by Szabolcs.As its doucument

IGLayoutSphere[graph] lays out vertices approximately uniformly distributed on a sphere.

It's exactly what I want.

Needs["IGraphM"]
g = IGLayoutSphere[CompleteGraph[50]]


Graphics3D[{Point[GraphEmbedding[g]], Sphere[]}]


Regular equidistribution can be achieved by choosing circles of latitude at constant intervals $$dθ$$ and on these circles points with distance $$dφ$$, such that $$dθ \approx dφ$$ and that $$dθdφ$$ equals the average area per point. This then gives the following algorithm:

(In short, the basic idea is to divide the surface area into small squares (each square effectively contains one point), then find the theta, phi. Convert them into Cartesian coordinates)

lst = {};
num = 100 (*  close to required number of points*);
r = 1.;
(*Discretise theta and phi*)
ar = (4.*Pi*r^2)/num;
len = Sqrt[ar]; ndthe = Round[Pi/len]; len2 = ar/len;
For[m = 0, m <= ndthe - 1, m = m + 1,
theta = (m + 1/2)*(Pi/ndthe);
nphi = Round[2*Pi*Sin[theta]/len2];
For[n = 0, n <= nphi - 1, n = n + 1,
phi = (n)*(2*Pi/nphi);
(* spherical to cart *)
x = r*Sin[theta]*Cos[phi];
y = r*Sin[theta]*Sin[phi];
z = r*Cos[theta];
cor = {x, y, z};
corr = AppendTo[lst, cor]]]
ListPointPlot3D[corr, BoxRatios -> 1]
ListPlot3D[{corr, -corr}, BoxRatios -> 1]


Ref-Short Discription

Ref-Detailed Discription

More precisely, These points are orthogonally equidistant.

• You really don't need to use For[]: Block[{r, num, ar, len, ndthe, theta, nphi, phi}, r = 1.; num = 100; ar = (4 Pi r^2)/num; len = Sqrt[ar]; ndthe = Round[Pi/len]; Flatten[Table[theta = (m + 1/2) (Pi/ndthe); nphi = Round[2 Pi Sin[theta]/len]; Table[phi = 2 Pi n/nphi; r {Sin[theta] Cos[phi], Sin[theta] Sin[phi], Cos[theta]}, {n, 0, nphi - 1}], {m, 0, ndthe - 1}], 1]] Mar 7, 2019 at 5:28
• @J.M., Thank you! Mar 7, 2019 at 6:37

Marsaglia's method: uniform on disk (here rejection method), transform by Lambert azimuthal equal-area projection:

uniformDisk[n_] :=
Take[Select[Table[RandomReal[{-1, 1}], {3 n}, {2}], #.# <= 1 &], n];

diskToSphere[p_] := {2 p[[1]] Sqrt[1 - Dot[p, p]],
2 p[[2]] Sqrt[1 - Dot[p, p]],
1 - 2 Dot[p, p]}

Graphics3D[{
Sphere[{0, 0, 0}, .999],
Specularity[Blue, 10],
Sphere[Map[diskToSphere, uniformDisk[1256]], .015]},
Boxed -> False]

• This is more of an answer to this question, but Marsaglia's method is already there. Jun 27, 2016 at 8:49
• Didn't see it. Rereading does the OP mean "well-spaced" then and not " equidistributed"... like some low discrepancy on $\mathbb{S}^2$? Jun 27, 2016 at 9:06
• I would say that the more appropriate term is "uniformly distributed" if you are talking about (quasi)-randomly generated points. Jun 27, 2016 at 9:08

Just in case what you were really asking was how to distribute $n$ points relatively evenly on an $d$-dimensional hypersphere, here's a code snippet that will work... :)

    genSphere[n_, d_] :=
Table[
Normalize[
RandomVariate[
MultinormalDistribution[
ConstantArray[0, d],IdentityMatrix[d]
]
]
], {n}]

• How to understand this?
– yode
Apr 7, 2017 at 20:09
• Oops, had to take out the "Abs[ ]", which was putting all points in the positive quadrant. Use this to generate data to test multi-dimensional geometric algorithms. Apr 7, 2017 at 20:17
• Thanks,and do you test Graphics3D[{Point[genSphere[20, 3]], Sphere[]}]?It's not a good result still.
– yode
Apr 7, 2017 at 20:32
• The points are distributed as per a uniform distribution over the surface (each patch has same probability of getting a point as any other patch of the same size) but that inevitably results in a clumpy appearance. This could be a good seed for the simulated annealing or other iterative algorithm. It also readily generalizes to arbitrary dimension, which I find fascinating, based on its simplicity. Apr 7, 2017 at 21:22
• This is more of an answer to this question. Apr 7, 2017 at 23:33