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

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@#}] &

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[{[email protected], 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.

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.