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