How to calculate the surface area of a sphere cutting off two intersecting perpendicular caps? For a theoretical design of a football helmet, I am modeling it as a sphere with radius r, with a bottom cap removed ($-y$ direction) (distance between center of helmet and center of the cutoff cross section is $ar$) for the neck. There is also the facemask area, which I am modeling by cutting off a cap pointed toward $+x$ direction whose center is a distance $br$ from the helmet center. Clearly, $0< a, b < 1$.
Since I don't have a picture, I'll provide an extreme case of what I'm looking for. If $a = b = 0$, then the $+y$, $-x$ hemisphere of the sphere will be left, and the SA in that case would be $\frac14 4\pi r^2$. 
Here is my attempt so far:
I can find the SA of the entire sphere and the individual caps, but I cannot subtract because the shared SA between the caps. I thought about looking at a surface integral, but I wasn't able to figure out the proper bounds. I would love an analytic solution, but I can solve the integral numerically if need be.
 A: For simplicity, we will only consider the case $r = 1$.
I have changed the association of parameters $a,b$ from directions $-y$ and $+x$ to directions $+x$ and $+y$. The end formula is the same as it is symmetric with respect to $a$ and $b$.
We will use following version of Gauss Bonnet theorem to compute the desired area.

For any "polygonal" region $\Omega \subset S^1$ whose boundary $\partial \Omega$ consists of vertices $v_1, v_2, \ldots, v_n$ joined by smooth curves $e_1, e_2, \ldots, e_n$ as edges ($e_i$ connects $v_i$ to $v_{i+1}$, $v_{n+1}$ is an alias of $v_1$). Its area is given by the formula:
  $$\verb/Area/(\Omega) = 2\pi - \sum_{i=1}^n \Delta(v_i) - \sum_{i=1}^n \int_{e_i} k_g ds$$
  where $\Delta(v_i)$ is the angle between the tangent vectors of $e_{i-1}$ and $e_{i}$ at vertex $v_i$ ($e_0$ is an alias of $e_n$). $k_g$ is the geodesic curvature along edge $e_i$.

For any $a,b \in \mathbb{R}$, let $a' = \frac{a}{\sqrt{1-a^2}}$, $b' = \frac{b}{\sqrt{1-b^2}}$ and $c = \sqrt{a^2+b^2}$. 
Let $X_a$ and $Y_b$ be the half spaces $x \ge a$ and $y \ge b$ respectively. 
When $a^2+b^2 < 1$, $a', b', c$ are real numbers and the intersection of above two halfspaces with unit sphere, $\Omega = S^1 \cap X_a \cap Y_b$, is non-empty. The boundary $\partial \Omega$ consists of two vertices $v_1 = (a,b,c)$, $v_2 = (a,b,-c)$ and having two circular arcs $e_1$, $e_2$ as edges.


*

*The circular arc $e_1$ connects $v_1$ to $v_2$. It lives on small circle $S^1 \cap \partial Y_b$ on the plane $y = b$. It subtends an angle $\pi - 2\tan^{-1}\frac{a}{c}$ with respect to its center $(0,b,0)$. Since $\int k_g ds$ over the whole small circle $S^1 \cap \partial Y_b$ is $2\pi b$, we find  $$\int_{e_1} k_g ds = b\left(\pi - 2\tan^{-1}\frac{a}{c}\right)$$

*The circular arc $e_2$ connects $v_2$ to $v_1$. It lives on small circle $S^1 \cap \partial X_a$ one the plane $x = a$. It subtends an angle $\pi - 2\tan^{-1}\frac{b}{c}$ with respect to its center $(a,0,0)$. Since $\int k_g ds$ over the whole small circle $S^1 \cap \partial X_a$ is $2\pi a$, we find  $$\int_{e_2} k_g ds = a\left(\pi - 2\tan^{-1}\frac{b}{c}\right)$$

*$e_1$ lies on plane $y = b$, its tangent vector at $v_1$ points toward $(0,1,0) \times (a,b,c) = (c,0,-a)$.
$e_2$ lies on plane $x = a$, its tangent vector at $v_1$ points toward $(1,0,0) \times (a,b,c) = (0,-c,b)$.
This means at $v_1$, the tangent vectors has turned for an angle $$\Delta(v_1) = \cos^{-1}\left[\frac{(c,0,-a)\cdot(0,-c,b)}{\sqrt{(c^2+a^2)(c^2+b^2)}}\right] = 
\cos^{-1}\left[-\frac{ab}{\sqrt{(1-a^2)(1-b^2)}}\right] =
\pi - \cos^{-1}(a'b')$$

*By symmetry, we have $\Delta(v_2) = \Delta(v_1)$. 
Combine all these, we find
$$\begin{align}\verb/Area/(\Omega) 
&= 2\pi - 2(\pi - \cos^{-1}(a'b')) 
- a\left(\pi - 2\tan^{-1}\frac{b}{c}\right)
- b\left(\pi - 2\tan^{-1}\frac{a}{c}\right)\\
&= 2\left(\cos^{-1}(a'b') + a\tan^{-1}\frac{b}{c} + b\tan^{-1}\frac{a}{c}\right) - \pi(a+b)
\end{align}
$$
Notice if $\theta = \cos^{-1}(a'b')$, then
$$\cos^2\theta = (a'b')^2 = \frac{(ab)^2}{(1-a^2)(1-b^2)} = \frac{(ab)^2}{c^2+(ab)^2}
\implies \cot^2\theta = \left(\frac{ab}{c}\right)^2$$
With a little bit of algebra, one can verify $\cos^{-1}(a'b') = \frac{\pi}{2} - \tan^{-1}\frac{ab}{c}$.
We can simplify above expression of area to
$$\bbox[padding: 1em;border:1px solid blue;]{
\verb/Area/(\Omega) =
\pi(1-a-b) + 2\left(a\tan^{-1}\frac{b}{c} + b\tan^{-1}\frac{a}{c} - \tan^{-1}\frac{ab}{c}\right)
}$$
A: Great question. This is not a complete answer, but it's too long to be a comment. 
For convenience, I'll let $r=1$, since the area scales up by $r^2$. It's possible to compute the area of one of the caps minus the intersection.
Let's say it's the $x$ cap, then you have a cross section circle with radius $\sqrt{1-a^2}$ parallel to the $yz$ plane, therefore this partial cap is bounded by
$$ a < x < \sqrt{1-y^2-z^2}, \ -\sqrt{1-a^2-y^2} < z < \sqrt{1-a^2-y^2},\ -b < y < \sqrt{1-a^2} $$
If you define this surface as $x = f(y,z) = \sqrt{1-y^2-z^2}$ then the surface area is given by the integral
$$ \iint \sqrt{\left(\frac{\partial x}{\partial y}\right)^2 + \left(\frac{\partial x}{\partial z} \right)^2+1}\ dz\ dy = \int_{-b}^{\sqrt{1-a^2}}\int_{-\sqrt{1-a^2-y^2}}^{\sqrt{1-a^2-y^2}} \frac{1}{\sqrt{1-y^2-z^2}}\ dz\  dy $$
I don't know if this has an analytic solution, so that's up to you to find out.
Alternatively, you may also compute the intersectional area instead by integrating $y$ from $-\sqrt{1-a^2}$ to $-b$
