[Proj] Problems with Pittman geodesic??

[Karney, Charles]

> From: Gerald I. Evenden [geraldi.evenden at gmail.com]
> Sent: Sunday, January 18, 2009 20:57
> To: PROJ.4 and general Projections Discussions
> Subject: [Proj] Problems with Pittman geodesic??
>
> While checking out the accuracy of Vincenty vs. Pittman I was fairly
> consistently getting agreement to micron or better.
>
> BUT not always!!!
>
> I would appreciate anyone who has a Pittman geodesic routine to test
> the following inverse (WGS84 ellipsoid):
>
> point 1 at 0 lat, 0 lon and point 2 at 45N 90E.
>
> With my version of Pittman I get a distance of 9993541.5348708 AND I
> get the same distance while changing the longitude of the second point
> from about 89.6 to a hair over 90.
>
> At 0,0/45,89.5:
> Pittman: 9970963.0100082
> Vincenty:9970963.01000801
>
> At 0,0/45,89.6
> Pittman: 9993541.5348708
> Vincenty:9978847.65947167
>
> Pittman remained constant in this longitude interval
>
> At 0,0/45,90
> Pittman:  9993541.5348708
> Vincenty:10010386.3610382
>
> Between 90.0000001 and 90.000001 Pittman finally came back into near
> agreement with Vincenty.
>
> Note: I did nothing to the Pittman FORTRAN loaned from Mugnier other
> than compile it with gfortran and link to a C driver.

Gerald:

I've coded up Pittman's algorithm (as given in his 1986 paper) in C++.
I left the order alone (N = 5), but increased the number of allowed
iterations in the inverse method to 100.  I followed the Fortran code
fairly slavishly, but used long doubles for double precision.  (Pittman
was on a machine where 1.0d0 - 1.0d-18 differed from 1.  In addition,
there are numerous places in his algorithm where the method is sensitive
to roundoff; so I figured I'd provide the extra 11 bits of precision to
be safe.)  Here are the results I get for the inverse calculation (for
the WGS84 ellipsoid).  Lat1 = lon1 = 0, lat2 = 45; the other columns are

lon2 azi1 azi2 s12 niter

89.0 45.093504806662311 89.443934945243591  9931540.5187563721 6
89.1 45.094149622603338 89.514653188105030  9939424.8761606768 6
89.2 45.094706860005074 89.585369717060000  9947309.3159994913 7
89.3 45.095176523124121 89.656084748302839  9955193.8262614792 7
89.4 45.095558615682521 89.726798498013506  9963078.3949346433 7
89.5 45.095853140867789 89.797511182358922  9970963.0100083940 7
89.6 45.096212150579780 90.000000000000000  9993541.5243879881 100
89.7 45.096212150579780 90.000000000000000  9993541.5243879881 100
89.8 45.096212150579780 90.000000000000000  9993541.5243879881 100
89.9 45.096212150579780 90.000000000000000  9993541.5243879881 100
90.0 45.096012330347754 90.151066188700727 10010386.3610382384 10
90.1 45.095781488302975 90.221777019794016 10018271.0023121920 7
90.2 45.095463086233842 90.292488298435018 10026155.6059209196 7
90.3 45.095057123052886 90.363200240733359 10034040.1598547278 7
90.4 45.094563597138410 90.433913062797377 10041924.6521030197 6
90.5 45.093982506334513 90.504626980735483 10049809.0706572333 6

This confirms your observation and offers an explanation.  Pittman's
inverse method fails to converge in the interval lon2 in [89.6, 89.9]
where the returned distance is constant.

I was going to spend another day checking over my transcription of the
code; but I figured that, given you see the same behavior, the problem
is likely to be in the algorithm.

The iterative method Pittman uses for the inverse is Newton's method
where he has substituted an approximate expression for the derivative.
Unsurprisingly, this sometimes fails to converge.  Another possible
limitation of Pittman's method is that the code needs to know the number
of vertices between the endpoints.  As the iteration proceeds, the
estimate of the number of vertices isn't updated but should be(?).  This
would explain why problems occur near azi2 = 90.

I've also found a couple of bugs in the the published code.
(Uninitialized variables are accessed in a couple of places; the guards
against taking square roots of negative numbers are not always correct.)
However, it's possible that the version of the code you got from Mugnier
has these fixed.

--
Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

URL: http://charles.karney.info
Tel: +1 609 734 2312
Fax: +1 609 734 2662