[Proj] Re: libproj4 stmerc = French Gauss-Laborde projection

[Strebe at aol.com]

I heartily recommend the AGM method for computing the elliptic integral. It's 
fast and accurate everywhere, and can be implemented to arbitrary precision. 
Polynomials seem a poor choice for this.

See Abramowitz & Stegun, p. 598 (17.6). Bulirsch (1965) supplies an excellent 
implementation of the AGM, but you can use it as described in Abramowitz & 
Stegun with great efficiency.

Regards
-- daan Strebe


In a message dated 6/22/06 00:39:37, martin.vermeer at hut.fi writes:


> On Thu, 2006-06-15 at 12:27 +0300, Martin Vermeer wrote:
> > On Thu, 2006-06-15 at 05:05 -0400, Strebe at aol.com wrote:
> > >
> > > I'm not lost. I've analyzed the method and programmed it. It works
> > > fine.
> > >
> > > Yes, p is the parametric co-latitude in the first formula. It's wrong
> > > to use the same variable in the later formulae because p refers to
> > > something else there -- in fact, it's a complex variable in the later
> > > instances.
> > >
> > > The "trick" is this: Wallis uses the polar stereographic because it's
> > > the simplest way to get a conformal mapping to the plane. Once the
> > > ellipsoid is mapped, he treats the plane as the complex plane and
> > > looks for a complex "co-latitude" which can be used with the polar
> > > stereographic, but this time treating the polar stereographic as
> > > function of a complex variable. The reason he does this is (a) to
> > > preserve conformality; and (b) so that the central meridian (in which
> > > the imaginary axis is 0) maps back to the parametric colatitude. At
> > > this point the ellipsoid is mapped conformally in such a way that
> > > leaves the central meridian effectively unmapped.
> > >
> > > Leaving the complex plane aside, using the colatitude as the parameter
> > > to the elliptic integral of the second kind gives the distance from
> > > the pole to the colatitude. Since this odd mapping Wallis contrived
> > > effectively leaves the central meridian unmapped, and since any
> > > analytic function applied to a conformal mapping results in a
> > > conformal mapping, and since the elliptic integral has an analytic
> > > form, all that is left is to push the mapping through the complex form
> > > of the elliptic integral of the second kind. This "straightens out"
> > > the central meridian to its true differential lengths whilst dragging
> > > the whole complex plane with it in a conformal fashion. The result
> > > must be the transverse Mercator, since the central meridian is
> > > projected with correct scale and since a conformal projection is
> > > unique except with respect to scale and rotation.
> > >
> > > Regards,
> > > -- daan Strebe
> > >
> >
> > Ah! But that's precisely what I have been doing numerically, using a
> > polynomial expansion rather than an elliptic integral! (And starting
> > from classical Mercator rather than stereographic, so it will run into
> > problems at high latitudes.) It's essentially solving a boundary value
> > problem, with the set of PDEs being the Cauchy-Riemann conformity
> > conditions and the central meridian the boundary.
> >
> > I suppose I have to get it written up in english ;-)
> 
> 
> Just found my Matlab routines for playing with this. Posted them here:
> 
> http://users.tkk.fi/~mvermeer/gk.zip
> 
> You run it as
> 
> >> gkmain(phi, lab, [phimin phimax delta maxk refphi]);
> 
> e.g.
> 
> >> gkmain(60, 5, [55 65 0.5 12 60])
> 
> (lab is the distance from the central meridian)
> 
> With this I have also obtained sub-millimeter precision (in typical
> situations). Yes, it breaks down at the poles due to using Mercator as
> the starting projection. Replacing this by stereographic as Dr Wallis
> did, shouldn't be hard.
> 
> Using a polynomial expansion (instead of the elliptic integral /
> Newton-Raphson thing) has the advantage of simplicity in implementation.
> One disadvantage that I noticed is, that the normal matrix for
> estimating the coefficients gets poorly conditioned for expansion degree
> maxk = 15 or more.
> 
> I used this code in an exercise for my students; in a production
> environment you would print out and hardwire the estimated coefficients
> for the area of interest into your code.
> 
> The exercise instruction (in Finnish :) at
> 
> http://users.tkk.fi/~mvermeer/kotiharjB.pdf
> 
> 
> - Martin
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20060622/25ebc4aa/attachment.html>