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@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@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/proj -- Charles Karney <kar...@alum.mit.edu> 702 Prospect Ave Princeton, NJ 08540 _______________________________________________ PROJ mailing list PROJ@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/proj
