[PROJ] Geodetic to Authalic latitude conversions

[Charles Karney]

On 9/12/24 9:18 AM, Even Rouault wrote:
> Jérôme,
> 
> feel free to make the code refactorings that make sense. All this 
> geodetic <--> authalic conversion is PROJ internals.
> 
> Regarding pj_auth2geodlat() signature, I would suggest that it is:
> 
> double pj_auth2geodlat(const double * cp, double sinphi, double cosphi)
> 
> since looking at the code, at least in cea, eqearth and laea, we 
> actually compute phi by doing a "asin(something)" before
> 
> So instead of doing:
> 
> phi = asin(something);
> phi = pj_auth2geodlat(cp, phi);
> 
> we could just do:
> 
> phi = pj_auth2geodlat(cp, something, cos(asin(something)))
> 
> And that would save us one sin() call
> 
> And potentially cos(asin(something)) = sqrt(1 - something^2)

Yes, this is a good idea, not just to save on a trig evaluation, but
because specifying cosphi lets you resolve angles very close to 90d.
However, the bird has already escaped the cage, if all you know is
sinphi -- you've lost 1/2 of the precision near the pole.

So, a task for another day, is to go back to the formulas for sinphi and
manipulate them to give accurate results for cosphi.

This also suggests that pj_auth2geodlat might instead profitably return
the sine and cosine of the authalic latitude.  But this, also, is
perhaps a task for another day.

On cos(asin(x)) vs sqrt(1 - x^2): a third and better(?) method of
calculating this is sqrt((1 - x) * (1 + x)), because, for x >= 0.5,
(1 - x) is computed no roundoff error.  For x = 1 - 1e-7, the relative
errors of the three methods are:

   2.02e-11, 2.00e-11, 5e-19