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).