The Tammes problem

I wanted to 3D print a ball with dimples like a golf ball, so I got to looking up how to evenly distribute points over the surface of a sphere. Thinking about this problem leads to a very natural optimization problem: given a natural number n, place n points on the surface of a sphere in a way that maximizes the shortest distance between any two points. This problem has a name: it is the Tammes problem. Of course, for my purposes, it really doesn't matter whether I have an exact solution to the problem--an approximate one will do.

A natural way to try to approximately solve the problem is to pretend that the points are particles that have repulsive forces between them, and then run a computer simulation of initially randomly distributed particles moving under the influence of these forces, with some frictional damping.

Inspired by this paper, I initially worked with a repulsive force inversely proportional to d^p, where d is the distance between the points and p is an exponent that is ramped up as the simulation progresses.

Experimenting with various parameters, I found it was helpful to start with p=1 and go up to p=4.5 and then stay at p=4.5 for a while before finishing the simulation. Velocity-dependent friction seems to work a little better than velocity-independent friction. The physical precision of the simulation, of course, doesn't matter at all, except as a means to getting a large minimum spacing. For 500 points, with 1000 steps of Euler-Cromer simulation and carefully tuned parameters (friction, dynamic step size, ramping schedule), I was able to get a minimum spacing of about 0.153 in one run. There is a theorem by Fejes-Tóth that implies one can't do better than 0.1702, so it's pretty close.

A hint from the above-linked paper helped along the way: one can arrange the initial positions of the particles to be symmetric around the origin (i.e., if we place one particle is at x, we place another at −x). Then we only need to simulate the motions of half of the particles, since the movements of the other half are just a reflection about the origin. Of course, this optimization only works if n is even (though if n is odd, it still may be worth arranging all but one particles symmetrically initially).

Then I had another idea. At each simulation time, we already calculate the current distance d_min between the two particles closest together, and what we want to do is to particularly strongly push apart those particles whose distance is close to that. After a fair amount of fine-tuning, I ended up modifying the repulsive force to (d−cd_min)^p, where now we ramp p from 1 to 4.5, and c from 0 to 0.9. The result was noticeably better: my best answer for n=500 went from 0.153 to about 0.162, and typical runs give me about 0.161 after only 500 steps.

In the videos, the diameters of the red ball are equal to d_min, so the point is to maximize the size of the balls without allowing them to collide. The code is written in C. It's been some time since I've programmed in C, so it was fun to go back to C. And it was also fun to go back to programming a numerical simulation, which I did a lot of back when I was a teenager. Since my teenage years, things have changed. Computers are so much faster that my ordinary laptop has no difficulty with handling n-body interaction for n around 500 or 1000: back as a teenager, the most I worked with was about n=64. Moreover, my laptop has multiple cores, and OpenMP makes it super easy to split n-body problems between cores. My Dell laptop does a 500 step calculation with 500 particles in 1.4 seconds (but for even n, we have a symmetry optimization that cuts computation time in about a half; 499 particles takes 2.6 seconds).

Enjoy the code. It's in messy but portable C, and there is a a 32-bit Windows binary that uses all the cores you have. Just give tammes one argument: the number of particles. Visualization is done by feeding an -animate option and piping to a Classic VPython script.

You can can generate an OpenSCAD golfball from the code by using a -scad option instead. Unfortunately, OpenSCAD is really slow in processing the output of tammes. The golfball on the right has n=336. I seem to have read that that's a pretty normal n for golfballs.