<HTML dir=ltr><HEAD><TITLE>Re: [Proj] Problems with Pittman geodesic??</TITLE>
<META http-equiv=Content-Type content="text/html; charset=unicode">
<META content="MSHTML 6.00.6000.16809" name=GENERATOR></HEAD>
<BODY>
<DIV id=idOWAReplyText70196 dir=ltr>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=2>I do recall that Pittman had several telephone conversations with Thaddeus Vincenty before T.V. passed away. I recall that Pittman indeed made some modifications to his published Fortran code as a result of his conversations. I do not know or recall what they were other than it was sometbing about singularities. The compiler Pittman used was Lahey Fortran 77; I don't recall which version, but it was the current one as of his paper.</FONT></DIV>
<DIV dir=ltr><FONT size=2></FONT> </DIV>
<DIV dir=ltr><FONT size=2>C. Mugnier</FONT></DIV>
<DIV dir=ltr><FONT size=2>LSU</FONT></DIV></DIV>
<DIV dir=ltr><BR>
<HR tabIndex=-1>
<FONT face=Tahoma size=2><B>From:</B> proj-bounces@lists.maptools.org on behalf of Karney, Charles<BR><B>Sent:</B> Mon 16-Mar-09 16:12<BR><B>To:</B> geraldi.evenden@gmail.com; PROJ.4 andgeneral Projections Discussions<BR><B>Subject:</B> Re: [Proj] Problems with Pittman geodesic??<BR></FONT><BR></DIV>
<DIV>
<P><FONT size=2>> From: Gerald I. Evenden [geraldi.evenden@gmail.com]<BR>> Sent: Sunday, January 18, 2009 20:57<BR>> To: PROJ.4 and general Projections Discussions<BR>> Subject: [Proj] Problems with Pittman geodesic??<BR>><BR>> While checking out the accuracy of Vincenty vs. Pittman I was fairly<BR>> consistently getting agreement to micron or better.<BR>><BR>> BUT not always!!!<BR>><BR>> I would appreciate anyone who has a Pittman geodesic routine to test<BR>> the following inverse (WGS84 ellipsoid):<BR>><BR>> point 1 at 0 lat, 0 lon and point 2 at 45N 90E.<BR>><BR>> With my version of Pittman I get a distance of 9993541.5348708 AND I<BR>> get the same distance while changing the longitude of the second point<BR>> from about 89.6 to a hair over 90.<BR>><BR>> At 0,0/45,89.5:<BR>> Pittman: 9970963.0100082<BR>> Vincenty:9970963.01000801<BR>><BR>> At 0,0/45,89.6<BR>> Pittman: 9993541.5348708<BR>> Vincenty:9978847.65947167<BR>><BR>> Pittman remained constant in this longitude interval<BR>><BR>> At 0,0/45,90<BR>> Pittman: 9993541.5348708<BR>> Vincenty:10010386.3610382<BR>><BR>> Between 90.0000001 and 90.000001 Pittman finally came back into near<BR>> agreement with Vincenty.<BR>><BR>> Note: I did nothing to the Pittman FORTRAN loaned from Mugnier other<BR>> than compile it with gfortran and link to a C driver.<BR><BR>Gerald:<BR><BR>I've coded up Pittman's algorithm (as given in his 1986 paper) in C++.<BR>I left the order alone (N = 5), but increased the number of allowed<BR>iterations in the inverse method to 100. I followed the Fortran code<BR>fairly slavishly, but used long doubles for double precision. (Pittman<BR>was on a machine where 1.0d0 - 1.0d-18 differed from 1. In addition,<BR>there are numerous places in his algorithm where the method is sensitive<BR>to roundoff; so I figured I'd provide the extra 11 bits of precision to<BR>be safe.) Here are the results I get for the inverse calculation (for<BR>the WGS84 ellipsoid). Lat1 = lon1 = 0, lat2 = 45; the other columns are<BR><BR>lon2 azi1 azi2 s12 niter<BR><BR>89.0 45.093504806662311 89.443934945243591 9931540.5187563721 6<BR>89.1 45.094149622603338 89.514653188105030 9939424.8761606768 6<BR>89.2 45.094706860005074 89.585369717060000 9947309.3159994913 7<BR>89.3 45.095176523124121 89.656084748302839 9955193.8262614792 7<BR>89.4 45.095558615682521 89.726798498013506 9963078.3949346433 7<BR>89.5 45.095853140867789 89.797511182358922 9970963.0100083940 7<BR>89.6 45.096212150579780 90.000000000000000 9993541.5243879881 100<BR>89.7 45.096212150579780 90.000000000000000 9993541.5243879881 100<BR>89.8 45.096212150579780 90.000000000000000 9993541.5243879881 100<BR>89.9 45.096212150579780 90.000000000000000 9993541.5243879881 100<BR>90.0 45.096012330347754 90.151066188700727 10010386.3610382384 10<BR>90.1 45.095781488302975 90.221777019794016 10018271.0023121920 7<BR>90.2 45.095463086233842 90.292488298435018 10026155.6059209196 7<BR>90.3 45.095057123052886 90.363200240733359 10034040.1598547278 7<BR>90.4 45.094563597138410 90.433913062797377 10041924.6521030197 6<BR>90.5 45.093982506334513 90.504626980735483 10049809.0706572333 6<BR><BR>This confirms your observation and offers an explanation. Pittman's<BR>inverse method fails to converge in the interval lon2 in [89.6, 89.9]<BR>where the returned distance is constant.<BR><BR>I was going to spend another day checking over my transcription of the<BR>code; but I figured that, given you see the same behavior, the problem<BR>is likely to be in the algorithm.<BR><BR>The iterative method Pittman uses for the inverse is Newton's method<BR>where he has substituted an approximate expression for the derivative.<BR>Unsurprisingly, this sometimes fails to converge. Another possible<BR>limitation of Pittman's method is that the code needs to know the number<BR>of vertices between the endpoints. As the iteration proceeds, the<BR>estimate of the number of vertices isn't updated but should be(?). This<BR>would explain why problems occur near azi2 = 90.<BR><BR>I've also found a couple of bugs in the the published code.<BR>(Uninitialized variables are accessed in a couple of places; the guards<BR>against taking square roots of negative numbers are not always correct.)<BR>However, it's possible that the version of the code you got from Mugnier<BR>has these fixed.<BR><BR>--<BR>Charles Karney <ckarney@sarnoff.com><BR>Sarnoff Corporation, Princeton, NJ 08543-5300<BR><BR>URL: <A href="http://charles.karney.info/">http://charles.karney.info</A><BR>Tel: +1 609 734 2312<BR>Fax: +1 609 734 2662<BR>_______________________________________________<BR>Proj mailing list<BR>Proj@lists.maptools.org<BR><A href="http://lists.maptools.org/mailman/listinfo/proj">http://lists.maptools.org/mailman/listinfo/proj</A><BR></FONT></P></DIV></BODY></HTML>