[PROJ] Geodetic to Authalic latitude conversions

[Charles Karney]

Surely you should use Horner to compute the polynomials.  This gives you
tighter code, lets you vary the order of the polynomial and is more
accurate (small terms accumulated first).

Likewise, use Clenshaw for the trig sum.  All the same advantages accrue
+ it only needs one or two evaluations of trig functions.

My arxiv paper has been published in Survey Review

  https://doi.org/10.1080/00396265.2023.2217604

However, retain the arxiv link since the journal article is behind a
paywall.

On Tue, Sep 10, 2024 at 9:19 PM Jérôme St-Louis via PROJ
<proj at lists.osgeo.org> wrote:
>
> Hi Even, all,
>
> As suggested in https://github.com/OSGeo/PROJ/pull/4211 , I've just written an implementation of the authalic / geodetic latitude conversion with the goal to automatically perform this conversion when using ISEA with an non-spherical ellipsoid, and was wondering whether we would want to improve the PROJ version (pj_authlat() for the geodetic ==> authalic), which is used for ISEA, but also by:
>
> - Equal Area Cylindrical (cea.cpp)
> - EqualEarth (eqearth.cpp)
> - HEALPix (healpix.cpp) -- which might also imply rHEALPix
> - Lambert Azimuthal Equal Area (laea.cpp)
>
> The current version of pj_authlat() uses a series expansion based on the eccentricity squared (I believe this is based on a Lagrange reversion as mentioned here), and Charles Karney in this paper establishes that:
>
> series expansions using the third flattening as the small parameter uniformly result in smaller truncation errors compared to using the eccentricity.
>
> The current version also only uses order 3, whereas GeographicLib would normally use order 6 for double precision, with order 4 being the minimum configuration.
> The current order 3, which requires 2 sin() calls, gives precision up to roughly 8 decimals.
> Order 6 would mean twice the number of sin() calls, and the 8 decimals should be enough to point to Waldo on a page, however C.K. argues that:
>
> Modern geodetic libraries should strive to achieve full double-precision accuracy (page 4)
>
> and I would much agree with that.
>
> Based on the direction we want to go, I would submit a PR including adjustments to avoid the separate heap allocation for the doubles (as discussed here), and adjustments for tests etc.
>
> PROJ currently does not use the full GeographicLib which already includes code for conversion between different auxiliary latitudes, does it?
>
>     ( https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp )
>
> See my proposed new version below, which I wrote based on equation 20 on page 3, and matrix A20 on page 13 from the paper, applying optimizations similar to the current pj_authlat().
>
> Kind regards,
>
> -Jerome
>
> void pj_authset(double a, double b, double cp[6])
> {
>    // https://arxiv.org/pdf/2212.05818
>    // ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1)    -- (20)
>    double n = (a - b) / (a + b);  // Third flattening
>    static const double C[21] = // (A20)
>    {
>       4 / 3.0,  4 / 45.0,   -16/35.0,  -2582 /14175.0,  60136 /467775.0,    28112932/ 212837625.0,
>                46 / 45.0,  152/945.0, -11966 /14175.0, -21016 / 51975.0,   251310128/ 638512875.0,
>                          3044/2835.0,   3802 /14175.0, -94388 / 66825.0,    -8797648/  10945935.0,
>                                         6059 / 4725.0,  41072 / 93555.0, -1472637812/ 638512875.0,
>                                                        768272 /467775.0,  -455935736/ 638512875.0,
>                                                                           4210684958/1915538625.0
>    };
>    double p;
>
>    cp[0]  = C[ 0] * n;
>
>    p = n * n;
>    cp[0] += C[ 1] * p;
>    cp[1]  = C[ 6] * p;
>
>    p *= n;
>    cp[0] += C[ 2] * p;
>    cp[1] += C[ 7] * p;
>    cp[2]  = C[11] * p;
>
>    p *= n;
>    cp[0] += C[ 3] * p;
>    cp[1] += C[ 8] * p;
>    cp[2] += C[12] * p;
>    cp[3]  = C[15] * p;
>
>    p *= n;
>    cp[0] += C[ 4] * p;
>    cp[1] += C[ 9] * p;
>    cp[2] += C[13] * p;
>    cp[3] += C[16] * p;
>    cp[4]  = C[18] * p;
>
>    p *= n;
>    cp[0] += C[ 5] * p;
>    cp[1] += C[10] * p;
>    cp[2] += C[14] * p;
>    cp[3] += C[17] * p;
>    cp[4] += C[19] * p;
>    cp[5]  = C[20] * p;
> }
>
> double pj_authlat(double auth, const double * cp)
> {
>    // https://arxiv.org/pdf/2212.05818
>    // ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1)    -- (20)
>    double a2x = auth + auth, a4x = a2x + a2x, a8x = a4x + a4x;
>    return auth +
>       cp[0] * sin(a2x) +
>       cp[1] * sin(a4x) +
>       cp[2] * sin(a4x + a2x) +
>       cp[3] * sin(a8x) +
>       cp[4] * sin(a8x + a2x) +
>       cp[5] * sin(a8x + a4x);
> }
>
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj



-- 
Charles Karney <karney at alum.mit.edu>
702 Prospect Ave
Princeton, NJ 08540