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 and to the theories of spherical designs and multivariate Gauss quadrature formulae. We propose a two-stage algorithm to compute 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}$.