> From: Gerald I. Evenden [

[hidden email]]

> 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 <

[hidden email]>

Sarnoff Corporation, Princeton, NJ 08543-5300

URL:

http://charles.karney.infoTel: +1 609 734 2312

Fax: +1 609 734 2662

_______________________________________________

Proj mailing list

[hidden email]
http://lists.maptools.org/mailman/listinfo/proj