How to get the cardinal direction from one location to another? [closed]

Consider the initial point to be $(\theta_1,\phi_1)$ and the final point to be $(\theta_2,\phi_2)$. The angle $\theta$ is chosen between $-\pi/2$ (South) and $+\pi/2$ (North), with $0$ being the equator, because that's how you wanted it. The coordinates on the unit sphere are $$ V_{i} =(\cos\theta_i\cos\phi_i,\cos\theta_i\sin\phi_i,\sin\theta_i), \qquad i=1,2 $$ where the angles should have the right indices $1$ or $2$ for the two points. Without a loss of generality, we may assume $\phi_1=0$ - we rotate the Earth around the pole-pole axis - and $\phi_2$ will represent what we originally called $\phi_2-\phi_1$; at the end, we will substitute $\phi_2-\phi_1$ for $\phi_2$ everywhere.

The shortest path between the points is along the maximal circle. It's useful to quantify the direction of the plane in which this maximal circle lies. Its normal direction the cross product of $V_1$ and $V_2$, i.e. $N = V_1\times V_2$. Note that this cross product is degenerate (zero) if the two points are antipodal i.e. if $\phi_2=\pi$ and $\theta_2=-\theta_1$. $$ N = (-\cos\theta_2 \sin\phi_2 \sin\theta_1, \cos\phi_2 \cos\theta_2 \sin\theta_1 - \cos\theta_1 \sin\theta_2, \cos\theta_1 \cos\theta_2 \sin\phi_2 ) $$ The coordinates belonging to the maximum circle are given by $N\cdot (x,y,z) = 0$ where $(x,y,z)$ is a point on the unit sphere. Needless to say, $(x,y,z)=V_i$ satisfy this condition both for $i=1$ and $i=2$.

This was really unnecessary but it could be useful for other tasks. Now, let's calculate the tangential vector $T$ along the unit sphere from the first point $V_1$ that points to the second point $V_2$. It has to be in the $V_1$-$V_2$ plane, i.e. orthogonal to $N$, but also orthogonal to $V_1$. The right solution is $T=N\times V_1$. The vector is

$$ T = (\sin\theta_1 (\cos\phi_2\cos\theta_2 \sin\theta_1 - \cos\theta_1 \sin\theta_2), \cos\theta_2 \sin\phi_2, \cos\theta_1 (-\cos\phi_2 \cos\theta_2 \sin\theta_1 + \cos\theta_1 \sin\theta_2)) $$

Now, the cardinal direction $\alpha$ is calculable from this vector easily if we rotate the vector above in the $xz$-plane by $\theta_1$ so that the original point would get mapped to the $(1,0,0)$ point at the equator. So we want to evaluate the vector $$ U = (t_1\cos \theta_1 + t_3\sin\theta_1, t_2, t_3\cos\theta_1-t_1\sin\theta_1 ) $$ and the value turns out to be simply $$ U = (0, \cos\theta_2\sin\phi_2, -\cos\phi_2 \cos\theta_2 \sin\theta_1 + \cos\theta_1 \sin\theta_2 ) $$ As expected, the first coordinate vanishes because the vector has to be orthogonal to $(1,0,0)$. The remaining two components determine the cardinal direction $\alpha$ via $$ \alpha = \arctan \left( \frac{\cos\theta_2 \sin\phi_2}{-\cos\phi_2 \cos\theta_2 \sin\theta_1 + \cos\theta_1 \sin\theta_2} \right) $$ This arctan is obviously defined mod $\pi$ only. One has to choose a value between $\pi/2$ and $3\pi/2$ if $\theta_2 < -\theta_1$ i.e. if we go more to the South than to the North. Don't forget that $\phi_2$ above is a shortcut for $\phi_2-\phi_1$.

As a check, note that for $\phi_2=0$, we obtain $\alpha=0$ (North) or $\alpha=\pi$ (South). For $\theta_1=\theta_2=0$ (both on equator), we obtain $\alpha =\arctan(1/0) = \pm \pi/2$ (East, West). As is usual with $\arctan$, the $\pm \pi$ ambiguity is unnecessary. The more accurate condition is $$ e^{i\alpha} = r_+ [(-\cos\phi_2 \cos\theta_2 \sin\theta_1 + \cos\theta_1 \sin\theta_2) + i (\cos\theta_2 \sin\phi_2)], \qquad r_+\in R^+ $$