Signed spherical angle between two great circle arcs

Question

I am trying to calculate the signed spherical angle between two intersecting great circle arcs. In my specific case, I am considering a unit sphere ($r=1$). The first great circle arc goes from $A(\lambda_A, \phi_A)$ to $B(\lambda_B, \phi_B)$, while the second great circle arc goes from $B$ (the intersection point) to $C(\lambda_C, \phi_C)$. Here, $\lambda$ is longitude, or the angle with the reference meridian, while $\phi$ is latitude, or the angle with respect to the equatorial plane, like is the case in the geographic coordinate system.

My spherical coordinates are related to cartesian coordinates as follows: $$ x = \cos(\lambda) \cos(\phi) $$ $$ y = \sin(\lambda) \cos(\phi) $$ $$ z = \sin(\phi) $$

The reference meridian is thus described by $y=0, x>0$ and the equatorial plane is defined as $z=0$.

I can calculate the angle between $AB$ and $BC$. First, I describe $A$, $B$, $C$ as cartesian vectors, e.g. $$\vec{A}= \big(\cos(\lambda_A) \cos(\phi_A),\ \sin(\lambda_A) \cos(\phi_A),\ \sin(\phi_A) \big).$$ Then I take the cross products $\vec{V_{ab}} = \vec{A}\times\vec{B}$ and $\vec{V_{bc}} = \vec{B}\times\vec{C}$. Then I take the angle between these two vectors (which describe the planes that the great circle arcs $AB$ and $BC$ respectively lie on):

$$\text{spherical angle}(AB, BC) = \frac{\vec{V_{ab}}\cdot \vec{V_{bc}}}{ ||\vec{V_{ab}}||\cdot ||\vec{V_{bc}}||}. $$

However, I want the signed spherical angle, so that the angle that $AB$ forms with $BC$ is defined clockwise between $[0, 2\pi]$.

Note: Perhaps it helps to say that $A$ in my specific case is the North Pole, i.e. $A = (0, 0, 1)$ or equivalently $\lambda_A = \frac{1}{2}\pi$. A solution that works for this particular case is fine as well.

The correct approach should thus for $A(0,\frac{1}{2}\pi)$, $B(0,0)$ and $C(\frac{1}{2}\pi, 0)$ output $\frac{1}{2}\pi$, while for $C(-\frac{1}{2}\pi, 0)$ it should output $\frac{3}{2}\pi$ (or equivalently $-\frac{1}{2}\pi$)

Answer 1

I figured out that I can check on which side $\vec{C}$ lies with respect to the plane formed by $\vec{V_{ab}}$ using the dot product $\vec{V_{ab}} \cdot \vec{C}$ and that I can make the output negative or positive depending on that.

Answer 2

For the purpose of implementing this in a computer program, it might be useful to apply some transformations and obtain a solution with a small number of operations and good numerical stability, which can be applied to a large collection of inputs.

I originally saw this formula in the very useful blog post by Brenton R S Recht: https://brsr.github.io/2021/05/01/vector-spherical-geometry.html. I borrow their notation for the scalar triple product, $\lvert \vec x, \vec y, \vec z \rvert$.

Denote the angle of interest as $\theta_{BAC}$. We can use some identities about cross products (1, 2), and the fact that our vectors are unit vectors (thus $\vec u \cdot \vec u = \lVert u \rVert = 1$), to obtain expressions for $\cos\left(\theta_{BAC}\right)$ and $\sin\left(\theta_{BAC}\right)$:

$$ \cos\left(\theta_{BAC}\right) = \frac { \left( \vec A \cdot \vec C \right) - \left( \vec B \cdot \vec A \right) \left( \vec B \cdot \vec C \right) } { \lVert \vec V_{BA} \rVert \lVert \vec V_{BC} \rVert } $$

$$ \sin\left(\theta_{BAC}\right) = \frac { \lvert \vec B, \vec A, \vec C \rvert } { \lVert \vec V_{BA} \rVert \lVert \vec V_{BC} \rVert } $$

Then we use $\tan\left(\theta_{BAC}\right) = \frac{\sin\left(\theta_{BAC}\right)}{\cos\left(\theta_{BAC}\right)}$ to obtain:

$$ \tan\left(\theta_{BAC}\right) = \frac { \lvert \vec B, \vec A, \vec C \rvert } { \left( \vec A \cdot \vec C \right) - \left( \vec B \cdot \vec A \right) \left( \vec B \cdot \vec C \right) } $$

and therefore

$$ \theta_{BAC} = \arctan \left( \frac { \lvert \vec B, \vec A, \vec C \rvert } { \left( \vec A \cdot \vec C \right) - \left( \vec B \cdot \vec A \right) \left( \vec B \cdot \vec C \right) } \right) $$

The above can be evaluated efficiently over a dataset consisting of triples of $\vec A$, $\vec B$, and $\vec C$ vectors in e.g. Python / Numpy:

import numpy as np


def scalar_triple_product_3d(v, x, y):
    assert 2 == v.ndim == x.ndim == y.ndim
    assert 3 == v.shape[1] == x.shape[1] == y.shape[1]

    c0 = x[:, 1] * y[:, 2] - x[:, 2] * y[:, 1]
    c1 = x[:, 2] * y[:, 0] - x[:, 0] * y[:, 2]
    c2 = x[:, 0] * y[:, 1] - x[:, 1] * y[:, 0]
    return v[:, 0] * c0 + v[:, 1] * c1 + v[:, 2] * c2


def dot_rows(x, y):
    assert 2 == v.ndim == x.ndim == y.ndim
    assert 3 == v.shape[1] == x.shape[1] == y.shape[1]

    return np.sum(x, y, axis=1)


def signed_spherical_angle_3d(v, x, y):
    assert 2 == v.ndim == x.ndim == y.ndim
    assert 3 == v.shape[1] == x.shape[1] == y.shape[1]

    numerator = scalar_triple_product(v, x, y)
    denominator = (
        dot_product_rows(x, y) -
        dot_product_rows(v, x) * dot_product_rows(v, y)
    )
    return np.arctan2(numerator, denominator)

This implementation of the scalar triple product was taken from https://stackoverflow.com/a/42158550/2954547. This implementation of the dot products of pairs of rows was taken from https://stackoverflow.com/a/31021818/2954547.