Problems with Pittman geodesic??

classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

Problems with Pittman geodesic??

Gerald I. Evenden
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.

--
The whole religious complexion of the modern world is due
to the absence from Jerusalem of a lunatic asylum.
-- Havelock Ellis (1859-1939) British psychologist
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Problems with Pittman geodesic??

hamish-2
Gerald wrote:

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


I would guess that some internal variable is defined as a REAL (float)
instead of a DOUBLE PRECISION (double), and so the effect you are seeing
is quantization due to variable precision.


Hamish



     

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Problems with Pittman geodesic??

Gerald I. Evenden
On Sunday 18 January 2009 11:00:36 pm Hamish wrote:

> Gerald wrote:
> > 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.
>
> I would guess that some internal variable is defined as a REAL (float)
> instead of a DOUBLE PRECISION (double), and so the effect you are seeing
> is quantization due to variable precision.

Both replies, so far seemed to think I was trying to debug the Pittman
procedures.

I am not interested in fixing the Pittman procedures.  I am asking for anyone
with the Pittman routines to check if their version has the same problem mine
has.  My only interest in Pittman is as a comparison to the Vincenty
procedures because Pittman has been touted as the standard for geodesic
calculations.

Thank-you.

--
The whole religious complexion of the modern world is due
to the absence from Jerusalem of a lunatic asylum.
-- Havelock Ellis (1859-1939) British psychologist
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Problems with Pittman geodesic??

OvV_HN
In reply to this post by Gerald I. Evenden
Some time ago I referred to 2 articles by Saito.
He computes the geodesic problems by numerical integration. There are some
errors in the formulae, but his second article gives correct results for
'normal' values. Saito attends much detail to borderline situations, but
there are probably some errors in his descriptions.
Anyhow, he gives a couple of test points which everybody with geodesics code
definitively should try to replicate. You might be surprised!

Saito article 2 examples.

Using the formulas and the procedure described in the article, with an
equatorial radius of 6378388 meters and an inverse flattening of 297, the
following results were obtained. The computation was made in double
precision arithmetic - through twenty significant digits.

Line lat1                lat2               delta lon
1. 37d 19' 54.95367"   26d 07' 42.83946"  041d 28' 35.50729"
2. 35  16  11.24862    67  22  14.77638   137  47  28.31435
3. 01  00  00.00000   -00  59  53.83076   179  17  48.02997
4. 01  00  00.00000    01  01  15.18952   179  46  17.84244
5. 41  41  45.88000    41  41  46.20000   000  00  00.56000
6. 30  00  00.00000    37  53  32.46584   116  19  16.68843
7. 30  19  54.95367   -30  11  50.15681   179  58  17.84244
8. 00  39  49.12586   -00  45  14.13112   179  58  17.84244
9. 00  00  54.95367    00  00  42.83946   179  28  17.84244

Line  dist.           azim.                back azim.
1.  4085966.7026 m   095d 27' 59.630899"  118d 05' 58.961609"
2.  8084823.8383     015  44  23.748498   114  55  39.921473
3. 19959999.9998     088  59  59.998971   091  00  06.118356
4. 19780006.5588     004  59  59.999957   174  59  59.884800
5.       16.2839751  052  40  39.390671   052  40  39.763172
6. 10002499.9999     045  00  00.000004   129  08  12.326009
7. 19989590.5480     002  23  52.108130   177  36  19.670109
8. 19994529.4446     177  39  39.010104   002  20  21.153036
9. 19977290.7711     054  08  27.731619   125  51  32.272327

From:
THE COMPUTATION OF LONG GEODESICS ON THE ELLIPSOID
THROUGH GAUSSlAN QUADRATURE
Tsutomu SAITO
Construction College,
Kodaira-shi Tokyo, 187 Japan.


Oscar van Vlijmen


----- Original Message -----
From: "Gerald I. Evenden" <[hidden email]>
To: "PROJ.4 and general Projections Discussions" <[hidden email]>
Sent: Monday, January 19, 2009 2:57 AM
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):
>

.....


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Problems with Pittman geodesic??

Gerald I. Evenden
On Monday 19 January 2009 1:16:38 pm OvV_HN wrote:
> Some time ago I referred to 2 articles by Saito.
> He computes the geodesic problems by numerical integration. There are some
> errors in the formulae, but his second article gives correct results for
> 'normal' values. Saito attends much detail to borderline situations, but
> there are probably some errors in his descriptions.
> Anyhow, he gives a couple of test points which everybody with geodesics
> code definitively should try to replicate. You might be surprised!

I did not find any surprises here because where problems do show up are for
the near apodal points which Vincenty (and others) usually flunk.

I have more problems with failures with Pittman and some non-apodal points.

> Saito article 2 examples.
>
> Using the formulas and the procedure described in the article, with an
> equatorial radius of 6378388 meters and an inverse flattening of 297, the
> following results were obtained. The computation was made in double
> precision arithmetic - through twenty significant digits.

"double precision" --- "through twenty significant digits." ???  What kind of
hardware is being used?  "double" on most 64bit systems is about 15 digits
+-.  With so much significance, why are the following tables showing 12 at
best?

> Line lat1                lat2               delta lon
> 1. 37d 19' 54.95367"   26d 07' 42.83946"  041d 28' 35.50729"
> 2. 35  16  11.24862    67  22  14.77638   137  47  28.31435
> 3. 01  00  00.00000   -00  59  53.83076   179  17  48.02997
> 4. 01  00  00.00000    01  01  15.18952   179  46  17.84244
> 5. 41  41  45.88000    41  41  46.20000   000  00  00.56000
> 6. 30  00  00.00000    37  53  32.46584   116  19  16.68843
> 7. 30  19  54.95367   -30  11  50.15681   179  58  17.84244
> 8. 00  39  49.12586   -00  45  14.13112   179  58  17.84244
> 9. 00  00  54.95367    00  00  42.83946   179  28  17.84244
>
> Line  dist.           azim.                back azim.
> 1.  4085966.7026 m   095d 27' 59.630899"  118d 05' 58.961609"
> 2.  8084823.8383     015  44  23.748498   114  55  39.921473
> 3. 19959999.9998     088  59  59.998971   091  00  06.118356
> 4. 19780006.5588     004  59  59.999957   174  59  59.884800
> 5.       16.2839751  052  40  39.390671   052  40  39.763172
> 6. 10002499.9999     045  00  00.000004   129  08  12.326009
> 7. 19989590.5480     002  23  52.108130   177  36  19.670109
> 8. 19994529.4446     177  39  39.010104   002  20  21.153036
> 9. 19977290.7711     054  08  27.731619   125  51  32.272327
>
> From:
> THE COMPUTATION OF LONG GEODESICS ON THE ELLIPSOID
> THROUGH GAUSSlAN QUADRATURE
> Tsutomu SAITO

Using Gaussian Quadrature is not a big thing as long as one has the integrals
to evaluate---one does not need to use gsl to do that. ;-)  I am surprised
that he can get away with just using that method---or does he?

Lastly, the Saito source is apparently available though "Springer --
something" that charges $35(presumably US) for a pdf copy of Vincenty that
you can get for free from NOAA.

--
The whole religious complexion of the modern world is due
to the absence from Jerusalem of a lunatic asylum.
-- Havelock Ellis (1859-1939) British psychologist
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Problems with Pittman geodesic??

OvV_HN
From: "Gerald I. Evenden" <[hidden email]>
Sent: Thursday, January 22, 2009 9:40 PM
Subject: Re: [Proj] Problems with Pittman geodesic??


> On Monday 19 January 2009 1:16:38 pm OvV_HN wrote:
>> Some time ago I referred to 2 articles by Saito.
....
>> Anyhow, he gives a couple of test points which everybody with geodesics
>> code definitively should try to replicate. You might be surprised!
> I did not find any surprises here because where problems do show up are
> for
> the near apodal points which Vincenty (and others) usually flunk.

Nice. Indeed my implementation of NOAAs Vincenty went wrong with the
antipodal examples 7, 8, 9. As did the earlier mentioned method of Xue-Lian
Zhang.

>> Saito article 2 examples.
...
>> following results were obtained. The computation was made in double
>> precision arithmetic - through twenty significant digits.
> "double precision" --- "through twenty significant digits." ???  What kind
> of
> hardware is being used?

The article is from 1978; only Fortran is mentioned.

> Using Gaussian Quadrature is not a big thing as long as one has the
> integrals
> to evaluate---one does not need to use gsl to do that. ;-)

Then you probably don't need GNU GSL for the Mayr-Tobler projection!

Oscar van Vlijmen




_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Problems with Pittman geodesic??

Karney, Charles
In reply to this post by Gerald I. Evenden
> 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.info
Tel: +1 609 734 2312
Fax: +1 609 734 2662
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Problems with Pittman geodesic??

Clifford J Mugnier
Re: [Proj] Problems with Pittman geodesic??
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.
 
C. Mugnier
LSU


From: [hidden email] on behalf of Karney, Charles
Sent: Mon 16-Mar-09 16:12
To: [hidden email]; PROJ.4 andgeneral Projections Discussions
Subject: Re: [Proj] Problems with Pittman geodesic??

> 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.info
Tel: +1 609 734 2312
Fax: +1 609 734 2662
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj