<!DOCTYPE html>
<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hi Even, all,</p>
<p>As suggested in <a class="moz-txt-link-freetext" href="https://github.com/OSGeo/PROJ/pull/4211">https://github.com/OSGeo/PROJ/pull/4211</a> , 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 (<i>pj_authlat()</i>
for the geodetic ==> authalic), which is used for ISEA, but
also by:</p>
<p>- Equal Area Cylindrical (cea.cpp)<br>
- EqualEarth (eqearth.cpp)<br>
- HEALPix (healpix.cpp) -- which might also imply rHEALPix<br>
- Lambert Azimuthal Equal Area (laea.cpp)</p>
<p>The current version of <i>pj_authlat()</i> uses a series
expansion based on the eccentricity squared (I believe this is
based on a <a
href="https://en.wikipedia.org/wiki/Lagrange_inversion_theorem">Lagrange
reversion</a> as <a
href="https://en.wikipedia.org/wiki/Latitude#Inverse_formulae_and_series">mentioned
here</a>), and Charles Karney in <a
href="https://arxiv.org/pdf/2212.05818">this paper</a>
establishes that:<br>
</p>
<p>
<blockquote type="cite"><i>series expansions using the third
flattening as the small parameter uniformly result in smaller
truncation errors compared to using the eccentricity.</i></blockquote>
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.<br>
The current order 3, which requires 2 <i>sin()</i> calls, gives
precision up to roughly 8 decimals.<br>
Order 6 would mean twice the number of <i>sin()</i> calls, and
the 8 decimals should be enough to <a
href="https://xkcd.com/2170/">point to Waldo on a page</a>,
however C.K. argues that:</p>
<p>
<blockquote type="cite"><i>Modern geodetic libraries should strive
to achieve full double-precision accuracy </i>(page 4)<i><br>
</i></blockquote>
</p>
<p>and I would much agree with that.<br>
</p>
<p>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 <a
href="https://github.com/OSGeo/PROJ/pull/4211#pullrequestreview-2233903246">discussed
here</a>), and adjustments for tests etc.</p>
<p>PROJ currently does not use the full GeographicLib which already
includes code for conversion between different auxiliary
latitudes, does it?</p>
<p> (
<a class="moz-txt-link-freetext" href="https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp">https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp</a>
)<br>
</p>
<p>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 <i>pj_authlat()</i>.</p>
<p>Kind regards,</p>
<p>-Jerome<br>
</p>
<blockquote type="cite"><font size="5" face="monospace">void
pj_authset(double a, double b, double cp[6]) <br>
{ <br>
// <a class="moz-txt-link-freetext"
href="https://arxiv.org/pdf/2212.05818">https://arxiv.org/pdf/2212.05818</a>
<br>
// ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1) -- (20) <br>
double n = (a - b) / (a + b); // Third flattening <br>
static const double C[21] = // (A20) <br>
{ <br>
4 / 3.0, 4 / 45.0, -16/35.0, -2582 /14175.0, 60136
/467775.0, 28112932/ 212837625.0, <br>
46 / 45.0, 152/945.0, -11966 /14175.0, -21016 /
51975.0, 251310128/ 638512875.0, <br>
3044/2835.0, 3802 /14175.0, -94388 /
66825.0, -8797648/ 10945935.0, <br>
6059 / 4725.0, 41072 /
93555.0, -1472637812/ 638512875.0, <br>
768272
/467775.0, -455935736/ 638512875.0, <br>
4210684958/1915538625.0 <br>
}; <br>
double p; <br>
<br>
cp[0] = C[ 0] * n; <br>
<br>
p = n * n; <br>
cp[0] += C[ 1] * p; <br>
cp[1] = C[ 6] * p; <br>
<br>
p *= n; <br>
cp[0] += C[ 2] * p; <br>
cp[1] += C[ 7] * p; <br>
cp[2] = C[11] * p; <br>
<br>
p *= n; <br>
cp[0] += C[ 3] * p; <br>
cp[1] += C[ 8] * p; <br>
cp[2] += C[12] * p; <br>
cp[3] = C[15] * p; <br>
<br>
p *= n; <br>
cp[0] += C[ 4] * p; <br>
cp[1] += C[ 9] * p; <br>
cp[2] += C[13] * p; <br>
cp[3] += C[16] * p; <br>
cp[4] = C[18] * p; <br>
<br>
p *= n; <br>
cp[0] += C[ 5] * p; <br>
cp[1] += C[10] * p; <br>
cp[2] += C[14] * p; <br>
cp[3] += C[17] * p; <br>
cp[4] += C[19] * p; <br>
cp[5] = C[20] * p; <br>
} <br>
<br>
double pj_authlat(double auth, const double * cp) <br>
{ <br>
// <a class="moz-txt-link-freetext"
href="https://arxiv.org/pdf/2212.05818">https://arxiv.org/pdf/2212.05818</a>
<br>
// ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1) -- (20) <br>
double a2x = auth + auth, a4x = a2x + a2x, a8x = a4x + a4x; <br>
return auth + <br>
cp[0] * sin(a2x) + <br>
cp[1] * sin(a4x) + <br>
cp[2] * sin(a4x + a2x) + <br>
cp[3] * sin(a8x) + <br>
cp[4] * sin(a8x + a2x) + <br>
cp[5] * sin(a8x + a4x); <br>
} </font><br>
</blockquote>
</body>
</html>