<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Sorry I should correct a few confusing mistakes in my previous
e-mail:</p>
<p> - I was talking mainly about authalic latitude ==>
geodetic latitude conversion, not the other way around.<br>
<i>pj_authlat()</i> returns the geodetic latitude for a
given authalic latitude (the name of the function is quite
confusing -- I would propose to rename to something less ambiguous
if it is only used internally)<br>
<br>
- The function currently uses 3 <i>sin()</i> calls for order
3 (not 2, made a typo)<br>
<br>
- The ISEA projection does not currently make use of <i>pj_authlat()</i>
or <i>pj_qsfn() </i>for automatically converting authalic /
geodetic latitudes, but I am proposing to introduce this when
using non-spherical ellipsoid (just like LAEA of which it is a
modified form).</p>
<p>Best,</p>
<p>-Jerome</p>
<div class="moz-cite-prefix">On 9/10/24 9:19 PM, Jérôme St-Louis via
PROJ wrote:<br>
</div>
<blockquote type="cite"
cite="mid:9ed4372a-f6a3-4420-b32b-232c3c4638b0@ecere.com">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<p>Hi Even, all,</p>
<p>As suggested in <a class="moz-txt-link-freetext"
href="https://github.com/OSGeo/PROJ/pull/4211"
moz-do-not-send="true">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"
moz-do-not-send="true">Lagrange reversion</a> as <a
href="https://en.wikipedia.org/wiki/Latitude#Inverse_formulae_and_series"
moz-do-not-send="true">mentioned here</a>), and Charles Karney
in <a href="https://arxiv.org/pdf/2212.05818"
moz-do-not-send="true">this paper</a> establishes that:<br>
</p>
<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/" moz-do-not-send="true">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>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"
moz-do-not-send="true">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"
moz-do-not-send="true">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"
moz-do-not-send="true">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"
moz-do-not-send="true">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>
<br>
<fieldset class="moz-mime-attachment-header"></fieldset>
<pre class="moz-quote-pre" wrap="">_______________________________________________
PROJ mailing list
<a class="moz-txt-link-abbreviated" href="mailto:PROJ@lists.osgeo.org">PROJ@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="https://lists.osgeo.org/mailman/listinfo/proj">https://lists.osgeo.org/mailman/listinfo/proj</a>
</pre>
</blockquote>
</body>
</html>