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