We develop a numerical quadrature formula for integrating a spherical homogeneous Bernstein-Bezier polynomial p over a spherical triangle \({\mathcal {T}}\) . The spherical polynomial is projected onto a planar triangle and is interpolated by a bivariate Bernstein-Bezier cubic polynomial q. Integrating q exactly allows us to express the original integral as a linear combination of the coefficients of p with the respective weights defined explicitly. The weights depend on the vertices of \({\mathcal {T}}\) and the degree of p. The approximating formula is symmetric with respect to the vertices of \({\mathcal {T}}\) . The method is easily adapted for integrating algebraically defined functions. For sufficiently smooth functions the error in the approximation is bounded by the power six of the size of the triangle. The approximation accuracy is further improved by splitting the domain triangle, effectively decreasing the error by a factor of 16 with each split. Additionally, the method conveniently couples with Richardson extrapolation, further reducing the error. Numerical experiments demonstrate the accuracy and computational effectiveness of the method.
