Intersection of circle and geodesic segment on sphere

I am trying to find an efficient way of computing the intersection point(s) of a circle and line segment on a spherical surface.

Say you have a sphere of radius R. On the surface of this sphere are

1. a circle with center ($\theta_c$,$\phi_c$) and radius r
2. a geodesic line segment defined by endpoints ($\theta_1$,$\phi_1$) and ($\theta_2$,$\phi_2$)

where $\theta$ is the colatitude in $[0,\pi]$, $\phi$ is the longitude in $[0,2\pi]$, and $r$ is measured by the geodesic distance on the sphere (not straight line distance in Euclidean space). How would you

1. determine whether the circle and line intersect at all, including whether the segment is contained by the circle?
2. compute the intersection point(s)?

We can assume there is nothing pathological going on. $r$ is not zero and is not so large that it's bigger than the sphere, the line's endpoints are not identical, etc.

Solution 1:

If everything is converted from geographic coordinates to Cartesian vectors, the calculation can be done by linear algebra and trigonometry.

Let $P_{0}$ denote the center of the circle $C$; $P_{1}$ and $P_{2}$ the endpoints of the geodesic segment; and assume

• $P_{1} \neq \pm P_{2}$, so there is a unique short arc of great circle from $P_{1}$ to $P_{2}$;
• $r < \pi R/2$, so $C$ is not a great circle.

The great circle through $P_{1}$ and $P_{2}$ is the intersection of the sphere with the plane $\Pi$ containing $P_{1}$, $P_{2}$, and the origin (the center of the sphere), which has unit normal vector $$ N = \frac{P_{1} \times P_{2}}{\|P_{1} \times P_{2}\|}. $$ The angle subtended at the center of the sphere by the center of $C$ and a point of $C$ is $\theta = r/R$. Further, $\Pi$ intersects $C$ if and only if $N$ makes angle no larger than $\theta$ with the "equator" perpendicular to $P_{0}$, if and only if $$ \biggl|\frac{\pi}{2} - \arccos \frac{P_{0} \cdot N}{R}\biggr| = \biggl|\arcsin \frac{P_{0} \cdot N}{R}\biggr| \leq \theta. $$ If this inequality fails to hold, the plane $\Pi$ (and hence the geodesic segment) does not intersect $C$. If equality holds, $C$ is tangent to $\Pi$, and if strict inequality holds $C$ crosses $\Pi$.

Suppose the preceding inequality is satisfied. The Euclidean center of $C$ is $c_{0} = (\cos\theta)P_{0}$.

Let $P_{2}'$ denote the vector on the short arc starting at $P_{1}$ and heading toward $P_{2}$ that is perpendicular to $P_{1}$; namely, scale the orthogonal projection $$ P_{2} - \frac{P_{1} \cdot P_{2}}{R^{2}} P_{1} $$ to have magnitude $R$. Points on the great circle have the form $$ Q = (\cos t)P_{1} + (\sin t)P_{2}' $$ for some real $t$. Such a point lies on $C$ if and only if $(Q - c_{0}) \cdot P_{0} = 0$, or $$ (\cos t)P_{1} \cdot P_{0} + (\sin t)P_{2}' \cdot P_{0} = c_{0} \cdot P_{0} = R^{2}\cos\theta. $$ If we set $$ A = P_{1} \cdot P_{0},\qquad B = P_{2}' \cdot P_{0},\qquad C = R^{2}\cos \theta, $$ this equation is written $A\cos t + B\sin t = C$. To solve, let $t_{0} = \operatorname{atan2}(B, A)$ be the unique angle in $[0, 2\pi)$ satisfying $$ \cos t_{0} = \frac{A}{\sqrt{A^{2} + B^{2}}},\qquad \sin t_{0} = \frac{B}{\sqrt{A^{2} + B^{2}}}. $$ By the addition formula for cosine, we have $$ \cos(t - t_{0}) = \frac{R^{2}\cos\theta}{\sqrt{(P_{1} \cdot P_{0})^{2} + (P_{2}' \cdot P_{0})^{2}}}, $$ which yields $$ t = t_{0} \pm \arccos\frac{R^{2}\cos\theta}{\sqrt{(P_{1} \cdot P_{0})^{2} + (P_{2}' \cdot P_{0})^{2}}}. $$ This gives two points on the great circle that lie on $C$.

It remains to check whether either point lies on the short arc of great circle from $P_{1}$ to $P_{2}$. But $Q = (\cos t)P_{1} + (\sin t)P_{2}'$ lies on the arc of great circle from $P_{1}$ to $P_{2}$ if and only if $0 \leq t \leq \arccos(P_{1} \cdot P_{2}/R^{2})$. Alternatively, a point $Q$ of the great circle lies between $P_{1}$ and $P_{2}$ (i.e., on the short arc) if and on if the angle from $P_{1}$ to $Q$ plus the angle from $Q$ to $P_{2}$ equals the angle from $P_{1}$ to $P_{2}$, i.e., $$ \arccos \frac{P_{1} \cdot Q}{R^{2}} + \arccos \frac{P_{2} \cdot Q}{R^{2}} = \arccos \frac{P_{1} \cdot P_{2}}{R^{2}}. $$ (In the diagram, these were checked numerically and colored accordingly, green if yes and red if no.)