Calculate the area on a sphere of the intersection of two spherical caps

Given a sphere of radius $r$ with two spherical caps on it defined by the radii ($a_1$ and $a_2$) of the bases of the spherical caps, given a separation of the two spherical caps by angle $\theta$, how do you calculate the surface area of that intersection?

To clarify, the area is that of the curved surface on the sphere defined by the intersection. At the extreme where both $a_1,a_2 = r$, we would be describing a spherical lune.

Alternatively define the spherical caps by the angles $\Phi_1 = \arcsin(a_1/r)$ and $\Phi_2 = \arcsin(a_2/r)$.

You probably don't use the site anymore since the question was asked a year ago, but I'll post an answer anyway because it was a good problem and there's an ever so slight chance someone else will take a look at this.

Before I solve this, let me get some nomenclature/assumptions of my own out of the way:

1. Instead of $\theta$ being the separating angle for the two spherical caps, let it be $\alpha$.
2. When I refer to the central axis of a spherical cap (SC) I imply the axis which intersects the center of the SC's base and the center of our sphere.
3. I'll assume none of the SCs are larger than hemispheres.
4. I'll assume the SCs intersect.

Therefore, the farthest any point from the intersection of the two SCs from the center of the SC can only be $\Delta\phi = \pi/2$, using the spherical coordinate system.

The problem states that there is are two SCs with base radii $a_1$ and $a_2$ on a sphere of radius $r$ with an angle of $\alpha$ between the central axes of the two SCs.

As you mentioned, the angle between the SC's central axis and a radius from the center of the sphere to the circumference of the SC's base are $\Phi_1 = \arcsin(a_1/r)$ and $\Phi_2 = \arcsin(a_2/r)$.

Align the sphere such that SC-1's central axis and the $x$-axis are parallel (in the $\mathbb{R^3}$ coordinate system). Translate the sphere to the origin. Rotate the sphere such that the central axis of SC-2 is on the $xz$-plane. Angle $\alpha$ now rises from the $x$-axis towards the $z$-axis. Angle $\Phi_2$ is contained in $\alpha$ in the same plane, rising in the same direction (but not starting from the $x$-axis). These actions are legal since they don't distort the shape or its surface area.

Now that somewhat rigidly defined our problem we can approach it.

Let the space curve which defines the intersection of the SC-2's base with the sphere be $\vec s_2(t)$. After some (read: a lot of) vector arithmetic and trigonometry, we arrive at $\vec s_2(t)$. $$ \vec s_2(t) = \left\{\sqrt{r^2-a_2^2}\cos\alpha+r\cos(\alpha-\Phi_2)\cos t, a_2\sin t, z_2(t)\right\} $$ for $0\le t\le 2\pi$. This expression of $\vec s_2(t)$ allows us to project onto the $xy$-plane. The projection of this base-sphere intersection looks is an ellipse. The parametrization of the space curve $\vec s_1(t)$ is a bit easier because of the way we set our axes up. $$ \vec s_1(t)=\left\{\sqrt{r^2-a_1^2},a_1\sin t, a_1\cos t\right\} $$ When projected onto the $xy$-axis this looks like a line segment. Using the two projected space curves, I can find the angular limits for $\theta$ in the spherical coordinate system for 3D by finding the intersection of the ellipse and the line segment. $$ \theta_0=\arccos\frac{\sqrt{r^2-a_1^2}}{r^2-a_1^2+(a_2\sin\arccos\frac{\sqrt{r^2-a_1^2}-\sqrt{r^2-a_2^2}\cos\alpha}{r\cos(\alpha-\Phi_2)})^2} $$ It's messy, but it's a constant at least! We know that for the surface area we are interested in is for values of $\theta$ such that $-\theta_0\le \theta\le\theta_0$

Now we need to find the angular limits for $\phi$ in our spherical coordinate system. These we know to change over $\theta$ since you can imagine that as we rotate our sphere along the $z$-axis the area of interest's width changes (the angle from the $z$-axis to the position vectors of the surface area we're looking at changes).

For the top limit of $\phi$, we look at the bottom spherical cap's surface along the sphere. We use the parametrization I provided earlier. Since $z_1(t) = a_1cos(t)$ and $t$ is analogous to $\theta$ in my parametrization, I can use $\cos \phi = \frac{z}{r}$ to find that $$ \phi_1(\theta) = \arccos\frac{a_1\cos\theta}{r} $$ Similarly, using $$ z_2(t)=r\sin(\alpha-\Phi_2)+a_2(1-\cos t)\sin(\arcsin(\frac{r}{a_2}\sin\Phi_2)+\alpha) $$ which I didn't post at $\vec s_2(t)$ because of spatial concerns, I can find $$ \phi_1(\theta) = \arccos\frac{z_2(\theta)}{r} $$ Thus, for the surface area $A$ of the intersection of two spherical caps $\Sigma$, $$ A=\int\int_\Sigma dS=\int_{-\theta_0}^{\theta_0}\int_{\phi_2(\theta)}^{\phi_1(\theta)}(r^2 \sin\phi) d\phi d\theta $$ I'll be first to admit this integral is probably only numerically solvable, but I couldn't find any elegant geometric approaches to this, so I went with the "brute force" method.

EDIT

Thanks to Christian Blatter's answer to this question, I can answer in a more concise manner.

If I take our sphere situated on the origin as described before, and rotate it in the positive direction about $\hat j$ until $\alpha$ is centered on the $z$-axis, then I say:

1. Project space curves $s_1(t)$ and $s_2(t)$ onto the $xy$-plane.
2. Said curves must be ellipses or circles, since the extreme case $\alpha=\pi$ is the only one which yields straight lines upon projection and the answer to that scenario is that the surface area $A=0$.

If an ellipse is contained within another, then the surface area can be found using the formula for finding the external surface area of a spherical cap (see Christian's answer).

If the ellipses intersect at two points (they cannot intersect in this case in three or four), then using the answer to the aforementioned question, we can set up a surface integral for the calculation of $A$, given that the intersection of the ellipses on the $xy$-plane can be described by the two inequalities $a\le u\le b$ and $g(u)\le v\le f(u)$ such that $\left\{a,b:a,b\in\mathbb{R}\land a<b\right\}$, $f:\mathbb{R}\rightarrow\mathbb{R}$, $g:\mathbb{R}\rightarrow\mathbb{R}$, and $f(u)>g(u)$ over $u:[a, b]$.

Let $x=\frac{u}{\sqrt{r_x}}+c_x$ and $y=v$, for $r_x$ being the semi-diameter in the $x$ direction and $c_x$ being the horizontal shift from the origin to the center of the ellipse from the projection of $s_2(t)$ onto the $xy$-plane.

Let $\vec F=\{x, y,\sqrt{r^2-x^2-y^2}\}$, the position vector. $$ A=\int\int_\Sigma dS=\int_{u=a}^b\int_{v=f(u)}^{g(u)}\left\|\frac{\partial\vec F}{\partial v}\times\frac{\partial\vec F}{\partial u}\right\|dvdu $$

This problem can be solved by elementary geometry, resp., spherical trigonometry.

We may assume $r=1$. Let $\rho$ and $\rho'$ be the spherical radii of the two caps.

When the caps are disjoint there is nothing to compute.

When the smaller cap (of radius $\rho$) is a subset of the larger then the area $A_\cap$ of the intersection is given by $$A_\cap=4\pi\sin^2{\rho\over2}\ .$$

Now the interesting case: Assume that the two boundary circles intersect in two points $P$ and $Q$. Connect $P$ and $Q$ by an arc of a great circle; furthermore connect the centers $C$, $C'$ of the two caps with both $P$ and $Q$. The two arcs emanating from $C$ have length $\rho$ and enclose an angle $\alpha>0$, and similarly there is an angle $\alpha'$ at $C'$. The two angles $\alpha$ and $\alpha'$ showing up at $C$ and $C'$ can be computed from the given data by spherical trigonometry. I leave that to the OP.

You will then see that the intersection of the two caps (also in cases not covered by the above figure) can be written as an algebraic (i.e. $\pm$) sum of shapes of the following kind:

(a) Sectors $S$ of spherical caps of radius $\rho$, resp. $\rho'$, and central angle $\alpha>0$, resp. $\alpha'>0$. The area of such a sector is given by $$A_S=2\alpha\sin^2{\rho\over2}\ .$$ (b) Isosceles spherical triangles $\Delta$ with two sides equal to $\rho$ (or $\rho'$) and central angle $\alpha>0$. The base angle $\beta$ of such a triangle can be computed from $$\tan\beta={\cot(\alpha/2)\over\cos\rho}\ ;$$ then its area $A_\Delta$ is given by $$A_\Delta=\alpha+2\beta-\pi\ .$$

Here's a simplified formula as a function of your 3 variables, $a_1$, $a_2$, and $\theta$:

$$ 2\cos(a_2)\arccos \left ( \frac{-\cos(a_1) + \cos(\theta)\cos(a_2)}{\sin(\theta)\sin(a_2)} \right ) \\ -2\cos(a_1)\arccos \left ( \frac{\cos(a_2) - \cos(\theta)\cos(a_1)}{\sin(\theta)\sin(a_1)} \right ) \\ -2\arccos \left ( \frac{-\cos(\theta) + \cos(a_1)\cos(a_2)}{\sin(a_1)\sin(a_2)} \right ) \\ -2\pi\cos(a_2) $$

As previously stated, the caps must intersect and one cannot entirely contain the other.

This solution is copied from a graphics presentation by AMD, and is originally from a biochemistry paper. I have no proof that it works, so take it for what it's worth.

TOVCHIGRECHKO, A. AND VAKSER, I.A. 2001. How common is the funnel-like energy landscape in protein-protein interactions? Protein Sci. 10:1572-1583

Introduction

This problem was recently brought to my attention and requested to make a contribution. I will answer this question in 4 different ways. I will also numerically compare the solutions via Python code, and I will provide some 2D and 3D graphics to help illustrate the concepts.

The paper by Tovchigrechko and Vakser (TV) (mentioned and linked in another entry) describes the solution geometrically and is similar to Christian Blatter’s solution, and it is similar to the spherical trigonometry approach shown in this entry. One note about the TV paper: their cap radii are measured along the surface of the sphere unlike the OP’s request to be straight-line distances constrained to the plane of the cap’s base.

Solution

Here are the steps to solve the intersection problem (all angles are in radians throughout this entry; for more context about the givens see the OP and Figure 1 (r is the radius, a1 and a2 are straight line distances, and $\theta$ is an angular separation)). We will assume each cap has the hemisphere as its limiting case; also, as noted in other entries, we are assuming an intersection occurs where neither cap is completely contained within the other:

Figure 1: 3D graph showing the givens.

$\text{Given a1, a2, r, and }\theta:$

$\text{let a}=\sin ^{-1}\left(\text{a1}/r\right),$

$\text{b}=\sin ^{-1}\left(\text{a2}/r\right),$

$\text{c}=\theta,$

$s=\frac{1}{2} (a+b+c),$

$k=\sqrt{\frac{\sin (s-a) \sin (s-b) \sin (s-c)}{\sin (s) }},$

$A=2 \tan ^{-1}\left(\frac{k}{\sin (s-a)}\right),$

$B=2 \tan ^{-1}\left(\frac{k}{\sin (s-b)}\right),$

$C=2 \tan ^{-1}\left(\frac{k}{\sin (s-c)}\right),$

$\Delta _1=r^2 (A+B+C-\pi ),$

$\Delta _2=r^2 (A+B+C-\pi ),$

$\mathcal{A}_{\text{s1}}=2 B r^2 (1-\cos (a) ),$

$\mathcal{A}_{\text{s2}}=2 A r^2 (1-\cos (b) )$

$$ \bbox[5px,border:2px solid red] { \mathcal{A}_{\cap }=\mathcal{A}_{\text{s1}}+\mathcal{A}_{\text{s2}}-\left(\Delta _1+\Delta _2\right) } $$

That’s the answer (this, and different spherical angle formulas, can also be used to derive variants of the formula including the one found in Tovchigrechko and Vakser). The rest of this entry just provides more context and explanation.

More Info

The following is equivalent to the above steps, and is the result found in the TV paper except the $r^2$ part (a1, a2, and $\theta$ are spherical distances, i.e., distances along the curved surface on the unit sphere; to make a1 and a2 above work here, convert them like so: t=arcsin(a1/r) and u=arcsin(a2/r)):

$2 r^2 (\pi$

$-\cos ^{-1}(\csc (\text{t}) \csc (\text{u}) \cos (\theta )-\cot (\text{t}) \cot (\text{u}))$

$-\cos ^{-1}(\csc (\text{t}) \cos (\text{u}) \csc (\theta )-\cot (\text{t}) \cot (\theta ))\cos (\text{t}) $

$-\cos ^{-1}(\cos (\text{t}) \csc (\text{u}) \csc (\theta )-\cot (\text{u}) \cot (\theta ))\cos (\text{u}) )$

Figure 2: Two examples of an intersection.

Figure 2a and 2b show an intersection of 2 spherical caps. Plane I and Plane II are the planar bases of their respective caps; the planes are used to help highlight the intersection, and they'll be used later to convert the problem into 2D. We'll use Figure 2a for the following discussion.

Figure 3: Spherical sectors.

Figure 3 shows the 2 separate spherical sectors, i.e., a red area and a cyan area. The summation of the red and cyan areas contains the intersection of the spherical caps, but notice the 2 spherical triangles that get included twice when adding the red area and the cyan area. The area of these 2 spherical triangles needs to be removed from the combined area of the 2 sectors. (In this context a and b have a radius component that needs to be divided out. In the intersection formula given above the radius component is removed by the time the sector area is computed, which is why the formulas in Figure 3 appear slightly different compared to $\mathcal{A}_{\text{s1}}$ and $\mathcal{A}_{\text{s2}}$ in the intersection formula.)

Figure 4: Spherical triangle close-up.

Figure 4 shows what’s involved to calculate the area of the 2 spherical triangles whose combined area must be removed from the combined area of the 2 spherical sectors.

The angles of the spherical triangle were calculated using the half-angle formulas for tangents (equations 27 and 34-38):

/SphericalTrigonometry

And the spherical triangle area formula appears here:

/SphericalTriangle

Both topics are at Wolfram's Mathworld, but I can't make these links because my account is limited to only 8 links at the moment.

Numerical Verification (or Additional Solutions)

Cartesian Coordinates

Next, we’ll look at comparing the formula with numerical integration over a region. Imagine rotating the objects in Figure 2 such that you are looking down the line where the planes intersect. The planes would then appear like lines. Here is a 2D version of what that looks like (the red lines suggest the continuation of the planes beyond the circle):

Figure 5: 2D perspective of an intersection.

In Figure 5, a1 and a2 are radii of their respective cap shown here for context (and the radius, r, is shown to avoid confusion with $\mathcal{R}$ of the region). The blue area is the 2D projected region where the caps intersect. We will rotate (i.e., change $\theta$, a given) one cap across the span of the other cap to test the area of several different intersections as the series of images in Figure 6a show. In 6b a few perspectives of the 2 intersecting caps are shown in relation to one of the regions. One reason for doing this sweep is to consider extreme points where we might encounter numerical instability.

Figure 6: Integration in Cartesian coordinates.

Figure 6c shows the concepts involved in computing the integral of the intersecting caps in Cartesian coordinates. The double integral, with respect to the function $r / \sqrt{r^2-x^2-y^2}$, is improper and likely for that reason it presented some problems for scipy. I was able to get it to work by changing the variables in both the double integral and single integral versions, but in the code I’ve chosen the single integral version as it makes the code a little easier to read. This corresponds with the set of integrals below the dashed horizontal line in Figure 6c, which applies to the inner integral the common technique of making the substitution $y=\lvert b\rvert \sin (\theta )$. The $\sin ^{-1}\left(\frac{\lvert \sqrt{r^2-x^2} \rvert}{\lvert \sqrt{r^2-x^2} \rvert}\right)$ can be reduced to $\pi /2$.

Monte-Carlo Integration

Next, a simple approximation of the integration will be setup for the purpose of being a guard rail for catching obvious errors that might be present. In particular, a Monte-Carlo integration will be used. Here we simply remove the integrand from the double integrals in Figure 6c and integrate the scalar constant of 1 over the same regions as shown in Figure 6d. If we trust the region definitions (i.e., if they’re simple) and if the functions are well behaved, then this presents a robust approach to compute the area; for these reasons it’s the approach I used to compute the region’s area needed in the Monte-Carlo function in the Python code example. We don’t need to multiply by 2 since this is the 2D region area, not the surface area on the sphere above and below the region. The value of this region area is multiplied by the average of the $r / \sqrt{r^2-x^2-y^2}$ evaluated at each point in a random set of points within the region, which yields the Monte-Carlo approximation of the surface above the region (multiply by 2 to include the area below the region).

Integration in Spherical Coordinates

The advantage in spherical coordinates is we’ll avoid an improper integral for the most part, but a disadvantage is the area formula will appear more complex. Here’s a parametric form of a spherical cap we can use to calculate its area (it was designed for this discussion and is not intended to be flexible enough for general use (the domain of $\phi$ will depend on the cap, but for this problem $\theta$ always runs from $-\pi \text{ to } \pi$)):

$x=r \cos (\phi )$

$y=r \sin (\phi ) \cos \left(\theta f(\phi )\right)$

$z=r \sin (\phi ) \sin \left(\theta f(\phi )\right)$

$\text{where} f(\phi )=\frac{1}{\pi}\tan ^{-1}\left(\frac{\sqrt{-(b\;+\;m\;r \cos (\phi ))^2\;+\;r^2 \left(-\cos ^2(\phi )\right)\;+\;r^2}}{b\;+\;m\;r \cos (\phi )}\right)$

$\text{and}-\pi \leq \theta \leq \pi$

Figure 7: Integration in spherical coordinates.

Figure 7a is a visualization showing how the point at (1,0,0) is transformed to a/the red arc along the spherical cap in accordance with the given parametric form. Multiplication of the matrices in Figure 7a yield the above parametric form. The $f(\phi )$ part; though cluttered, is just computing the intersection at the x-value of $\phi$ along the line of a spherical cap’s base (viewing it edgewise), and then projecting that intersection upwards and downwards onto the sphere, which then yields points along the red arc in Figure 7a with respect to $\theta$ as it sweeps across from $-\pi \text{ to } \pi$ between those 2 projections, and finally sweeping it across $\phi$ yields the cap’s surface. The images in Figure 7b were computed via the parametric form above and serve as a visual verification that the parametric form is correct.

The area of a parameterized surface is equal to the magnitude of infinitesimal tangent plane normals integrated over the region, i.e.:

$\underset{\mathcal{R}}{\int }\;\lvert\lvert \mathbb{T}_u\times \mathbb{T}_v\rvert\rvert \;du\;dv$

and in our context of having x, y and z given in terms of two parameters is equivalent to (this is how it is done in the Python code in the appendix):

$\underset{\mathcal{R}}{\int }\;\sqrt{\left[\frac{\partial (x,y)}{\partial (u,v)}\right]^2+\left[\frac{\partial (y,z)}{\partial (u,v)}\right]^2+\left[\frac{\partial (x,z)}{\partial (u,v)}\right]^2}\;du\;dv=\underset{\mathcal{R}}{\int }\;\sqrt{\left\lvert \begin{array}{cc} \frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} \\ \frac{\partial y}{\partial u} & \frac{\partial y}{\partial v} \\ \end{array} \right\rvert ^2+\left\lvert \begin{array}{cc} \frac{\partial y}{\partial u} & \frac{\partial y}{\partial v} \\ \frac{\partial z}{\partial u} & \frac{\partial z}{\partial v} \\ \end{array} \right\rvert ^2+\left\lvert \begin{array}{cc} \frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} \\ \frac{\partial z}{\partial u} & \frac{\partial z}{\partial v} \\ \end{array} \right\rvert ^2}\;du\;dv$

(see section 7.4, Vector Calculus, 3ed. Jerrold E. Marsden and Anthony J. Tromba.)

(The use of u and v here instead of $\phi$ and $\theta$ is just to avoid implicitly suggesting that the formula has something to do with angles.)

See Figure 7c for more details about the regions of integration for our parametric form.

Results

Figure 8: Test results.

Figure 8 shows the results of running calculations of the first formula presented (Main), integration in Cartesian coordinates (Cart.Coords), integration in spherical coordinates (Sph.Coords), integration via Monte-Carlo (M.Carlo), and via the formula found in Tovchigrechko and Vakser (TV) with the added radius aspect, and these are for various angles of separation for radius equal to 1/3, 1, and 3, along with a1 and a2 values. The nan results in the Monte-Carlo are likely just caused by not having enough sample points fall within the tiny region, but the Monte-Carlo was just there as a sanity check, so not a big deal if we’re missing some values.

Appendix

Code

import numpy as np
from numpy import cos, arcsin, arccos, sin, tan, arctan, arctan2
from numpy import sqrt, pi, array, dot
from scipy import integrate, optimize
from math import isinf

# Note: I don't follow a naming convention here. Also, no attempt at
# being efficient was made. The purpose is to use alternative approaches
# to find discrepancies with the main function, i.e., AreaOfCapsIntersection.

# Here is our main process for computing the area of 
# intersecting spherical caps. a1 and a2 are the straight line
# distance of the radius of each cap's respective base. r is the sphere's
# radius. theta is the angular separation between the caps (in radians).
def AreaOfCapsIntersection(a1, a2, r, theta):
    a = arcsin(a1 / r)
    b = arcsin(a2 / r)
    c = theta
    s = (a + b + c) / 2
    k = sqrt((sin(s - a) * sin(s - b) * sin(s - c) / sin(s)))
    A = 2 * arctan(k / sin(s - a))
    B = 2 * arctan(k / sin(s - b))
    C = 2 * arctan(k / sin(s - c))
    T1 = r**2 * (A + B + C - pi)
    T2 = T1
    S1 = 2 * B * r**2 * (1 - cos(a))
    S2 = 2 * A * r**2 * (1 - cos(b))
    return S1 + S2 - (T1 + T2)

# This is from the Tovchigrechko and Vakser paper mentioned in another entry
# and modified to account for a radius other than 1. a1 and a2 have been
# adjusted to the distance on the sphere's surface by the calling code.
# Theta is also expected to be the angular separation of the 2 caps.
def TovchigrechkoAndVakser(a1, a2, r, theta):
    def csc(t):
        return 1 / sin(t)
    def cot(t):
        return 1 / tan(t)
    return r**2 * 2 * (pi 
    - arccos(cos(theta) * csc(a1) * csc(a2) - cot(a1) * cot(a2))
    - arccos(cos(a2) * csc(theta) * csc(a1) - cot(theta) * cot(a1)) * cos(a1)
    - arccos(cos(a1) * csc(theta) * csc(a2) - cot(theta) * cot(a2)) * cos(a2))

# This will return the distance from the origin to the base of the
# spherical cap given its baseRadius and the sphere's radius.
def FindCapDistance(baseRadius, r):
    return abs(optimize.brentq(lambda h: 
        sqrt(abs(baseRadius)**2 + abs(h)**2) - r, 0, 1000))

# Rotate a point about the origin
def Rotate(theta, x, y):
    c, s = cos(theta), sin(theta)
    R = [[c, -s], 
         [s, c]]
    pt = array([x, y]).reshape(2, 1)
    return dot(R, pt).flatten()

# Find the angle that rotates [xFrom, yFrom] onto [xTo, yTo].
def FindRotation(xFrom, yFrom, xTo, yTo):
    u0, u1, v0, v1 = xFrom, yFrom, xTo, yTo
    return arctan2(u0 * v1 - u1 * v0, u0 * v0 + u1 * v1)

# Rescales x along the interval [a,b] to the interval [c,d].
def Rescale(x, a, b, c, d):
    return (-(b * c) + a * d + (c - d) * x) / (a - b)

# Given the vectors a=fromPt and b=toPt from the origin find a point
# along the vector a+t*(b-a) at parameter t where t will be rescaled
# as [fromPt[0], toPt[0]] to [0,1].
def PointAlongLine(fromPt, toPt, t):
    a, b = array(fromPt), array(toPt)
    return a + Rescale(t, a[0], b[0], 0, 1) * (b - a)

# Finds the point where two lines intersect. Each parameter is a 2D array.
def FindIntersectionOfLines(l1From, l1To, l2From, l2To):
    def f(x):
        pt1 = PointAlongLine(l1From, l1To, x[0])
        pt2 = PointAlongLine(l2From, l2To, x[0])
        # The 2-norm is the same thing as Euclidean distance in our context:
        return np.linalg.norm(pt1 - pt2)
    # Because we're rescaling t in a+t(b-a) to coincide with the bounds
    # both functions are evaluated at the same x, so we just need 1 dimension.
    x0 = array([0])
    # The function is well behaved, keep it simple, nelder-mead will be fine.
    r = optimize.minimize(f, x0, method='nelder-mead', options={'xatol': 1e-8})
    return PointAlongLine(l1From, l1To, r.x[0])

# Integrate in Cartesian coordiantes given the radius, r, the
# Cartesian 2D region start point (p11, p12), the intersecting 
# point (p21, p22), and the end point (p31, p32). (See article for graphics.)
def IntegrateInCartesianCoordinates(r, p11, p12, p21, p22, p31, p32):
    r2 = r**2
    def commonIntegrand(a, b, x):
        d = sqrt(r2 - x**2)
        return r * (arcsin(b / d) - arcsin(a / d))
    def f(x):
        return PointAlongLine([p11, p12], [p21, p22], x)[1]
    def g(x):
        return PointAlongLine([p21, p22], [p31, p32], x)[1]
    areaType = 1 if p31 > p21 else (2 if p31 == p21 else 3)
    if areaType == 1:
        i1 = integrate.quad(lambda x: 
            commonIntegrand(f(x), sqrt(r2 - x**2), x), p11, p21)[0]
        i2 = integrate.quad(lambda x: 
            commonIntegrand(g(x), sqrt(r2 - x**2), x), p21, p31)[0]
        return 2 * (i1 + i2)
    elif areaType == 2:      
        return 2 * integrate.quad(lambda x: 
            commonIntegrand(f(x), sqrt(r2 - x**2), x), p11, p21)[0]
    elif areaType == 3:
        def h(x):
            # p31 may look out of place; it's the projection onto the horz line
            return PointAlongLine([p31, p22], [p21, p22], x)[1]
        def i(x):
            return PointAlongLine([p31, p32], [p21, p22], x)[1]        
        i1 = integrate.quad(lambda x: 
            commonIntegrand(f(x), sqrt(r2 - x**2), x), p11, p31)[0]
        i2 = integrate.quad(lambda x: 
            commonIntegrand(h(x), i(x), x), p31, p21)[0]
        return 2 * (i1 + i2)
    else:
        raise Exception("Area type not recognized.")
      
# Monte-Carlo integration in Cartesian coordiantes given the radius, r, the
# Cartesian 2D region start point (p11, p12), the intersecting 
# point (p21, p22), and the end point (p31, p32). (See article for graphics.)
def MonteCarloIntegration(r, p11, p12, p21, p22, p31, p32):
    sampleSize = 100000
    r2 = r**2
    def cross(p1, p2):
        return p1[0] * p2[1] - p2[0] * p1[1]
    def leftOf(p1, p2, p3): # is point p3 left of line [p1, p2]
        a1, a2, a3 = map(array, [p1, p2, p3])
        return cross(a3 - a2, a2 - a1) <= 0 # keep on the line too i.e. == 0
    def f(x):
        return PointAlongLine([p11, p12], [p21, p22], x)[1]
    def g(x):
        return PointAlongLine([p21, p22], [p31, p32], x)[1]
    def h(x):
        # p31 may look out of place; it's the projection onto the horz. line.
        return PointAlongLine([p31, p22], [p21, p22], x)[1]
    def i(x):
        return PointAlongLine([p31, p32], [p21, p22], x)[1]
    def interval1(x):
        return [f(x), sqrt(r2 - x**2)]
    def interval2(x):
        return [g(x), sqrt(r2 - x**2)]
    def interval3(x):
        return [h(x), i(x)]
    # Given line1 [p1, p2] and line2 [p2, p3] get random points within the
    # disk of radius r and to the left of l1 and l2. Note, p2 is the end
    # point of l1 and the start point of l2.
    def mcIntegrand(p1, p2, p3):
        u1 = np.random.uniform(-r, r, sampleSize)
        u2 = np.random.uniform(-r, r, sampleSize)
        pts = []
        for i in zip(u1, u2):
            x2 = i[0]**2
            y2 = i[1]**2
            if x2 + y2 < r2 \
                and leftOf(p1, p2, [i[0], i[1]]) \
                and leftOf(p2, p3, [i[0], i[1]]):
                pts.append(r / sqrt(r2 - i[0]**2 - i[1]**2))
        return np.mean(list(pts))
    areaType = 1 if p31 > p21 else (2 if p31 == p21 else 3)
    if areaType == 1:
        i1 = integrate.nquad(lambda y, x: 1, [interval1, [p11, p21]])[0]
        i2 = integrate.nquad(lambda y, x: 1, [interval2, [p21, p31]])[0]
        area = i1 + i2
        return 2 * area * mcIntegrand([p11, p12], [p21, p22], [p31, p32])
    elif areaType == 2:      
        area = integrate.nquad(lambda y, x: 1, [interval1, [p11, p21]])[0]
        return 2 * area * mcIntegrand([p11, p12], [p21, p22], [p31, p32])
    elif areaType == 3:
        i1 = integrate.nquad(lambda y, x: 1, [interval1, [p11, p31]])[0]
        i2 = integrate.nquad(lambda y, x: 1, [interval3, [p31, p21]])[0]
        area = i1 + i2
        return 2 * area * mcIntegrand([p11, p12], [p21, p22], [p31, p32])
    else:
        raise Exception("Area type not recognized.")
        
# Integration in spherical coordiantes given the radius, r, the
# Cartesian 2D region start point (p11, p12), the intersecting 
# point (p21, p22), and the end point (p31, p32). (See article for graphics.)
# Too slow; should be rewritten.
def SphericalCoordinateIntegration(r, p11, p12, p21, p22, p31, p32):
    # Just in case these names used by sympy may collide with the outer
    # code, they'll be scoped here as a precaution.
    from sympy import sin, cos, Matrix, atan, pi, sqrt
    from sympy.abc import rho, phi, theta, m, b
    from sympy import lambdify
    
    y = b + m * rho * cos(phi)
    fphi = atan(sqrt(-(y)**2 + rho**2 * (-(cos(phi)**2)) + rho**2) / y) / pi
    
    xyz = Matrix([rho * cos(phi),
                  rho * sin(phi) * cos(theta * fphi),
                  rho * sin(phi) * sin(theta * fphi)])
    uv = Matrix([phi, theta])
    
    mat = xyz.jacobian(uv)
    
    f = lambdify("phi, theta, rho, m, b", 
                 Matrix([[mat[0,0], mat[0,1]], [mat[1,0], mat[1,1]]]).det()**2)
    g = lambdify("phi, theta, rho, m, b", 
                 Matrix([[mat[1,0], mat[1,1]], [mat[2,0], mat[2,1]]]).det()**2)
    h = lambdify("phi, theta, rho, m, b", 
                 Matrix([[mat[0,0], mat[0,1]], [mat[2,0], mat[2,1]]]).det()**2)
    
    t1 = FindRotation(1, 0, p11, p12) # pt is on the circle
    t2 = FindRotation(1, 0, p21, np.sqrt(r**2 - p21**2)) # project onto circle
    t3 = FindRotation(1, 0, p31, p32) # pt is on the circle
    
    try: slope1 = (p22 - p12) / (p21 - p11)
    except ZeroDivisionError: slope1 = float('+inf')
    b1 = PointAlongLine([p11, p12], [p21, p22], 0)[1]
    try: slope2 = (p32 - p22) / (p31 - p21)
    except ZeroDivisionError: slope2 = float('+inf')
    b2 = PointAlongLine([p21, p22], [p31, p32], 0)[1]
    
    # The integrals are negated b/c the outer integral is done clockwise.
    i1 = -integrate.nquad(lambda p,t: np.sqrt(f(p, t, r, slope1, b1) 
        + g(p, t, r, slope1, b1) + h(p, t, r, slope1, b1))
        , [[t1, t2], [-pi, pi]])[0]    
    if isinf(slope2):
        return i1
    else:
        i2 = -integrate.nquad(lambda p,t: np.sqrt(f(p, t, r, slope2, b2) 
            + g(p, t, r, slope2, b2) + h(p, t, r, slope2, b2))
            , [[t2, t3], [-pi, pi]])[0]
        return i1 + i2

################################################
## Check the functions
################################################

def WithinTol(value):
    assert abs(value) <= 1e-8
    
WithinTol(abs(FindCapDistance(1, 1)))
WithinTol(abs(FindCapDistance(0.7453559924999299, 1) - 2/3))
assert Rescale(15, 10, 20, 0, 1) == 0.5
assert PointAlongLine([1,2], [11,3], 6)[1] == 2.5
WithinTol(abs(FindRotation(1, 0, -0.5, 0.86602540378443) - 2 * pi / 3))
WithinTol(abs(Rotate(2 * pi / 3, 1, 0)[0] - -0.5))
WithinTol(abs(Rotate(2 * pi / 3, 1, 0)[1] - 0.8660254037844386))
_r = FindIntersectionOfLines([5, 1], [1, 10], [1, 2], [10, 5])
WithinTol(abs(_r[0] - 4.09677419))
WithinTol(abs(_r[1] - 3.03225807))
print("Tests passed.")

################################################
## Driver
################################################

print("MonteCarloIntegration is the slowest function here.")
print("Disable it or set its sample size smaller to speed things up.")
print(".....")
print("SphericalCoordinateIntegration is also slow, but would need to")
print("be rewritten to account for its slowness.")
print(".....")
print("The calls to MonteCarloIntegration and SphericalCoordinateIntegration")
print("have been commented out of the Driver.")
print(".....")

# Each row holds a value respectively for these elements:
# [radius, a1, a2].
rows = [
    [1/3, 0.235702, 0.321975],
    [1, 0.707107, 0.965926],
    [3, 2.12132, 2.89778]]

for row in rows:
    [r, a1, a2] = row
    p11 = [-a1, FindCapDistance(a1, r)]
    p12 = [a1, p11[1]]
    p21 = [-a2, FindCapDistance(a2, r)]
    p22 = [a2, p21[1]]
    # Calculate the angle step size to sweep across:
#    fromAngle = FindRotation(1, 0, p12[0], p12[1]) \
#        - FindRotation(1, 0, p22[0], p22[1]) + 0.000174533
#    toAngle = FindRotation(1, 0, p11[0], p11[1]) \
#        - FindRotation(1, 0, p22[0], p22[1])- 0.000174533
#    stepSize = (toAngle - fromAngle) / 4
#    atAngle = fromAngle
    # I am inserting a pi/2 angle to the above angles to try a vertical line.
    angles = [0.5237767175432327,0.916388451764179,1.3090001859851252,
              1.5707963267948966,1.7016119202060715,2.094223654427018]
    # Generate the report for the current arguments:
    #while atAngle <= toAngle:
    for atAngle in angles:
        print([r, a1, a2, atAngle])
        rp21 = Rotate(atAngle, p21[0], p21[1])
        rp22 = Rotate(atAngle, p22[0], p22[1])
        intersection = FindIntersectionOfLines(p11, p12, rp21, rp22)
        x11, x12 = p11[0], p11[1]
        x21, x22 = intersection[0], intersection[1]
        x31, x32 = rp22[0], rp22[1]
        print(AreaOfCapsIntersection(a1, a2, r, atAngle))
        print(IntegrateInCartesianCoordinates(r, x11, x12, x21, x22, x31, x32))
        #print(SphericalCoordinateIntegration(r, x11, x12, x21, x22, x31, x32))
        #print(MonteCarloIntegration(r, x11, x12, x21, x22, x31, x32))
        print(TovchigrechkoAndVakser(arcsin(a1 / r), arcsin(a2 / r), r, atAngle))
        #atAngle += stepSize