Distributing points across the surface of a sphere

[Mark R]

I was thinking about this problem, as it comes up obliquely at work, albeit as a solved problem as part of a piece of equipment.

So, I set about thinking how to solve this. There's not a lot to find on the internet about this, so I actually had to do this from scratch.

Anyway, I had a fun couple of hours coding up a solution to this problem, and it was rather harder than I had anticipated at the start.

If anyone is interested, I'll post the code. I used c#, which is what I'm most familiar with, but python or matlab would have been a much more appropriate tool, so this little project has given me some motivation to learn these.

Spoiler

I looked to nature and how atoms are arranged in molecules and modelled this using electrostatics as a system of like-charged particles constrained to the surface of a sphere

 

[Mark R]

Anyway. I thought I'd post the code and an explanation of how it works.

By modelling the system as electrostatically repelling charges constrained to a spherical surface, it is possible to find the optimal solution. I wasn't able to find an analytical solution, but the problem is numerically solvable by steepest descent (or other algorithm).

By parameterizing the location of the points as spherical coordinates with r=1, it becomes possible to perform an unconstrained optimisation, and the resulting energy function is continuous and doubly differentiable, which greatly assists solving.

The energy function is simply the sum of 1/r where r is the straight line distance between any two points. Because the problem can be differentiated, computing the partial derivatives dE/dtheta_i and dE/dphi_i for each point (i) makes it possible to use a gradient based solver (such as steepest descent), rather than a non-gradient based technique such as simplex.

I ended up solving a slightly different problem, by requiring each point to have a mirror pair, which required minor changes to the energy function and derivative. Additionally, as I had a quasi-newton solver (limited memory Broyden–Fletcher–Goldfarb–Shanno) immediately to hand, I ended up using this which works much better than steepest descent.

 static double [] Main(int nPoints)
        {
            Random r = new Random();
            //The state vector represents concatenation of the location of the individual points in polar coordinates.
            //Odd index represent phi; even index represent theta
            double[] InitialState = new double[nPoints * 2];
            for (int i = 0; i < InitialState.Length; i++)
                InitialState[i] = r.NextDouble() * PI * 2.0 - PI;

            return new DotNumerics.Optimization.L_BFGS_B
            { AccuracyFactor = 1e1, Tolerance = 1e-12, MaxFunEvaluations = 1000 }.ComputeMin(EnergyFunction, DelEnergy, InitialState); 
        }

        //Electrostatic energy function. In this case, we assume each point has a pair opposite (i.e. a total of 2x nVectors points)
        static double EnergyFunction(double[] v)
        {
            double Energy = 0;
            for (int i = 0; i < nPoints; i++)
                for (int j = 0; j < nPoints; j++)
                {
                    if (i == j) continue;
                    //Spherical coordinates in this case are (r, theta, phi) where theta is the azimuthal coordinate
                    //and phi is the polar coordinate such that zero represents the Z axis.
                    Vector3 v1 = Vector3.FromSph(1.0, v[i * 2], v[i * 2 + 1]);
                    Vector3 v2 = Vector3.FromSph(1.0, v[j * 2], v[j * 2 + 1]);
                    Energy += 1.0 / (v1 - v2).Magnitude + 1.0 / (v1 + v2).Magnitude; // mirrored point solution
                    //Energy += 1.0 / (v1 - v2).Magnitude; // individual point solution
                }
            return Energy;
        }

        //Calculate the partial derivatives dE/dtheta_i and dE/dphi_i by summing each over j
        static double[] DelEnergy(double[] v)
        {
            const double k = 0.21022410381342863576; // 2 ^ (-9/4)
            double[] del = new double[nPoints * 2];
            for (int i = 0; i < nPoints; i++)
            {
                double theta1 = v[i * 2]; double phi1 = v[i * 2 + 1];
                double sp1 = Sin(phi1); double st1 = Sin(theta1); double ct1 = Cos(theta1);
                for (int j = 0; j < nPoints; j++)
                {
                    if(i == j) continue;
                    double theta2 = v[j * 2]; double phi2 = v[j * 2 + 1];
                    double sp2 = Sin(phi2); double cp1 = Cos(phi1); double cp2 = Cos(phi2);
                    double a = sp1 * sp2 * (st1 * Sin(theta2) + ct1 * Cos(theta2)) + cp1 * cp2;
                    double b = k * (Pow(1.0 + a, -1.5) - Pow(1.0 - a, -1.5)); // mirrored point solution
                    //double b = -k * Pow(1.0 - a, -1.5); // individual point solution
                    del[i * 2] += sp1 * sp2 * Sin(theta1 - theta2) * b;
                    del[i * 2 + 1] += (sp1 * cp2 - cp1 * sp2 * Cos(theta1 - theta2)) * b;
                }
            }
            return del;
        }

 

[Mr Evil]

Using something physical to solve a mathematical problem - very clever. It reminds me of the way you can find the ideal shape for an arch by looking at the shape of a hanging chain.

I wanted to evenly distribute some points on a sphere once. I used the vertices of a geodesic sphere, but that doesn't work for an arbitrary number of points like your solution.