Calcular el área en una esfera de la intersección de dos tapas esféricas

Question

Dada una esfera de radio $r$ con dos casquetes esféricos en ella definidos por los radios ($a_1$ y $a_2$) de las bases de los casquetes esféricos, dado un separación de los dos casquetes esféricos por el ángulo $\theta$, ¿cómo se calcula el área de la intersección?

Para aclarar, el área es la de la superficie curva en la esfera definida por la intersección. En el caso extremo donde ambos $a_1,a_2 = r$, estaríamos describiendo un casquete esférico.

Alternativamente, define los casquetes esféricos por los ángulos $\Phi_1 = \arcsin(a_1/r)$ y $\Phi_2 = \arcsin(a_2/r)$.

Answer 1

Probablemente ya no uses el sitio ya que la pregunta se realizó hace un año, pero de todas formas publicaré una respuesta porque era un buen problema y hay una leve posibilidad de que alguien más lo vea.

Antes de resolver esto, permíteme aclarar algunas nomenclaturas/suposiciones:

1. En lugar de $\theta$ como el ángulo separador de las dos capas esféricas, sea $\alpha$.
2. Cuando me refiera al eje central de una capa esférica (SC) implica el eje que intersecta el centro de la base del SC y el centro de nuestra esfera.
3. Supondré que ninguna de las SCs es más grande que hemisferios.
4. Supondré que las SCs se intersectan.

Por lo tanto, lo más lejos que puede estar un punto de la intersección de las dos SCs desde el centro de la SC solo puede ser $\Delta\phi = \pi/2$, usando el sistema de coordenadas esféricas.

El problema establece que hay dos SCs con radios de base $a_1$ y $a_2$ en una esfera de radio $r$ con un ángulo de $\alpha$ entre los ejes centrales de las dos SCs.

Como mencionaste, el ángulo entre el eje central de la SC y un radio desde el centro de la esfera hasta la circunferencia de la base del SC son $\Phi_1 = \arcsin(a_1/r)$ y $\Phi_2 = \arcsin(a_2/r)$.

Alinea la esfera de modo que el eje central de la SC-1 y el eje $x$ sean paralelos (en el sistema de coordenadas $\mathbb{R^3$}). Traduce la esfera al origen. Rota la esfera de modo que el eje central de la SC-2 esté sobre el plano $xz$. El ángulo $\alpha$ ahora sube desde el eje $x$ hacia el eje $z$. El ángulo $\Phi_2$ está contenido en $\alpha$ en el mismo plano, subiendo en la misma dirección (pero no comenzando desde el eje $x$). Estas acciones son legales ya que no distorsionan la forma ni su área superficial.

Ahora que hemos definido algo rígidamente nuestro problema, podemos abordarlo.

Permite que la curva espacial que define la intersección de la base de la SC-2 con la esfera sea $\vec s_2(t)$. Después de algo (léase: mucho) de aritmética vectorial y trigonometría, llegamos a $\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\} $$ para $0\le t\le 2\pi$. Esta expresión de $\vec s_2(t)$ nos permite proyectar en el plano $xy$. La proyección de esta intersección base-esfera se ve como una elipse. La parametrización de la curva espacial $\vec s_1(t)$ es un poco más fácil debido a la forma en que configuramos nuestros ejes. $$ \vec s_1(t)=\left\{\sqrt{r^2-a_1^2},a_1\sin t, a_1\cos t\right\} $$ Cuando se proyecta en el eje $xy$ esto se ve como un segmento de línea. Usando las dos curvas espaciales proyectadas, puedo encontrar los límites angulares para $\theta$ en el sistema de coordenadas esféricas para 3D encontrando la intersección de la elipse y el segmento de línea. $$ \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} $$ ¡Es confuso, pero al menos es una constante! Sabemos que para el área superficial que nos interesa son los valores de $\theta$ tales que $-\theta_0\le \theta\le\theta_0$

Ahora necesitamos encontrar los límites angulares para $\phi$ en nuestro sistema de coordenadas esféricas. Sabemos que estos cambian sobre $\theta$ ya que puedes imaginar que a medida que giramos nuestra esfera a lo largo del eje $z$, el ancho del área de interés cambia (el ángulo desde el eje $z$ a los vectores de posición del área superficial que estamos mirando cambia).

Para el límite superior de $\phi$, observamos la superficie de la capa esférica inferior a lo largo de la esfera. Usamos la parametrización que proporcioné anteriormente. Dado que $z_1(t) = a_1cos(t)$ y $t$ es análogo a $\theta$ en mi parametrización, puedo usar $\cos \phi = \frac{z}{r}$ para encontrar que $$ \phi_1(\theta) = \arccos\frac{a_1\cos\theta}{r} $$ De manera similar, usando $$ z_2(t)=r\sin(\alpha-\Phi_2)+a_2(1-\cos t)\sin(\arcsin(\frac{r}{a_2}\sin\Phi_2)+\alpha) $$ lo cual no posteé en $\vec s_2(t)$ debido a preocupaciones de espacio, puedo encontrar $$ \phi_1(\theta) = \arccos\frac{z_2(\theta)}{r} $$ Así, para el área superficial $A$ de la intersección de dos capas esféricas $\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 $$ Seré el primero en admitir que esta integral probablemente solo es soluble numéricamente, pero no pude encontrar enfoques geométricos elegantes para esto, así que opté por el método "fuerza bruta".

EDICIÓN

Gracias a la respuesta de Christian Blatter a esta pregunta, puedo responder de una manera más concisa.

Si tomo nuestra esfera situada en el origen como se describió anteriormente, y la roto en la dirección positiva sobre $\hat j$ hasta que $\alpha$ esté centrada sobre el eje $z$, entonces digo:

1. Proyectar curvas espaciales $s_1(t)$ y $s_2(t)$ en el plano $xy$.
2. Dichas curvas deben ser elipses o círculos, ya que el caso extremo $\alpha=\pi$ es el único que produce líneas rectas al proyectar y la respuesta a ese escenario es que el área superficial $A=0$.

Si una elipse está contenida dentro de otra, entonces el área superficial se puede encontrar usando la fórmula para encontrar el área superficial externa de una capa esférica (ver respuesta de Christian).

Si las elipses se intersectan en dos puntos (no pueden intersectarse en este caso en tres o cuatro), entonces usando la respuesta a la pregunta mencionada anteriormente, podemos configurar una integral superficial para el cálculo de $A$, dado que la intersección de las elipses en el plano $xy$ se puede describir por las dos desigualdades $a\le u\le b$ y $g(u)\le v\le f(u)$ de modo que $\left\{a,b:a,b\in\mathbb{R}\land ag(u)$ sobre $u:[a, b]$.

Sea $x=\frac{u}{\sqrt{r_x}}+c_x$ y $y=v$, para $r_x$ siendo el semidiámetro en la dirección $x$ y $c_x$ siendo el desplazamiento horizontal desde el origen al centro de la elipse de la proyección de $s_2(t)$ en el plano $xy$.

Sea $\vec F=\{x, y,\sqrt{r^2-x^2-y^2}\}$, el vector de posición. $$ 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 $$

Answer 2

Este problema puede ser resuelto por geometría elemental, respectivamente, trigonometría esférica.

Podemos asumir $r=1$. Deje que $\rho$ y $\rho'$ sean los radios esféricos de las dos tapas.

Cuando las tapas están disjuntas no hay nada que calcular.

Cuando la tapa más pequeña (de radio $\rho$) es un subconjunto de la más grande, entonces el área $A_\cap$ de la intersección está dada por $$A_\cap=4\pi\sin^2{\rho\over2}\ .$$

Ahora el caso interesante: Asuma que los dos círculos límite se intersecan en dos puntos $P$ y $Q$. Conecte $P$ y $Q$ por un arco de un gran círculo; además conecte los centros $C$ y $C'$ de las dos tapas tanto con $P$ como con $Q$. Los dos arcos que emanan de $C$ tienen una longitud de $\rho$ y encierran un ángulo $\alpha>0$, y de manera similar hay un ángulo $\alpha'$ en $C'$. Los dos ángulos $\alpha$ y $\alpha'$ que aparecen en $C$ y $C'$ pueden ser calculados a partir de los datos dados por trigonometría esférica. Dejo eso al OP.

Entonces verá que la intersección de las dos tapas (también en casos no cubiertos por la figura anterior) puede ser escrita como una suma algebraica (es decir, $\pm$) de formas de los siguientes tipos:

(a) Sectores $S$ de tapas esféricas de radio $\rho$, respectivamente $\rho'$, y ángulo central $\alpha>0$, respectivamente $\alpha'>0$. El área de dicho sector está dada por $$A_S=2\alpha\sin^2{\rho\over2}\ .$$ (b) Triángulos esféricos isósceles $\Delta$ con dos lados iguales a $\rho$ (o $\rho'$) y ángulo central $\alpha>0$. El ángulo de base $\beta$ de dicho triángulo puede ser calculado a partir de $$\tan\beta={\cot(\alpha/2)\over\cos\rho}\ ;$$ entonces su área $A_\Delta$ está dada por $$A_\Delta=\alpha+2\beta-\pi\ .$$

Answer 3

Aquí tienes una fórmula simplificada como una función de tus 3 variables, $a_1$, $a_2$ y $\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) $$

Como se mencionó anteriormente, las tapas deben intersectar y una no puede contener completamente a la otra.

Esta solución está copiada de una presentación gráfica de AMD, y es originalmente de un artículo de bioquímica. No tengo pruebas de que funcione, así que tómalo por lo que vale.

TOVCHIGRECHKO, A. AND VAKSER, I.A. 2001. ¿Qué tan común es el paisaje de energía en forma de embudo en las interacciones proteína-proteína? Protein Sci. 10:1572-1583

Answer 4

Introducción

Hace poco me llamaron la atención sobre este problema y me pidieron que hiciera una contribución. Responderé a esta pregunta de 4 maneras diferentes. También voy a comparar numéricamente las soluciones a través de código Python, y voy a proporcionar algunos gráficos 2D y 3D para ayudar a ilustrar los conceptos.

El artículo de Tovchigrechko y Vakser (TV) (mencionado y enlazado en otra entrada) describe la solución geométricamente y es similar a la solución de Christian Blatter, y es similar al enfoque de trigonometría esférica mostrado en esta entrada. Una nota sobre el artículo de TV: los radios de los casquetes se miden a lo largo de la superficie de la esfera, a diferencia de lo que pide la OP, que son distancias rectilíneas limitadas al plano de la base del casquete.

Solución

Estos son los pasos para resolver el problema de intersección (todos los ángulos están en radianes a lo largo de esta entrada; para más contexto sobre los dados ver el PO y la Figura 1 (r es el radio, a1 y a2 son distancias en línea recta, y $\theta$ es una separación angular)). Supondremos que cada casquete tiene la semiesfera como caso límite; también, como se ha señalado en otras entradas, supondremos que se produce una intersección en la que ninguno de los casquetes está completamente contenido dentro del otro:

Figura 1: Gráfico 3D que muestra los datos.

$\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) } $$

Esa es la respuesta (ésta, y diferentes fórmulas de ángulos esféricos, también pueden utilizarse para derivar variantes de la fórmula, incluida la que se encuentra en Tovchigrechko y Vakser). El resto de esta entrada sólo ofrece más contexto y explicaciones.

Más información

Lo siguiente es equivalente a los pasos anteriores, y es el resultado que se encuentra en el documento de TV excepto el $r^2$ parte (a1, a2, y $\theta$ son distancias esféricas, es decir, distancias a lo largo de la superficie curva en la esfera unidad; para que a1 y a2 de arriba funcionen aquí, conviértelas así: t=arcsin(a1/r) y 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}) )$

Figura 2: Dos ejemplos de intersección.

Las figuras 2a y 2b muestran una intersección de 2 casquetes esféricos. El Plano I y el Plano II son las bases planas de sus respectivos casquetes; los planos se utilizan para ayudar a resaltar la intersección, y se utilizarán más adelante para convertir el problema a 2D. Utilizaremos la figura 2a para la discusión siguiente.

Figura 3: Sectores esféricos.

La figura 3 muestra los 2 sectores esféricos separados, es decir, un área roja y un área cian. La suma de las áreas roja y cian contiene la intersección de los casquetes esféricos, pero observe los 2 triángulos esféricos que se incluyen dos veces al sumar el área roja y el área cian. El área de estos 2 triángulos esféricos debe eliminarse del área combinada de los 2 sectores. (En este contexto a y b tienen un componente de radio que hay que dividir. En la fórmula de intersección anterior, el componente de radio se elimina en el momento en que se calcula el área del sector, por lo que las fórmulas de la figura 3 parecen ligeramente diferentes en comparación con las fórmulas de la figura 4. $\mathcal{A}_{\text{s1}}$ y $\mathcal{A}_{\text{s2}}$ en la fórmula de intersección).

Figura 4: Primer plano de un triángulo esférico.

La figura 4 muestra lo que hay que hacer para calcular el área de los 2 triángulos esféricos cuya área combinada hay que quitar del área combinada de los 2 sectores esféricos.

Los ángulos del triángulo esférico se han calculado utilizando las fórmulas de medio ángulo para tangentes (ecuaciones 27 y 34-38):

/Trigonometría esférica

Y la fórmula del área del triángulo esférico aparece aquí:

/TriánguloEsférico

Ambos temas están en Wolfram's Mathworld, pero no puedo hacer estos enlaces porque mi cuenta está limitada a sólo 8 enlaces en este momento.

Verificación numérica (o soluciones adicionales)

Coordenadas cartesianas

A continuación, compararemos la fórmula con la integración numérica sobre una región. Imagina que giras los objetos de la figura 2 de tal manera que miras hacia la línea de intersección de los planos. Los planos aparecerían entonces como líneas. Aquí tienes una versión en 2D (las líneas rojas sugieren la continuación de los planos más allá del círculo):

Figura 5: Perspectiva 2D de una intersección.

En la figura 5, a1 y a2 son los radios de sus respectivas tapas que se muestran aquí para contextualizar (y el radio, r, se muestra para evitar confusiones con $\mathcal{R}$ de la región). La zona azul es la región proyectada en 2D donde se cruzan las tapas. Giraremos (es decir, cambiaremos $\theta$ (un determinado) un tapón a través de la luz del otro tapón para comprobar el área de varias intersecciones diferentes, como muestra la serie de imágenes de la figura 6a. En la 6b se muestran algunas perspectivas de las 2 tapas que se intersecan en relación con una de las regiones. Una de las razones para hacer este barrido es considerar los puntos extremos en los que podríamos encontrar inestabilidad numérica.

Figura 6: Integración en coordenadas cartesianas.

La figura 6c muestra los conceptos implicados en el cálculo de la integral de las tapas de intersección en coordenadas cartesianas. La integral doble, con respecto a la función $r / \sqrt{r^2-x^2-y^2}$ , es impropio y probablemente por esa razón presentó algunos problemas para scipy. Conseguí que funcionara cambiando las variables tanto en la versión integral doble como en la integral simple, pero en el código he elegido la versión integral simple porque hace que el código sea un poco más fácil de leer. Esto se corresponde con el conjunto de integrales bajo la línea horizontal discontinua de la Figura 6c, que aplica a la integral interior la técnica común de hacer la sustitución $y=\lvert b\rvert \sin (\theta )$ . En $\sin ^{-1}\left(\frac{\lvert \sqrt{r^2-x^2} \rvert}{\lvert \sqrt{r^2-x^2} \rvert}\right)$ puede reducirse a $\pi /2$ .

Integración Monte-Carlo

A continuación, se establecerá una aproximación sencilla de la integración con el fin de que sirva de barandilla para detectar errores evidentes que puedan estar presentes. En concreto, se utilizará una integración de Monte-Carlo. Aquí simplemente eliminamos el integrando de las integrales dobles de la Figura 6c e integramos la constante escalar de 1 sobre las mismas regiones como se muestra en la Figura 6d. Si confiamos en las definiciones de las regiones (es decir, si son sencillas) y si las funciones se comportan bien, entonces esto presenta un enfoque robusto para calcular el área; por estas razones es el enfoque que utilicé para calcular el área de la región necesaria en la función de Monte-Carlo en el ejemplo del código Python. No es necesario multiplicar por 2, ya que se trata del área de la región 2D, no de la superficie de la esfera situada por encima y por debajo de la región. El valor de esta área de región se multiplica por la media de los $r / \sqrt{r^2-x^2-y^2}$ evaluado en cada punto de un conjunto aleatorio de puntos dentro de la región, que produce la aproximación de Monte-Carlo de la superficie por encima de la región (multiplicar por 2 para incluir el área por debajo de la región).

Integración en coordenadas esféricas

La ventaja en coordenadas esféricas es que evitaremos una integral impropia en su mayor parte, pero una desventaja es que la fórmula del área parecerá más compleja. He aquí una forma paramétrica de un casquete esférico que podemos utilizar para calcular su área (se diseñó para esta discusión y no pretende ser lo suficientemente flexible para un uso general (el dominio de $\phi$ dependerá de la tapa, pero para este problema $\theta$ siempre se ejecuta desde $-\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$

Figura 7: Integración en coordenadas esféricas.

La figura 7a es una visualización que muestra cómo el punto en (1,0,0) se transforma en a/el arco rojo a lo largo del casquete esférico de acuerdo con la forma paramétrica dada. La multiplicación de las matrices de la figura 7a da como resultado la forma paramétrica anterior. La dirección $f(\phi )$ parte; aunque desordenada, es sólo calcular la intersección en el valor x de $\phi$ a lo largo de la línea de la base de un casquete esférico (viéndolo de canto), y proyectando luego esa intersección hacia arriba y hacia abajo sobre la esfera, lo que da lugar a puntos a lo largo del arco rojo de la figura 7a con respecto a $\theta$ como barre a través de $-\pi \text{ to } \pi$ entre esas 2 proyecciones, y finalmente barriéndolo a través de $\phi$ rinde la superficie de la tapa. Las imágenes de la Figura 7b se calcularon mediante la forma paramétrica anterior y sirven como verificación visual de que la forma paramétrica es correcta.

El área de una superficie parametrizada es igual a la magnitud de las normales infinitesimales del plano tangente integradas sobre la región, es decir:

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

y en nuestro contexto de tener x , y y z dado en términos de dos parámetros es equivalente a (así es como se hace en el código Python del apéndice):

$\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$

(véase el apartado 7.4, Cálculo Vectorial, 3ed. Jerrold E. Marsden y Anthony J. Tromba).

(El uso de u y v en lugar de $\phi$ y $\theta$ es sólo para evitar sugerir implícitamente que la fórmula tiene algo que ver con los ángulos).

Véase la Figura 7c para más detalles sobre las regiones de integración de nuestra forma paramétrica.

Resultados

Figura 8: Resultados de las pruebas.

La figura 8 muestra los resultados de la ejecución de los cálculos de la primera fórmula presentada (Principal), la integración en coordenadas cartesianas (Cart.Coords), la integración en coordenadas esféricas (Sph.Coords), la integración vía Monte-Carlo (M.Carlo), y vía la fórmula encontrada en Tovchigrechko y Vakser (TV) con el aspecto del radio añadido, y estos son para varios ángulos de separación para radios iguales a 1/3, 1, y 3, junto con los valores a1 y a2. La dirección nan en el Monte-Carlo probablemente se deban a que no hay suficientes puntos de muestra dentro de la pequeña región, pero el Monte-Carlo sólo estaba ahí para comprobar la cordura, así que no pasa nada si nos faltan algunos valores.

Anexo

Código

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

Answer 5

El penúltimo signo negativo en la fórmula anterior debería ser un signo más. La fórmula en el papel de Tovchigrechko y Vakser es correcta. El error con el signo aquí proviene de la presentación citada anteriormente.