This is the answer to @Dan's question about using the auxiliary sphere
to solve intersection problems.

No, the auxiliary sphere doesn't let you solve for intersections
directly.  The problem is that the mapping from the ellipsoid to the
sphere depends on the geodesic (e.g., its azimuth at the equator).  Thus
the auxiliary sphere is good for solving for a single geodesic but not
for solving any problem involving more that one geodesic.

My recommended way to do geodesic intersections and interceptions
(shortest path between a point and a geodesic) is to use the gnomonic
projection.  Geodesics on a sphere map to straight lines in the gnomonic
projection and so, provided that your problem is confined to a
hemisphere, the gnomonic projection converts these problems to 2d ones.

There is no projection which preserves the straightness of geodesics for
an ellipsoid (because its curvature is not constant).  However,
GeographicLib provides a generalization of the gnomonic projection in
which geodesics are very nearly straight.  This leads to rapidly
converging algorithms for geodesic intersection and interceptions (again
provided that the points are all well within a hemisphere).  See my
answers (with code!) to gpesquero at

  https://sourceforge.net/projects/geographiclib/forums/forum/1026621/topic/4085561

Finally, I'd like to point out that I've recently converted
GeographicLib's geodesic routines to Javascript, so that you can play
around with them in Google Maps.  See

http://geographiclib.sourceforge.net/scripts/geod-google.html

http://geographiclib.sourceforge.net/scripts/geod-calc.html

(I didn't convert the gnomonic projection to Javascript yet.  That would
be reasonably straight forward.  I would also convert the azimuthal
equidistant projection at the same time, since that's a convenient way
to solve another class of geodesic problems involving "median lines".)

**ADDENDUM (2014-08-19)**

It's also possible to solve for the time of closest approach for two vessels traveling at constant speed along geodesics.  Because we know the differential properties of geodesics, it's possible to use Newton's method to get an accurate solution in a few iterations.  Code to implement this is posted at

https://sourceforge.net/p/geographiclib/discussion/1026620/thread/33ce09e0