# Derivations for `constrsphtrig.py`: _Constructible spherical trigonometry in Python_

This document outlines the mathematical derivations underlying a number of algorithms in `constrsphtrig.py`.
We represent a point on the unit 2-sphere using spherical coordinates $p=(\phi_p,\theta_p)$, where $-\pi/2\leq\phi\leq\pi/2$ is its latitude and $-\pi<\theta\leq\pi$ its longitude.
$d_{ab}$ gives the distance between two points, and $\angle_{abc}$ the angle between three points.

For a spherical triangle with vertices $a,b,c$, the
[spherical cosine rule](https://en.wikipedia.org/wiki/Spherical_law_of_cosines)
gives

$$
    \cos d_{ab} = \cos d_{ac} \cos d_{bc} + \sin d_{ac} \sin d_{bc} \cos\angle _{acb} ~.
$$

Given any two points $a,b$ on the unit sphere, the
[geodesic distance](https://en.wikipedia.org/wiki/Great-circle_distance)
is

$$
    \cos d_{ab} = \sin\phi_a \sin\phi_b + \cos\phi_a \cos\phi_b \cos(\theta_a-\theta_b) ~,
$$

while the
[Cartesian coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system#Cartesian_coordinates)
of any point $p$ on the sphere are

$$
    \mathbf{p} = (x_p, y_p, z_p)
    = (\cos\phi_p \cos\theta_p, \cos\phi_p \sin\theta_p, \sin\phi_p) ~.
$$

## `ConstrSphPoint.ang(self, lon_e=zero)`: _Returns the angle between `self` and the equator at a given longitude._

Given a point $p$ and a point on the equator $q$, compute $\angle_{pqq'} := \angle_{(\phi_q,\theta_q) (0,\theta_q) (0,\theta_q+\pi/2)}$.
Assume that $p\not\in\{q,(0,\theta_q+\pi)\}$.
Since $d_{q'q}=\pi/2$, the spherical cosine rule yields

$$
    \cos\angle_{pqq'}
    = \frac{\cos d_{pq'}}{\sin d_{pq}} ~.
$$

The geodesic distance formula then gives the necessary data to obtain the numerator and denominator, noting that $\phi_q=\phi_{q'}=0$:

\begin{align}
    \cos d_{pq'} &= \cos\phi_p \sin(\theta_p-\theta_q) ~, \\
    \cos d_{pq} &= \cos\phi_p \cos(\theta_p-\theta_q) ~.
\end{align}

## `ConstrSphPoint.linear(a, b, d_ac, d_bc)`: _Solves the linear problem in spherical geometry._

_Given two points of a spherical triangle and distances from these points to
        the third point, this function returns the list of possible third points $c$._

Here, we solve for $x_c, y_c$ linearly in terms of $z_c$ and then apply the nonlinear condition $\mathbf{c}\cdot\mathbf{c}=1$.
The geodesic distance formula expressed in Cartesian coordinates succinctly yields $\cos d_{ab} = \mathbf{a}\cdot\mathbf{b} = \mathbf{a}^T \mathbf{b}$, so that

$$
\begin{pmatrix}
    \cos d_{ac} \\ \cos d_{bc}
\end{pmatrix}
= \begin{pmatrix}
    \mathbf{a}^T \\ \mathbf{b}^T
\end{pmatrix} \mathbf{c}
= \begin{pmatrix}
    x_a & y_a \\
    x_b & y_b
\end{pmatrix}
\begin{pmatrix}
    x_c \\ y_c
\end{pmatrix}
+
\begin{pmatrix}
    z_a \\ z_b
\end{pmatrix}z_c
~.
$$

Thus

$$
\begin{pmatrix}
    x_c \\ y_c
\end{pmatrix}
=
\begin{pmatrix}
    x_a & y_a \\
    x_b & y_b
\end{pmatrix}^{-1}
\left(
\begin{pmatrix}
    \cos d_{ac} \\ \cos d_{bc}
\end{pmatrix}
- \begin{pmatrix}
    z_a \\ z_b
\end{pmatrix}z_c
\right)
$$

Denoting

$$
    M :=
\begin{pmatrix}
    x_a & y_a \\
    x_b & y_b
\end{pmatrix}, ~
\mathbf{d'} := 
\begin{pmatrix}
    \cos d_{ac} \\ \cos d_{bc}
\end{pmatrix}, ~
\mathbf{z'} := 
\begin{pmatrix}
    z_a \\ z_b
\end{pmatrix}
$$

allows one to express the above as

$$\begin{equation}
\begin{pmatrix}
    x_c \\ y_c
\end{pmatrix}
=
\frac{1}{\det M}
\begin{pmatrix}
    y_b & -y_a \\
    -x_b & x_a
\end{pmatrix}^{-1}
(\mathbf{d'} - \mathbf{z'}z_c) ~. \tag{1}
\end{equation}$$

Applying the unit spherical condition $\mathbf{c}\cdot\mathbf{c}=1$ gives

$$
    \frac{1}{\det^2 M} \left(
        (y_b,-y_a)\cdot(\mathbf{d'}-\mathbf{z'}z_c)~^2 + 
        (-x_b,x_a)\cdot(\mathbf{d'}-\mathbf{z'}z_c)~^2
    \right)
    + z_c^2 = 1 ~,
$$

which is a quadratic equation in $z_c$.
Collecting its coefficients, this quadratic equation has the form $\alpha_2 z_c^2 - 2 \alpha_1 z_c + \alpha_0 = 0$, where

$$
\begin{align}
    \alpha_2 &:= \alpha_{yz}^2 + \alpha_{xz}^2 + \det\phantom{}^2 M \\
    \alpha_1 &:= \alpha_{yd}\alpha_{yz} + \alpha_{xd}\alpha_{xz} \\
    \alpha_0 &:= \alpha_{yd}^2 + \alpha_{xd}^2 - \det\phantom{}^2 M \\
    \alpha_{yu} &:= (y_b, -y_a) \cdot \mathbf{u'} \\
    \alpha_{xu} &:= (-x_b, x_a) \cdot \mathbf{u'}
\end{align}
$$

and $u\in\{d,z\}$.
The solution to this quadratic equation yields

$$
z_c = \frac{\alpha_1 \pm \sqrt{\alpha_1^2 - \alpha_2\alpha_0}}{\alpha_2} \tag{2} ~.
$$

Substituting $(2)$ into $(1)$ recovers the Cartesian coordinates of $c$.
Since the latitude always satisfies $\cos\phi\geq0$, it can be recovered directly from the last Cartesian component:
$$
    \phi_c = \verb|ConstrAngle|\left(\mathrm{sgn}~z_c, \sqrt{1-z_c^2}\right) ~.
$$

The longitude satisfies $\cos\theta=\frac{x}{\cos\phi}, \sin\theta=\frac{y}{\cos\phi}$ so that

$$
\theta_c = \verb|ConstrAngle|\left(\mathrm{sgn}~y_c, \frac{x_c}{\sqrt{1-z_c^2}}\right)
$$

### Edge cases

The above derivation is not applicable when $\det M=0$, which occurs when $\cos\phi_p=0$ for some $p\in\{a,b\}$ (one of the points is on a pole) or $\sin(\theta_b-\theta_a)=0$ (both points lie on the same or opposite line of longitude).
<!--This edge case is collectively defined by the property

$$
\cos\phi_a\cos\phi_b\sin(\theta_b-\theta_a)=0 ~. \tag{3}
$$
-->

Without loss of generality, assume that $\cos\phi_a=0$ or $\sin(\theta_b-\theta_a)=0$.
The geodesic distance formula gives

$$
\begin{pmatrix}
    \cos d_{ac} \\ \cos d_{bc}
\end{pmatrix}
= \underbrace{
\begin{pmatrix}
    \sin\phi_a & m\cos\phi_a \\
    \sin\phi_b & \cos\phi_b
\end{pmatrix}
}_{N}
\begin{pmatrix}
    \sin\phi_c \\ \cos\phi_c\cos(\theta_b-\theta_c)
\end{pmatrix}
$$

where $m=\cos(\theta_a-\theta_c)/\cos(\theta_b-\theta_c)$ is relevant only when $\sin(\theta_b-\theta_a)=0$, and attains the value $1$ if $\theta_a=\theta_b$, and $-1$ if $\theta_a=\pi+\theta_b$.
The matrix $N$ is guaranteed to have nonzero determinant under our assumptions and is thus invertible, so that

$$
\begin{pmatrix}
    \sin\phi_c \\ \cos\phi_c\cos(\theta_b-\theta_c)
\end{pmatrix}
=
\frac{1}{\det N}
\begin{pmatrix}
    \cos\phi_b & -m\cos\phi_a \\
    -\sin\phi_b & \sin\phi_a
\end{pmatrix}
\begin{pmatrix}
    \cos d_{ac} \\ \cos d_{bc}
\end{pmatrix} ~,
$$

from which one can read off $\sin\phi_c$ directly.
Since the cosine of a latitude is always non-negative,

$$
\phi_c = \verb|ConstrAngle|\left(\mathrm{sgn}\sin\phi_c, \sqrt{1-\sin^2\phi_c}\right) ~.
$$

The latter matrix equation also allows one to read off $\cos(\theta_b-\theta_c)$ directly, so that

$$
\theta_b-\theta_c = \begin{cases}
\verb|ConstrAngle|\left(\pm1, (0,1)\cdot N^{-1}\mathbf{d'}\right) \qquad& \text{if $\cos(\theta_b-\theta_c)\not\in\{-1,1\}$,} \\
\verb|ConstrAngle|\left(1, (0,1)\cdot N^{-1}\mathbf{d'}\right) & \text{otherwise.}
\end{cases}
$$

One can then unambiguously obtain the possible values of $\theta_c$ by subtracting the previous expression from $\theta_b$.


```python

```

5b07f35eebf8d33a6ddc9b106e3d869965b9a318
ipynb
Jupyter Notebook
constrsphtrig.ipynb
asemptote/constrsphtr
2ceffc80a748b9184b44f2659e88e516552decb5
MIT
2021-08-08T12:13:47.000Z
2021-08-08T12:13:47.000Z
constrsphtrig.ipynb
asemptote/constrsphtrig
2ceffc80a748b9184b44f2659e88e516552decb5
MIT
constrsphtrig.ipynb
asemptote/constrsphtrig
2ceffc80a748b9184b44f2659e88e516552decb5
MIT
Qwen/Qwen-72B
1. YES
2. YES