# How accurate is GeoDistance?

Two years ago I tried an evaluation copy of Mathematica. I reported various inaccuracies in the GeoDistance routine to which I never received any fix. Could someone check whether GeoDistance is any more accurate now? Thanks.

Here are the problems I reported on June 20, 2011:

                                     Mathematica   Correct     Error
GeoDistance[{30, 0}, {-30, 180}]  -> 19928486.7   20003931.5  -75444
GeoDistance[{0.4, 0}, {0.4, 180}] -> 19915473.0   19915472.0       1
GeoDistance[{0, 0}, {0.5, 179.5}] -> 19936426.5   19936288.6     138

-
For completeness: which reference frame are you using for the "correct" versions? By default, Mathematica uses "ITRF00" as the reference frame, and GeoDistance[] uses Vincenty's method (Method -> "Vincenty75") for distances. –  Ｊ. Ｍ. Jun 6 '13 at 3:56
The "correct" results were calculated for the GRS80 ellipsoid (a = 6378137, f = 1/298.25722210088271), i.e., the reference ellipsoid used by ITRF00. Using WGS84 gives the same results (when rounded to the number of digits given here). –  cffk Jun 6 '13 at 11:30
The "Vincenty75" method has known limitations; so if GeoDistance is still using this method, the results will be suspect for nearly antipodal points. Incidentally, Vincenty himself recognized the problems with his method and supplied a fix that went some way to curing the problems. –  cffk Jun 6 '13 at 12:54
Thanks, this is good to know. However, it does make you wish that Mathematica would get rid of its three slightly flaky methods for GeoDistance and instead use one that just worked! –  cffk Jun 6 '13 at 15:48
I submitted it Wolfram support on 2011-06-27 and got a reply (TS 22876) saying it had been forwarded to "the appropriate people in our development group". I did a follow-up a few days later with a similar reply. That was the last I heard about the matter (and I couldn't pursue it further since my trial license ran out). –  cffk Jun 7 '13 at 12:59

From the information provided by @0x4A4D, the answer to my question is that GeoDistance has the same problems as it did 2 years ago. I implemented (in C++) some improved algorithms for geodesics in 2009 which are described in this paper. The benefits relative to GeoDistance are:

• accurate to round-off for $|f| < 1/50$;
• always converges;
• computes ellipsoidal areas (and some other properties).

Subsequently I've ported the algorithms to C, Fortran, Python, MATLAB, and JavaScript (see http://geographiclib.sf.net). The port to Java is in progress (it'll be ready in less than a month). I'm not in a position to write a Mathematica version because I no longer have access to Mathematica. However, if someone else wants to take the lead on this, I would be happy to assist. Alternatively, perhaps it's easier to interface Mathematica to the C or Java version (Java would be better since it presents a class interface)?

-
It hopefully shouldn't be too hard. OTOH, since Mathematica has explicit support for elliptic integrals, one might be able to use those instead of the truncated Fourier approximations in your paper. –  Ｊ. Ｍ. Jun 9 '13 at 13:53
The formulation in terms of elliptic integrals is given here and implemented in the C++ library (no other languages at this point). For WGS84 with C++, the elliptic integral method is roughly 2 times less accurate (more accumulated round-off error) and 2 times slower. On the other hand, elliptic integrals are good for any flattening (I've checked $1/100 < b/a < 100$). –  cffk Jun 9 '13 at 14:20
"roughly 2 times less accurate" - depends on the implementation, I'd say. I note that you're using the imaginary modulus formulation (the Earth's oblate); I'm sure you're aware that those expressions can be reformulated so that the moduli are real and within $[0,1)$, which should make for more stable evaluation. OTOH, the formulae also involve $D(\phi,k)$ and the elliptic integral of the third kind, and those are indeed very difficult to deal with numerically, so I'd guess that's why you're not getting full accuracy from the elliptic integral formulation. –  Ｊ. Ｍ. Jun 9 '13 at 14:27
If push comes to shove, it might be more expedient to reformulate in terms of Carlson's symmetric integrals; at least, in that formulation, computing the two special functions I mentioned earlier is less bothersome. –  Ｊ. Ｍ. Jun 9 '13 at 14:31
I use Carlson's symmetric integrals; this is important because the expressions in terms of Legendre's integrals involve big cancellations. The use of an imaginary modulus is a non-issue once the integrals have been converted to Carlson's forms. I don't find the increased error too surprising. With the series solution, round-off can be very well controlled. It is difficult to maintain this accuracy with the more general formulation (however, I admit that I didn't try really hard). In any case, the accuracy is still very good -- about 25 nanometers for WGS84. –  cffk Jun 9 '13 at 17:19

I've completed the implementation of my geodesic algorithms in Java. Documentation (including download information) is available at

http://geographiclib.sf.net/html/java

I gather that it's easy to call Java code from Mathematica. If someone would be kind enough to post instructions on how to do this, I would be happy to include this in the documentation.

This implementation uses series expanded to 6th order in the flattening. This gives results accurate to round-off (using double precision) for |f| < 1/50. There are two possible extensions which may be of interest to Mathematica users:

1. The solution in terms of elliptic integrals (to deal with highly eccentric ellipsoids of revolution); this is currently implemented in the C++ and Maxima versions.

2. The use of arbitrary precision arbitrary precision arithmetic; this is currently implemented in the Maxima version.

-