[Proj] Re: Dozier's TM method---my summary

[Gerald I. Evenden]

On Wednesday 05 July 2006 8:31 am, Oscar van Vlijmen wrote:
> > From: "Gerald I. Evenden"
> >
> > Stating it another way, my surrender was to the Dozier article that while
> > the method sounds interesting the execution falls short.  The code had
> > several errors and some methodology is in question.  I have spent an
> > additional week rooting out a complex Bulirsch routine from the old Bell
> > lab site and spent several days relating Jacobi's form to Legendre in
> > order to test against Abramowitz tables and GSL software for at least
> > real arguments---they now agree.  But Dozier's code does not agree to
> > better than the third or forth place (I am not saying Dozier's code is
> > wrong yet).
> >
> > And I still have no handle on a replacement for Newton-Raphson.
>
> What surprises me is that I gave an alternative twice, once in a posting on
> this list dated 2006-06-13, once in a private email dated 2006-06-23,
> including a web reference.
> It's about ACM TOMS algorithm 365 from H. Bach. Fortran code can be found
> at netlib.org and one can buy a copy of the peer reviewed TOMS article at
> ACM (acm.org).
> I am not happy with the method because it is horribly slow, but it is a
> very robust method, at least if one knows how to tweak its control
> parameters.
>
> As for a replacement of complex Jacobi elliptic functions and a complex
> elliptic integral: it can be done without too much effort. I provided in
> the same manner (list posting and private email) references.
> I do admit that working with complex elliptics is no trivial endeavor
> without knowledge or experience.

The Bulirsch code seems to work and compares in the real case (x+0i) with the 
GSL library *and* Abramowich tables--- after a lot of footwork dealing with 
figuring out how to match up the (&*^(&^ parameters to each.  The input 
parameters were what was really blowing my mind.

I ran into another, perhaps the same source, FORTRAN source but found out I 
did not have a g77 compiler (I had recently  changed versions of Linux) and 
solving that problem introduced more headaches than I wanted at the moment.

In my little bench program where I compared GSL, Bulirsch and Dozier I ran 
into consistant disagreement with Dozier code in the 3rd or 4th decimal 
place.  At the moment, this does not make sense because Dozier's code is at a 
significant point in Dozier's second phase of the forward solution I do not 
understand how it could be incorrect and yet yield correct answers in many 
areas for the overall projection.

The following is a short example of a test run:

e^2 used: 0.006722670022

x + y*i input: 0.1 0
Dozier e2ci: 0.09999998496533813  0*i
Bulirsch elco2 ret: OK
elco2: 0.0999988817825081  0*i
gsl_sf_ellint_E_e x: 0.09999888178250819 --- err: 4.44104e-17

x + y*i input: 0.2 0
Dozier e2ci: 0.19999988044241  0*i
Bulirsch elco2 ret: OK
elco2: 0.1999911075211161  0*i
gsl_sf_ellint_E_e x: 0.1999911075211167 --- err: 8.88297e-17

x + y*i input: 0.5 0
Dozier e2ci: 0.4999982088519777  0*i
Bulirsch elco2 ret: OK
elco2: 0.4998667513591754  0*i
gsl_sf_ellint_E_e x: 0.4998667513591855 --- err: 2.22222e-16

x + y*i input: 0.3 0.3
Dozier e2ci: 0.3000007837129038  0.2999991577227006*i
Bulirsch elco2 ret: OK
elco2: 0.3000583008002545  0.2999373655430532*i

In the last case, GSL not run because y != 0.

That gives you an idea where I now stand with the complex elliptical integral 
of the second kind (incomplete).

The nasty part is that even if this is solved, the problem with the first 
phase and the root finding looms as an even bigger issue.

Sigh.  :-(
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj

-- 
Jerry and the low-riders: Daisy Mae and Joshua
"Cogito cogito ergo cogito sum"
   Ambrose Bierce, The Devil's Dictionary