[Proj] wrong computation of meridinal distance / pj mlfn.c

[Gerald I. Evenden]

A bit more on the meridional distance (pardon my continual misspelling).

The kernel of the integral in the manual is simply the radius of the
earth's surface in
the plane of the meridian:

rho = a(1-e^2)/(1-e^2sin^2(phi))^(3/2)

where a:earth's semi-major axis, e: eccentricity and phi: geodetic
latitude.
Note that if eccentricity becomes 0 then then rho = a and integrating
for
meridional distance simply becomes a*phi.

The traditional solution has been to do a binomial expansion of the
denominator
and follow with a term-by-term integration of the series. Very messy and
keeps
one out of the bars for a while.  The difference between the multiple
angle
series (sin(n*phi)) versus the power series ((sin(phi))^n) depends upon
which set
of integral tables one wishes to use.  The power series is
computationally
desirable in minimizing computer time by drastically cutting time spent
calling
the sine function (another power series).

The elliptic integral method (which is also a power series) came about
when
I happened to locate the integral in Gradsbteyn & Rysbik along with a
method
to evaluate the elliptic integral that gave solution which converged to
machine
precision fairly rapidly.  The method is slightly slower than the last
power
series method in PROJ4 but much more accurate.

Testing was made by numerical integration with GP-PARI multiple
precision package
although any multiple precision package ought to work.
-- 
_____________________________________________________________
Jerry and the Low Riders: Daisy Mae and Joshua
"The most certain test by which we judge whether a country is
really free is the amount of security enjoyed by minorities"
---Lord Acton, 1907 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20050713/dc4bbb13/attachment.html