Replaced the calculation of the complete elliptic integral of the third
kind with an implementation of the Bulirsch (1965) method, suggested by
Eric Agol at coffee last week.  This method is much simpler, but also
faster and seems to be more numerically stable.  The old R_J based
implementation has been left in place for now to aid comparisons but may
eventually be removed.

This change alleviates many of the numerical problems with using R_J and
has obviated the need for the tricks using identities, so those are also
gone.  Also tweaked the limits on the parameters y and p to (hopefully)
do a better job of preventing underflows / overflows.  While these don't
affect the results they could affect the performance.

After this change the calculations for lambda seem to be good to a few
ppm in single precision, but this has revealed some issues with the area
equations which limit the ultimate precision (especially in eta when the
discs overlap) to quite a bit worse than this.  Therefore I still
recommend using double precision for now until this can be fixed.

I'm still not completely happy with it, this doesn't yet solve the
"infinity times zero needs to equal pi/2" problem.  Some experiments
were done on this, but I have yet to come up with a good implementation,
there were problems with reducing the headroom by more than I wanted.
Work continues but for now this seems to produce more accurate results
than the original so it's getting committed anyway.