Understanding Great Circle Arcs Intersection Algorithm

The following material is the result of my attempt to understand the nice example from Jason Davies. I was puzzled about the origin of the algorithm used to find the intersection of two great circle arcs. Google helped and I discovered Roger Stafford’s post in Matlab newsgroup and the relevant Python’s implementation in the Spherical Geometry Toolkit. The algorithm You have two great circle arcs on a sphere, $a$ from point $\mathbf{a_0}$ to $\mathbf{a_1}$, and $b$ from $\mathbf{b_0}$ to $\mathbf{b_1}$, whose coordinates are expressed as longitude $\theta$ (positive going East of Greenwich) and latitude $\phi$ (positive going North).