In simulating interstellar dust clouds astrophysicists need high accuracy integration formulae for functions on the sphere. To construct better formulae than previously used, almost equidistantly spaced nodes on the sphere and weights belonging to these nodes are required. This problem is closely related to a facility dispersion problem on the sphere.

We propose a two-stage algorithm to compute locally optimal facility locations on the unit sphere with high accuracy and an appropriate algorithm to calculate the corresponding weights of the cubature formulae. These algorithms can be extended to other facility dispersion problems. Numerical examples show that the constructed formulae yield impressive small integration errors of approximately $10^{-12}$.

Keywords: facility dispersion, sphere, integration formulae, simulated annealing, quasi-Newton method