# 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..

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.

• Thank you, but I think there are some problems with this answer. The z-coordinate should be sin(a), and the integrand is missing a factor cos(a). More seriously I don't think that $f_1(a), f_2(a)$ are linear. Q lies on the great circle connecting P and R if OQ,OP,OR are coplanar and the determinant of these vectors is not a linear function. Maybe I did not understand your intentions? – Philip Roe Apr 1 '17 at 6:45
• No you are right, $f_1, f_2$ will not be linear functions of the polar coordinates. I am updating my answer. – Pushpendre Apr 1 '17 at 18:15
• Thanks, you took your coat off on that one! However, it is more or less what I had already done myself, and I was hoping that somebody knew of, or would discover, something a bit more slick. As I said, this is an ingredient in a numerical method for the linear wave equation. The method works nicely in two dimensions, and is very fast. This integral, and the corresponding one for quadratic P(x,y,z) would be the foundation of the method in three dimensions; they are more complicated than I had hoped. I need to take my own coat off now, but it is most useful that you have done this – Philip Roe Apr 2 '17 at 3:53
• Yeah, I too wished that I could come up with something more slick, probably using a better parameterization. From an algorithmic perspective however, the approach is pretty straightforward. It may even be possible to automate the derivation of $f_2$ and $f_1$ and the derivation of the expressions for all the integrals that need to be evaluated in Mathematica. I have never had to use complicated integrals on manifolds as part of research but I enjoyed this problem. Best of luck. – Pushpendre Apr 2 '17 at 8:45

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