vec_vs_angle.md

Aviation Formulary Revisited

I do not come in peace; I am here on a mission. I wish to vanquish the angles and trig which run amok
To my knowledge, this seems to stem from the thought:

The user provides and requests angles, therefore my formulas should input and output angles

At the very least, this is untrue in that the user provides degrees and we compute with radians
But more pervasively, I argue that the formulas are first and foremost Geometric and not Trigonometric

To make my case, I will go through common formulas both in terms of latitude, longitude ($\phi$, $\lambda$), and unit vectors of $ℝ^3$
Even if you remain unconvinced that you should prefer the use of the vectors, over angles and trigonometry,
I hope it provides a more intuitive form and derivation for why they work and where they come from

---

Conversions between latitude, longitude $({\phi \atop \lambda})$ and unit vectors of $ℝ^3$

We'll use $N$ to denote the North pole, defined as $\ \phi = \frac{\pi}{2} \ \ \equiv \ \hat{z}$
and $L$ where the Prime meridian meets the Equator $\ ({\phi = 0 \atop \lambda = 0}) \ \ \equiv \hat{x}$

$$ \begin{align*} \begin{pmatrix} \phi \\ \lambda \end{pmatrix} \ \overset{V}{\longrightarrow} \begin{pmatrix} \cos\phi\ \cos\lambda \\ \cos\phi\ \sin\lambda \\ \sin\phi \end{pmatrix} \\ \begin{pmatrix} x \\ y \\ z \end{pmatrix} \ \overset{𝒜}{\longrightarrow} \begin{pmatrix} \text{atan2} \left( \frac{z}{\sqrt{\smash[b]{x^2 + y^2}}} \right) \\ \text{atan2}\left(\frac{y}{x}\right) \end{pmatrix} \end{align*} $$

Important that from here on out, $\text{atan2} ( \thinspace \tfrac{y}{x} \thinspace )$ denotes atan2(y, x) as two parameters without division
And, while $V$ always outputs a unit vector, the input of $𝒜$ need not be normalized

---

Rotations

We'll start by doing a head-count of length preserving linear transforms. Namely, reflections and rotations
Reflection through the origin is simple enough: $P \to {-P}$
Reflection through plane $\hat{n}$ is similarly familiar: $P \to P - 2 \thinspace (P \cdot \hat{n})\hat{n}$

Rotations, on the other hand, deserve a bit more attention. We'll compare matrices and quaternions, which for our purposes can be thought of as a convenient and lightweight representation of rotations around arbitrary axes
While I won't be giving an in-depth exposition, I'll briefly list the notation, descriptions, and costs

$$ \begin{align*} &\dim(ℝ^{3 \times 3}) = 9 &\text{matrix size} &\text{9 floats} \\ &T = T_b T_a \ \ \in ℝ^{3 \times 3} &\text{matrix composition} &\text{18 adds, 27 muls} \\ &v\prime = Tv \ \ \in ℝ^3 &\text{matrix application} &\text{ 6 adds, 9 muls} \\ \\ \\ &\dim(ℍ) = 4 &\text{quaternion size} &\text{4 floats} \\ &q = \exp \left( \tfrac{\theta}{2}\hat{n} \right) = \cos\tfrac{\theta}{2} + \sin\tfrac{\theta}{2}\hat{n} \ \ \in ℍ &\text{quaternion construction} &\text{1 sincos, 4 muls} \\ &q = q_b q_a \ \ \in ℍ &\text{quaternion composition} &\text{12 adds, 16 muls} \\ &v\prime = q[v] = q v q^* \ \in ℝ^3 &\text{quaternion application}^\text{[1]} &\text{12 adds, 18 muls} \\ &T = \left[ \ q[\hat{x}] \ \ q[\hat{y}] \ \ q[\hat{z}] \ \right] \ \in ℝ^{3 \times 3} &\text{quaternion to matrix}^\text{[2,3]} &\text{12 adds, 12 muls} \end{align*} $$

[1] Fabian "ryg" Giesen - Rotating a single vector using a quaternion
[2] Marc B. Reynolds - On quaternion/rotation matrix conversions and errors
[3] Matrix construction is to construct a quaternion and convert to matrix, $\text{1 sincos, 12 adds, 16 muls}$

So if we need to rotate a single vector $v$ around $\hat{n}$ by $\theta$

• Construct then apply matrix: $\text{1 sincos, 18 adds, 25 muls}$
• Construct then apply quat: $\text{1 sincos, 12 adds, 22 muls}$

A couple things to note:

• To transform more than 1 vector by the same rotation, matrix counts are increasingly cheaper
• These are scalar op counts, they don't account for fma or SIMD usage, where matrices fare better
• Op counts for quaternion application and conversion presume unit-quaternions.
  All quaternions of this post are normalized by construction so I'll make no further note of this

While the rest of this post treats quaternions as a drop-in for rotations, it would be remiss of me not to mention, at least in passing, the much richer structure underlying the quaternions beyond just their algebraic and computational properties. While not the topic of this post, I may one day feel motivated to expand on this later

---

Antipode of $P$

$$ \begin{align*} &\phi\prime = {- \thinspace \phi} \\ &\lambda\prime = \begin{cases} \lambda + \pi \quad \text{if }\ \lambda \lt 0 \\ \lambda - \pi \end{cases} \end{align*} $$

$$ \begin{align*} P\prime = {-P} \end{align*} $$

The need to ensure $\lambda$ is normalized within $\pm\pi$ will be a common pattern for angles but thankfully absent for vectors

---

Great Circle Distance $\theta$ from $A$ to $B$

$$ \cos\theta = \sin\phi_A\ \sin\phi_B \ + \ \cos\phi_A\ \cos\phi_B\ \cos\ \Delta \lambda $$ $$ \cos\theta = A \cdot B $$

It's important to note that this computes the $\cos\theta$ rather than $\theta$ itself
The reason is simple, for most computations $\cos\theta$ and/or $\sin\theta$ are the values you actually end up using

There are, however, some computations which require $\theta$, though I claim that $\text{acos}$ and $\text{asin}$ are still a poor fit
So while I feel the need to stress the importance in avoiding $\text{acos}$ and $\text{asin}$, it doesn't provide a full picture

Numerically, $\text{acos}(v)$ is ill-conditioned in near-parallel cases, while $\text{asin}(v)$ is ill-conditioned in near-perpendicular cases
The role of $\text{atan2}(y,x)$ is to be well-conditioned for all cases without requiring clamping to [-1,1] to avoid NaNs

$$\theta = \text{atan2} \left( \frac{\sqrt{\smash[b]{(\cos\phi_B\ \sin\ \Delta \lambda)^2 + (\cos\phi_A\ \sin\phi_B - \sin\phi_A\ \cos\phi_B\ \cos \Delta \lambda)^2}}}{\sin\phi_A\ \sin\phi_B \ + \ \cos\phi_A\ \cos\phi_B\ \cos\ \Delta \lambda} \right) $$ $$\theta \ = \ \text{acos}\ A \cdot B \ = \ \text{atan2} \left( \frac{\lVert{A \times B}\rVert}{A \cdot B} \right) \ = \ 2 \cdot \text{atan2} \left( \frac{\lVert{A-B}\rVert}{\lVert{A+B}\rVert} \right) $$

First is the haversine formula, where the denominator is $\cos\theta$ shown above, hinting the numerator to be $\sin\theta$
The vector expressions are refreshingly clear. The last one returning twice the bisected angle

We can similarly improve the robustness of $a \times b$ by considering

$$ (a - b) \times (a+b) = 2 \ (a \times b) $$

Which allows us to always compute the cross-product of orthogonal vectors
And use Difference of Products if even more precision is required (and fma is available)

---

True Course $\gamma$ from $A$ to $B$ (when $A \neq \pm N$)

$$ \gamma = \begin{cases} h \quad \text{if } \sin(\lambda_A-\lambda_B) < 0 \\ \tau - h \end{cases} \quad h = \text{acos} \left( \frac{\sin\phi_B - \sin\phi_A\ \cos\theta}{\cos\phi_A\ \sin\theta} \right) $$ $$ \gamma = \text{atan2} \left( \frac{\cos\phi_B\ \sin(\lambda_B-\lambda_A)}{\cos\phi_A\sin\phi_B \ - \ \sin\phi_A\ \cos\phi_B\ \cos(\lambda_B-\lambda_A)} \right) $$

$$ \begin{align*} &{\atop \cos\ \gamma =} {\over{(A \times B)}} {\atop \cdot} {\over{(A \times N)}} \\ &{\atop \sin\ \gamma =} {\over{(A \times B)}} {\atop \cdot} {\atop(} {\over{(A \times N)}}{\atop \times A)} \end{align*} $$

---

Find point $B$ a given distance $\theta$ from $A$, and course $\gamma$ from North

$$ \begin{align*} &\phi_B = \text{asin}\left(\sin\phi_A\ \cos\theta \ + \ \cos\phi_A\ \sin\theta\ \cos\gamma\right) \\ &\lambda_B = \text{mod}\left(\lambda_A \ - \ \text{atan2}\left(\frac{\cos\phi_A\ \sin\theta\ \sin\gamma}{\cos\ \theta \ - \ \sin\phi_A\cos\phi_B}\right) \ + \ \pi,\ \ \tau\right) \ - \ \pi \end{align*} $$

$$ \begin{align*} &{\atop \text{let }\ q_\theta = \exp(\ \tfrac{\theta}{2} } \ {\over{A \times N}} {\atop \ ) } \\ & \text{let }\ q_\gamma = \exp({\text{-} \thinspace \tfrac{\gamma}{2}} \ A \ ) \\ &B = (q_\gamma \ q_\theta) \ [A] \end{align*} $$

The last expression is doing a quaternion product of $q_\theta$ and $q_\gamma$ which is applied to $A$
Reading right-to-left, we start at $A$, move "up" (towards north) until $\theta$, then rotating clockwise around $A$ by $\gamma$

---

Point $P$ given $\lambda$ on Great Circle of $AB$

$$ \begin{align*} \phi_0 = \text{atan2}\left( \frac{\sin\phi_A\ \cos\phi_B\ \sin(\lambda-\lambda_B) \ - \ \sin\phi_B\ \cos\phi_A\ \sin(\lambda-\lambda_A)} {\cos\phi_A\ \cos\phi_B\ \sin(\lambda_A-\lambda_B)} \right) \\ \phi_P = \begin{cases} \phi_0 - \pi \quad \text{if }\ \phi_0 \ > \ {+\tfrac{\pi}{2}} \\ \phi_0 + \pi \quad \text{if }\ \phi_0 \ < \ {-\tfrac{\pi}{2}} \\ \phi_0 \end{cases} \end{align*} $$

$$ \begin{align*} & \text{let }\ q_\lambda = \exp(\ \tfrac{\lambda}{2}\ N ) \\ &{\atop \text{let }\ Q \ = \ }{\over q_\lambda[ L \times N ] \ \times \ (A \times B)} \\ &P = \begin{cases} {+Q} \quad \text{if }\ Q \cdot q_\lambda[ L ] \ge 0 \\ {-Q} \end{cases} \end{align*} $$

We start by considering the $\lambda=0$ case first, then later apply a rotation by $\lambda$ to generalize
We form $L \times N$, a plane through the origin containing $\lambda=0$ then intersect it with $A \times B$, the plane of the great circle
The intersection of two distinct planes through the origin is a line through the origin, so the correct result is either $\pm Q$
We must now consider that our plane contains both $\lambda=0$ and $\lambda = \pm \tau / 2$, so we select $L$ which bisects $\lambda=0$
Now test whether $+Q$ is within a quarter-turn of $L$. If it's not, then $-Q$ must be within a quarter-turn of $L$

Alternatively, consider the plane with normal $\hat{n} = L$ whose upper half-space, determined by $\forall P\ (P \cdot \hat{n} \ge 0)$, contains $\lambda=0$
While distance to a point may initially be more intuitive, planes appear often enough that it's worthwhile to build familiarity

Now, that we've solidified our understanding we can apply some optimizations

$$ \begin{gather*} q_\lambda[L \times N] && q_\lambda[L] \\ q_\lambda[{-\hat{y}}] && q_\lambda[\hat{x}] \\ {-T_y} && T_x \\ \end{gather*} $$ $$ \text{where $T$ is the matrix of } q_\lambda $$

So instead of doing 2 quaternion applications, costing $\text{24 adds, 36 muls}$
We need only convert to a matrix, costing merely $\text{12 adds, 12 muls}$
And since the axis of rotation is $N = \hat{z}$ this cost becomes $\text{0 adds, 0 muls, 1 sincos}$

I have no issue with applying these optimizations. In fact, I enjoy and encourage doing so
But it should come only after the principles are understood

---

Midpoint $M$ of $A$ to $B$ on Line of constant $\lambda$

$$ \begin{align*} \phi_M = \tfrac{1}{2}(\phi_A + \phi_B) \end{align*} $$

$$ \begin{align*} {{\atop M = } {\over{A + B}}} \end{align*} $$

---

Midpoint $M$ of $A$ to $B$ on Line of constant $\phi$

$$ \begin{align*} &\text{let }\ \Delta\lambda = \lvert \lambda_B - \lambda_A \rvert \\ &\text{let }\ \lambda_0 = \text{min}(\lambda_A,\ \lambda_B) \ + \ \tfrac{\Delta\lambda}{2} \\ &\lambda_M = \begin{cases} \lambda_0 \quad \text{if }\ \Delta\lambda \leq \pi \\ \text{mod}(\lambda_0 + \tau,\ \tau) - \pi \end{cases} \end{align*} $$

$$ \begin{align*} &\text{let }\ h = A \cdot N \ \text{ and }\ \ell = \sqrt{1 - h²} \\ &{{\atop M =} {\over{(A-hN) \ + \ (B-hN)}}}{\atop \ell + hN} \end{align*} $$

---

Midpoint $M$ of $A$ to $B$ on Great Circle

$$ \begin{align*} \text{let }\ B_x = \cos\phi_B\ \cos(\lambda_B - \lambda_A) \\ \text{let }\ B_y = \cos\phi_B\ \sin(\lambda_B - \lambda_A) \end{align*} $$ $$ \begin{align*} &\phi_M = \text{atan2}\left(\frac{\sin\phi_A\ +\ \sin\phi_B}{\sqrt{\smash[b]{(\cos\phi_A + B_x)²\ +\ {B_y}²}}}\right) \\ \\ &\lambda_M = \lambda_A + \text{atan2}\left(\frac{B_y}{\cos\phi_A + B_x}\right) \end{align*} $$

$$ \begin{align*} {{\atop M = } {\over{A + B}}} \end{align*} $$

Lines of constant $\lambda$ are merely Great Circles which pass through $\pm N$

---

Midpoint $M$ of $A$ to $B$ on Small Circle with pole $\hat{n}$

$$ \text{Finding } \binom{\phi_M}{\lambda_M} \text{ is left as an exercise for the reader...} $$

$$ \begin{align*} &\text{let }\ h = A \cdot \hat{n} \ \text{ and }\ \ell = \sqrt{1 - h²} \\ &{{\atop M =} {\over{(A-h\hat{n}) \ + \ (B-h\hat{n})}}}{\atop \ell + h\hat{n}} \end{align*} $$

Lines of constant $\phi$ are merely Small Circles with pole $\hat{n} = N$

---

Interpolate $A$ to $B$ by $t \in [0, 1]$

$$ 𝒜(\text{slerp}(V(A),\ V(B),\ t))$$ $$ \text{slerp}( A,\ B,\ t)$$

While the others had some plausible deniability, this just converts angles to vectors, calls slerp, and converts back
Ironically, the vector $\text{slerp}$ is doing something similar. Computing an angle, and converting the angle back to a vector

At a high level, we first compute $\theta_t = \text{lerp}(0, \theta, t) = \theta \cdot t$ where $\theta$ is the angle between $A$ and $B$
Then take $\{A,\ A^{\perp}\}$ as a 2D Graham-Schimdt of $\{A,\ B\}$ and return $(\cos\theta_t \cdot A \ + \ \sin\theta_t \cdot A^{\perp})$

We may be tempted to use $\text{normalize(lerp}(A,\ B,\ t))$ but this fails the requirement $\forall t \left(A \cdot \text{interpolate}(A,\ B,\ t) = \cos\theta_t \right)$
It will still be geodesic, monotonic, and symmetric, with smaller arcs yielding smaller errors; so worthwhile to consider
Furthermore, t can be refit to approximate $\text{slerp}$. (this approximates quaternion slerp, with slightly different considerations)

Exercise: Subdivide $A$, $B$ into arcs no longer than $\theta$
Find $n = \lceil \theta_{AB} / \theta \rceil$ then compute $P_i = \text{slerp}(A,\ B,\ i/n)$ for $i \in [0, n]$
This is a good enough starting point, but it can be optimized further :)

---

Distance $\theta$ from $P$ to Great Circle segment $AB$

Starting from this StackOverflow post (with modifications and corrections of my own)

$$ \begin{align*} &\text{let }\ \Delta\gamma = \begin{cases} 0 + \left|\ \gamma_{AP} - \gamma_{AB} \right| \ &\text{if }\ \left|\ \gamma_{AP} - \gamma_{AB} \right| < \pi \\ \tau - \left|\ \gamma_{AP} - \gamma_{AB} \right| \end{cases} \\ &\text{let }\ \cos\theta_{PQ} = \sqrt{\smash[b]{1 - \sin^2\theta_{AP}\ \sin^2(\gamma_{AP} - \gamma_{AB})}} \\ \\ &\text{let }\ \cos\theta_{AQ} = \frac{\cos\theta_{AP}}{\cos\theta_{PQ}} \\ &\cos\theta = \begin{cases} \text{max}(\cos\theta_{AP},\ \cos\theta_{BP}) \quad &\text{if }\ \Delta\gamma > \tfrac{\pi}{2} \\ \cos\theta_{BP} \quad &\text{if }\ \cos\theta_{AQ} < \cos\theta_{AB} \\ \cos\theta_{PQ} \end{cases} \end{align*} $$

$$ \begin{align*} &{\atop \text{let }\ \hat{n} = }{\over{A \times B}} {\atop \ \text{ and }\ h = P \cdot \hat{n}}\\ &\text{let }\ Q = P - h\ \hat{n} \\ &\cos\theta = \begin{cases} P \cdot \overline{Q} \quad \text{if both } {{Q \cdot (\hat{n} \times A) > 0} \atop {Q \cdot (B \times \hat{n}) > 0}} \\ \text{max}(P \cdot A,\ P \cdot B) \end{cases} \end{align*} $$

When we try to implement this in code:

// cosTheta from P to arc AB using angles
let bearAB = Bearing(latA, lonA, latB, lonB);          // 1*atan2 + 3*sincos
let bearAP = Bearing(latA, lonA, latP, lonP);          // 1*atan2 + 3*sincos
let [cosAP, sinAP] = CS_Dist(latA, lonA, latP, lonP);  // 1*sqrt  + 3*sincos
let  cosBP         =  C_Dist(latB, lonB, latP, lonP);  //           3*sincos

var diff = abs(bearAP - bearAB);
if(diff > .pi){ diff = .tau - diff; }

if(diff > .pi/2){
  return max(cosAP, cosBP);                  // min_total: 2*atan2 + 1*sqrt + 12*sincos
}

let cosAB = C_Dist(latA, lonA, latB, lonB);  // 3*sincos
let sinPQ = sinAP*sin(bearAP - bearAB);      // 1*sin/cos
let cosPQ = sqrt(1 - sinPQ*sinPQ);           // 1*sqrt
let cosAQ = cosAP / cosPQ;
return (cosAQ < cosAB ? cosBP : cosPQ);      // max_total: 2*atan2 + 2*sqrt + 16*sincos

// cosTheta from P to arc AB using vectors
let n = normalize(cross(A, B));              // 1*sqrt
let h = dot(P, n);
let Q = P - h*n;
let inA = dot(Q, cross(n, A)) > 0;
let inB = dot(Q, cross(B, n)) > 0;
return (inA && inB ? dot(P, normalize(Q))    // 1*sqrt
                   : max(dot(P, A), dot(P, B)));

While the pure vector cost is 2 sqrt's, if we consider the cost of $V:\ (\phi, \lambda) \to ℝ^3$ there's an additional 6 sincos
Still much more efficient than the proposed method using angles, but let's try and do better

While projecting $P$ to $Q$ is helpful for our initial formulation, it turns out that this is unecessary
In order to avoid computing $Q$ and $\overline{Q}$, we'll start by observing:

$$ \begin{align*} Q \cdot (\hat{n} \times A) &= (P - h\ \hat{n}) \cdot (\hat{n} \times A) \\ &= P \cdot (\hat{n} \times A) \ - \ (h\ \hat{n} \cdot (\hat{n} \times A)) \\ &= P \cdot (\hat{n} \times A) \qquad \text{since }\ \hat{n} \perp \hat{n} \times A \\ \\ Q \cdot (B \times \hat{n}) &= P \cdot (B \times \hat{n}) \qquad \text{likewise} \end{align*} $$

Since the projection is orthogonal, it forms a right triangle $\triangle OQP$ allowing us to substitute:

$$ \cos\theta = P \cdot \overline{Q} = \sqrt{1-h²} $$

The use of $h^2$ instead of $h$ allows us to go even further:

$$ \begin{gather*} &h &= &P \cdot \hat{n} &= &P \cdot \tfrac{\overset{⇀}n}{\lVert \overset{⇀}n \rVert} & \\ \Rightarrow &h^2 &= &\frac{(P \cdot \overset{⇀}n)^2}{\lVert \overset{⇀}n \rVert^2} &= &\frac{d \cdot d}{\overset{⇀}n \cdot \overset{⇀}n} & \text{where }\ d = P \cdot \overset{⇀}n \end{gather*} $$

$$ \begin{align*} \cos\theta \ = \ P \cdot \overline{Q} \ = \ \sqrt{1 \ - \ \tfrac{d \cdot d}{\overset{⇀}n \cdot \overset{⇀}n}\ } \end{align*} $$

Which removes the need for normalizing $\overset{⇀}n$ bringing the cost to a single sqrt

$$ \begin{align*} &\text{let }\ {\overset{⇀} n} = A \times B \ \text{ and }\ d = P \cdot {\overset{⇀} n} \\ &\cos\theta = \begin{cases} \sqrt{\smash[b]{1 \ - \ \tfrac{d \cdot d}{{\overset{⇀} n} \cdot {\overset{⇀} n}}}} \quad \text{if both } {{P \cdot ({\overset{⇀} n} \times A) > 0} \atop {P \cdot (B \times {\overset{⇀} n}) > 0}} \\ \text{max}(P \cdot A,\ P \cdot B) \end{cases} \end{align*} $$

let n = cross(A, B);
let d = dot(P, n);
let inA = dot(P, cross(n, A)) > 0;  // 5 adds, 9 muls
let inB = dot(P, cross(B, n)) > 0;  // 5 adds, 9 muls
return (inA && inB ? sqrt(1 - d*d / dot(n, n))   // 1*sqrt
                   : max(dot(P, A), dot(P, B)));

While we could call it quits here, lets see if we can take one last stab at it
If we inline ${\overset{⇀} n} = A \times B$ we can apply the vector triple product and hoist common subexpressions

$$ \begin{align} {\overset{⇀} n} \times A &= \ \ (A \times B) \times A &= \ \ B(A \cdot A) - A(A \cdot B) &= \ \ B - A(A \cdot B) \\ B \times {\overset{⇀} n} &= \ \ B \times (A \times B) &= \ \ A(B \cdot B) - B(B \cdot A) &= \ \ A - B(A \cdot B) \\ \\ &\quad &\quad \text{since }\ A \cdot A = B \cdot B = 1 & \end{align} $$

$$ \begin{align*} P \cdot ({\overset{⇀} n} \times A) > 0 &\equiv \quad P \cdot (B - A(A \cdot B)) > 0 &\equiv \quad P \cdot B - (A \cdot B) (P \cdot A) > 0 &\equiv \quad P \cdot B > (A \cdot B) (P \cdot A) \\ P \cdot (B \times {\overset{⇀} n}) > 0 &\equiv \quad P \cdot (A - B(A \cdot B)) > 0 &\equiv \quad P \cdot A - (A \cdot B) (P \cdot B) > 0 &\equiv \quad P \cdot A > (A \cdot B) (P \cdot B) \\ \end{align*} $$

let PA = dot(P, A);
let PB = dot(P, B);
let AB = dot(A, B);   // 2 adds, 3 muls
let inA = PB > AB*PA; // 1 add,  1 mul
let inB = PA > AB*PB; // 1 add,  1 mul
if inA && inB {
  let n = cross(A, B);
  let d = dot(P, n);
  return sqrt(1 - d*d / dot(n, n));
}
return max(PA, PB);

Takes us from $\text{10 adds, 18 muls}$ down to $\text{4 adds, 5 muls}$ since we already take dot(P, A) and dot(P, B)
Breaking it down by case: interior saves $\text{2 adds, 7 muls}$ while exterior saves $\text{8 adds, 18 muls}$

---

Distance $\theta$ from $P$ to major Great Circle segment $AB$

Now that we know how to do minor arcs, we could extend this to major arcs as such:

Code	Angle Cost	Vector Cost
dist_major(P, A, B){ let M = antipode(midpoint(A, B)); let da = dist_minor(P, A, M); let db = dist_minor(P, B, M); return min_dist(da, db); }	2*atan2 + 1*sqrt + 3*sincos 2*atan2 + 2*sqrt + 16*sincos 2*atan2 + 2*sqrt + 16*sincos = 6*atan2 + 5*sqrt + 35*sincos	1*sqrt 1*sqrt 1*sqrt = 3*sqrt

Which seems innocent enough, but this is where the individual costs start adding up
While the vector code is cheaper, we're also in a better position to reason about it
Since we're detecting interior vs exterior, we can simply flip those cases and be done

Minor Arc	Major Arc
let PA = dot(P, A); let PB = dot(P, B); let AB = dot(A, B); let inA = PB > AB*PA; let inB = PA > AB*PB; if inA && inB { let n = cross(A, B); let d = dot(P, n); return sqrt(1 - d*d / dot(n, n)); } return max(PA, PB);	let PA = dot(P, A); let PB = dot(P, B); let AB = dot(A, B); let inA = PB >= AB*PA; let inB = PA >= AB*PB; if inA && inB { return max(PA, PB); } let n = cross(A, B); let d = dot(P, n); return sqrt(1 - d*d / dot(n, n));

Bringing us back down to 1*sqrt

---

Distance $\theta$ from $P$ to Small Circle segment $AB$ with pole $\hat{n}$

For Great Circles, $\hat{n} = A \times B$ is taken to mean “the minor arc of the Great Circle from $A$ to $B$ around $\hat{n}$”
It's a baked in assumption that Great Circles are used as geodesic/shortest paths, so of course we want the minor arc, and if someone wants the major arc, then they should be explicit about that

So when we're instead given a pole $\hat{n}$ equidistant from $A$ and $B$, it raises an interesting question:
Do we go from A to B “along the minor arc” or “around $\hat{n}$”? (since they're no longer one and the same)
I argue for “from $A$ to $B$ around $\hat{n}$” allows for more flexibility since if $\hat{n}$ is a minor arc, then ${-\hat{n}}$ is a major arc

We can check whether $\hat{n}$ and $A \times B$ point in the same direction, in which case we use the minor arc within the intersection of two half-spaces, but if they point in opposite directions, we use the major arc within the union of those half-spaces

$$ \begin{align*} &s = {-\text{sgn}(\hat{n} \cdot (A \times B))}\\ &s_A = s\ P \cdot (\hat{n} \times A) \\ &s_B = s\ P \cdot (B \times \hat{n}) \\ &\text{inside}\ \equiv\ s \cdot \text{max}(s_A,\ s_B) > 0 \end{align*} $$

The inspiration here comes from SDF operations and min/max identities (also a kind of De Morgan's law)
Since SDFs use the convention $d&lt;0$ for inside, $s = {-\text{sgn(...)}}$ flips the signs which were previously tested for $d&gt;0$
This way, $\text{max($s_A$, $s_B$)}$ is read “the intersection of $s_A$ and $s_B$” to match the SDF convention for $\text{max}$
When the sign of $s$ flips, it turns $\text{max}(s_A,\ s_B) \to {-\text{max}}({-s_A},\ {-s_B}) \equiv \text{min}(s_A,\ s_B)$ which is now “the union of $s_A$ and $s_B$”

So with that out of the way, let's take a shot at computing some distances

$$ \begin{align*} &\text{let }\ \delta = \hat{n} \cdot A = \hat{n} \cdot B \ \text{ and }\ O = \delta \hat{n} \\ &\text{let }\ \ell = \sqrt{1 - \delta²} = \Vert A - O \Vert = \Vert B - O \Vert \\ &\text{let }\ Q = P - \hat{n}(\hat{n} \cdot P - \delta) \\ &{\atop \text{let }\ \overline{Q} = }{\over{Q-O}} {\atop \ \ell + O} \\ &\cos\theta = \begin{cases} P \cdot \overline{Q} \quad \text{if inside} \\ \text{max}(P \cdot A,\ P \cdot B) \end{cases} \end{align*} $$

As you may have guessed, our use of $\overline{Q}$ is only as a stepping stone towards a better solution
We'll start by working backwards from $\cos\theta = P \cdot \overline{Q}$ by solving in terms of the squared chord of the triangle

$$ \begin{align*} &\delta = \hat{n} \cdot A = \hat{n} \cdot B \ \text{ and }\ O = \delta \hat{n} \\ &\ell = \sqrt{1 - \delta²} = \Vert \overline{Q} - O \Vert = (b+c) \\ &m = \hat{n} \cdot P \\ &h = m - \delta \\ &h² = (m - \delta)² = m² + \delta² - 2 m \delta \\ \Rightarrow &h² - m² - \delta² = {-2 m \delta} \end{align*} $$

$$ \begin{align*} d² = \Vert P - O \Vert² &= \Vert P \Vert² + \Vert O \Vert² - 2\ P \cdot O \\ &= 1² + \delta² - 2 \delta m \\ &= 1 + \delta² + (h² - m² - \delta²) \\ &= 1 + h² - m² \end{align*} $$

$$ \begin{align*} &c² \ = \ d² - h² \ = \ 1 - m² \\ &b² \ = \ (\ell - c)² \ = \ \ell² + c² - 2 \ell c \end{align*} $$

$$ \begin{align*} \quad \alpha² \ = \ h² + b² &= h² + \ell² + c² - 2 \ell c \\ &= h² + (1-\delta²) + (1-m²) - 2 \ell c \\ &= 2 - 2 \delta m - 2 \ell c \\ \\ \Rightarrow \cos\theta = 1- \tfrac{\alpha²}{2} &= \ell c + \delta m \\ &= \sqrt{1-\delta²} \sqrt{1-m²} + \delta m \\ &= \sqrt{\smash[b]{(}1-\delta²\smash[b]{)(}1-m²\smash[b]{)}} + \delta m \\ \end{align*} $$

So boiling this all down, we get to a simpler, and notably familiar:

$$ \begin{align*} &\text{let }\ \delta = \hat{n} \cdot A = \hat{n} \cdot B \ \text{ and }\ m = \hat{n} \cdot P \\ &\cos\theta = \begin{cases} \sqrt{\smash[b]{(1-\delta²)(1-m²)}} + \delta m \quad \text{if inside} \\ \text{max}(P \cdot A,\ P \cdot B) \end{cases} \end{align*} $$

This familiarity serves as an intruiging sanity check
Small Circles are formed by any plane which cuts the sphere, including those through the origin forming Great Circles
When planes contain the origin, $\hat{n} \cdot A = \hat{n} \cdot B = \delta = 0 \Rightarrow \hat{n} \perp \text{span}\{A,\ B\} \text{ and } h = m = P \cdot \hat{n} - 0$
This gives an end result of $\cos\theta = \sqrt{1-h²} + 0$ which coincides with our previous formula