# 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. – J. M. will be back soon Jun 25 '16 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.) – Michael E2 Jun 25 '16 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.) – Michael E2 Jun 25 '16 at 14:54
• In any event: equidistributed would be the more appropriate term than equispaced in this context. – J. M. will be back soon Jun 25 '16 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. – Andrew Mar 18 '18 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 '16 at 14:01
• but i still cannot get the idea of upvoters~ :P – Wjx Jun 25 '16 at 14:02
• This is precisely the Thomson problem mentioned by kirma. – J. M. will be back soon Jun 27 '16 at 3:57
• The generated set of vectors have no unit length, so the sphere is not a unit sphere. – plasmacel Jul 3 '16 at 9:36
• @plasmacel er, please Normalize them – Wjx Jul 3 '16 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] – J. M. will be back soon Jun 27 '16 at 7:38
• @J.M. I have a strange tendency to miss these Automatic variations... I'll improve the answer. – kirma Jun 27 '16 at 7:52

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],
Opacity[1/2], White, Sphere[{0, 0, 0}, 0.99]}] For non-Platonic configurations, function stays unevaluated:

spherePoints


spherePoints

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 '16 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. – kirma Jun 25 '16 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... – kirma Jun 25 '16 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]]. – kirma Jun 25 '16 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. – J. M. will be back soon Jun 27 '16 at 3:47

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. – J. M. will be back soon Jun 25 '16 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! – Neil_UK Jun 26 '16 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). – J. M. will be back soon Jun 27 '16 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. – Neil_UK Jun 27 '16 at 8:43
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, 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;
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[]}}] // 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 '16 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 '16 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,
RandomVariate@CircularRealMatrixDistribution, 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 '17 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 '17 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 which provides uniform samples from random 3D rotation matrices - which you can use to rotate the point set. – kirma Mar 16 '17 at 11:17
• Yes,that is what I mean. :) – yode Mar 16 '17 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[]];
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 '16 at 8:06
• They are uniformly spaced in angular variable. Are you looking for uniform spacing in cartesian coordinates? – Sumit Jun 25 '16 at 8:20
• Yes,uniform spacing in cartesian coordinates.Can you help me update your solution? – yode Jun 25 '16 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. ;) – J. M. will be back soon Jun 25 '16 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 - 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 + Outer[Times, c, c]/(1 + d) - LeviCivitaTensor.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] Graphics3D[{Point[GraphEmbedding[g]], Sphere[]}] 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[] Sqrt[1 - Dot[p, p]],
2 p[] Sqrt[1 - Dot[p, p]],
1 - 2 Dot[p, p]}

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

• This is more of an answer to this question, but Marsaglia's method is already there. – J. M. will be back soon Jun 27 '16 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$? – MB Reynolds Jun 27 '16 at 9:06
• I would say that the more appropriate term is "uniformly distributed" if you are talking about (quasi)-randomly generated points. – J. M. will be back soon Jun 27 '16 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 '17 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. – MikeY Apr 7 '17 at 20:17
• Thanks,and do you test Graphics3D[{Point[genSphere[20, 3]], Sphere[]}]?It's not a good result still. – yode Apr 7 '17 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. – MikeY Apr 7 '17 at 21:22
• This is more of an answer to this question. – J. M. will be back soon Apr 7 '17 at 23:33

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]] – J. M. will be back soon Mar 7 at 5:28
• @J.M.iscomputer-less, Thank you! – math Mar 7 at 6:37