Improved series for the meridian length and its inverse

[cffk]

This PR addresses #2387 point 3.

Replace the series for the meridian length in the functions

pj_enfn
pj_mlfn
pj_inv_mlfn

in mlfn.cpp with improved versions. The old version of pj_mlfn used Snyder's Eq. (3-21) extended from 3rd to 4th order in e^2 and pj_inv_mlfn inverted this by an iterative method.

The new version uses 6th order expansions in the third flattening, n, for both pj_mlfn and pj_inv_mlfn. Expansions in n have a smaller truncation error that the corresponding expansion in e^2. Increasing the order to 6th (from 4th) means that the routines give full double precision accuracy for |f| <= 1/150. Because the meridian length is used in several projections, it's important that it be accurate as possible. For example, now the azimuthal equidistant projection is equally accurate in the polar (uses pj_mlfn) and oblique (uses the geodesic routines) aspects.

The increased accuracy tripped one of the tests: an instance of the projection +proj=utm +approx. The fix is to use the non-approx results for the expected values.

Presumably this change results in pj_mlfn being a little slower while pj_inv_mlfn is considerably faster (since there's no iteration involved). Incidentally the new method for pj_inv_mlfn uses Snyder's recommended series Eq. (3-23) extended from 4th to 6th order (Snyder writes n as e_1).

More discussion of the series used here is given in my recent paper

On auxiliary latitudes, https://arxiv.org/abs/2212.05818

• Partially closes Some suggestions for code cleanup #2387
• Tests adjusted
• Added clear title that can be used to generate release notes

[rouault]

Can you try running the new bench_proj_trans binary before and after your changes ?
I get for example the following with current master:

$ bin/bench_proj_trans --pipeline "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone=31 +ellps=GRS80 +approx" 2 49
2 49 -> 426857.98771659 5427937.52334353
Duration: 407 ms
Throughput: 12.29 million coordinates/s

$ bin/bench_proj_trans --pipeline "+proj=pipeline +step +inv +proj=utm +zone=31 +ellps=GRS80 +approx +step +proj=unitconvert +xy_in=rad +xy_out=deg" 426857.98771659 5427937.52334353
426857.98771659 5427937.52334353 -> 1.99999999999991 48.9999999999999
Duration: 714 ms
Throughput: 7.00 million coordinates/s

[cffk]

Here you go:

BEFORE
2 49 -> 426857.98771659 5427937.52334353
Duration: 305 ms
Throughput: 16.39 million coordinates/s
426857.98771659 5427937.52334353 -> 1.99999999999991 48.9999999999999
Duration: 546 ms
Throughput: 9.16 million coordinates/s

AFTER
2 49 -> 426857.98771659 5427937.52334229
Duration: 333 ms
Throughput: 15.02 million coordinates/s
426857.98771659 5427937.52334353 -> 1.99999999999969 49.000000000011
Duration: 496 ms
Throughput: 10.08 million coordinates/s

These timings are roughly in line with my expectations.

[rouault]

I've given this a try against GDAL autotest suite whose one test is sensitive to +proj=utm +approx results, and things are good.
It is also very cool to have a non-iterative implementation of inv_mlfn, in the perspective of PROJ having a day vectorized implementations of coordinate transformations, transforming several coordinates at once in SIMD registers.
So +1 from me

[kbevers]

Looks good to me. Great addition!