integration on spherical triangle P(x,y,z) is a low-order (linear or quadratic) polynomial function of the Cartesian coordinates (x,y,z). T is a spherical triangle drawn on the sphere of radius r centered at the origin. T is not particularly small. Typically it covers about 5% of the sphere. I want to find the integral of P over T. An exact expression would be nice, and a numerical quadrature should be accurate to order $r^2$ or better.
The problem arises in connection with a numerical method for solving the linear wave equation..
 A: Assume that P is linear. Let $r,a,b$ be the polar coordinates corresponding to $x,y,z$, i.e. $(x,y,z) = (r\cos(a)\cos(b), r\cos(a)\sin(b), r\sin(a))$. 
Let us define, at each angle $a$, the bounds on angle $b$ for points inside a spherical triangle via functions $f_1, f_2$, i.e. the spherical triangle $\mathcal{T}$ can be defined as all angles $\mathcal{T} = \{(a,b) \mid a \in [0, A], b \in [f_1(a), f_2(a)]\}$. 
Since $P$ is linear therefore, $P(x,y,z) = P(r\cos(a)\cos(b), r\cos(a)\sin(b), r\sin(a)) = rP(\cos(a)\cos(b), \cos(a)\sin(b), \sin(a))$, and $P$ can be represented as a simple linear combination so that the integral decomposes. If $P(x,y,z) = p_1 x + p_2 y + p_3 z$ then the integral will become $r\int_0^A \int_{f_1(a)}^{f_2(a)} (p_1 \cos(a)\cos(b) + p_2 \cos(a)\sin(b) + p_3 \sin(b))\cos(a)$. So the trickiest integral will just be be $\int_0^A \cos^2(a)  (\sin(f_1(a)) - \sin(f_2(a))) da$. This integral can be numerically computed easily since it's a one dimensional integral.
Finally let us figure out what $f_1, f_2$ are:
Let $x_1, x_2$ be two vertices of $\mathcal{T}$ corresponding to polar corrdinates $(r, 0, f_1(0))$ and $(r, A, f_1(A))$. The cartesian coordinates of a point, $x$ on the arc between $x_1, x_2$ can be parameterized via an interpolation factor $\lambda \in [0, 1]$. I.e. if the polar coordinates of $x$ are $(r, \theta, f_1(\theta))$ then 
$\frac{x}{r} = \begin{bmatrix}\cos(\theta)\cos(f_1(\theta)) \\
\cos(\theta)\sin(f_1(\theta)) \\
\sin(\theta)\end{bmatrix}
=
x_\lambda / ||x_\lambda||
$
where $x_\lambda = \begin{bmatrix}
\lambda c_{11} + (1-\lambda)c_{12}\\
 \lambda c_{21} + (1-\lambda) c_{22}\\
\lambda \sin(0) + (1-\lambda) \sin(A)\end{bmatrix}$
and 
$\begin{aligned} c_{11} &= \cos(0)\cos(f_1(0)) &c_{12} &= \cos(A)\cos(f_1(A))\\ c_{21} &= \cos(0)\sin(f_1(0)) &c_{22} &= \cos(A)\sin(f_1(A))\end{aligned}$
So $f_1$ can be written as the composition of two functions $f_1^{\lambda} \circ f_1^{\theta \to \lambda}$. $f_1^{\lambda}(\lambda) = \arctan(\frac{\lambda \cos(0)\sin(f_1(0)) + (1-\lambda) \cos(A)\sin(f_1(A))}{\lambda \cos(0)\cos(f_1(0)) + (1-\lambda) \cos(A)\cos(f_1(A))})$ 
and $f_1^{\theta \to \lambda}(\theta)$ can be defined as the solution of a quadratic equation that is defined in terms of $\theta$. I.e. $\lambda$ must satisfy the following quadratic equation for any fixed $\theta$:
$$ (\frac{\sin(A)}{\sin(\theta)}(1-\lambda) \cos(\theta))^2 = (c_{12} + \lambda (c_{11} - c_{12}))^2 + (c_{22} + \lambda (c_{11} - c_{22}))^2$$ 
The appropriate root can be chosen simply by checking. So we now have a method for computing $f_1$, similarly we can write down the method for computing $f_2$ and then evaluate the whole integral.
A: This has been educational. I learn that the art of spherical geometry resides in drawing the right diagram very carefully. First try, a not very clever diagram leading to an uninspired answer. Second, an erroneous diagram leading to a beautiful but totally incorrect answer. Finally the right diagram (drawn on the outside of an orange) and a fairly pleasing answer.
The clever thing to do is the place the circumcenter of the triangle at one of the poles of the spherical coordinates. This splits the triangle into three isoceles triangles, each of which is fairly easy, although I may not yet have found the nicest way to do it.
I will not post details here unless somebody expresses interest, but if they do, I will do so with pleasure
