Implement uniform conversions between auxiliary latitudes #4446
Conversation
This follows on the discussion in OSGeo#4441. In particular... latitudes.cpp now offers methods for converting between auxiliary laitudes using trigonometric series whose coefficients are power series in n. The series are evaluated using Clenshaw summation and Horner's method. For small flattening, these often offer better accuracy than the analytical expressions. Now mlfn.cpp and tmerc.cpp use these conversions, simplifying the implementations in these files. In particular tmerc's functions gatg and clens (which puzzlingly seemed to do exactly the same thing!) have been removed. tmerc's clenS (with a capital "S") stays -- it does complex Clenshaw summation Along the way, I've simplified pj_conformal_lat and pj_conform_lat_inverse in latitudes.cpp. For the former, I've substituted a better conditioned formula; for the latter, I used the existing (and more efficent) approach implemented in pj_sinhpsi2tanphi. The computation of authalic latitudes in latitudes.cpp now uses the series method if the flattening is small. However some of the test cases of the spillhaus projection use large flattening (is this really necessary?); so I've had to retain the poorly conditioned analytical formulas as well. If we really wanted to support large flattening in this case, it's necessary to evaluate the authalic latitude using Eqs. (52), (53), (54) of "On auxiliary latitudes" (https://arxiv.org/abs/2212.05818) and using Newton's method on the tangents of the angles for the inverse, using Eq. (57). I don't see the need for these at present. I've simplified the interface somewhat: passing P instead of P->e, P->onees, etc. separately; dropping "_coeffs" off the q function name and the tagging of the inverse function name with "_exact". I had hoped to retire the external use of Snyder's q(sinphi). However, some equal area projections, e.g., the cylindrical version, really do need this function instead of the authalic latitude directly.
|
There is no hurry on this -- it doesn't fix any bugs. However, it removes code duplication and makes some improvements in accuracy. |
Are we changing the flattening in the spilhaus tests? I am not aware of that. We use either the normal ellipsoid or a sphere. |
I misspoke. It was healpix and not spillhaus that used large flattening in its tests. In any case, I changed no tests, and they all still pass. |
Also auxlat.mac is now checked in on devel branch of git@github.com:geographiclib/geographiclib.git
| // Instead of calling tan and sin, call sin and cos which the compiler | ||
| // optimizes to a single call to sincos. | ||
| double sphi = sin(phi), cphi = cos(phi); | ||
| return atan(sinh(asinh(sphi / cphi) - P->e * atanh(P->e * sphi))); |
There was a problem hiding this comment.
Isn't if throwing for phi = 90 deg?
In Mercator projection it is not defined. But in spilhaus it is (near the center of the projection). Maybe se should add a test case in spilhaus for latitudes -90 and 90 if it is not already there.
There was a problem hiding this comment.
There's no error. cos(pi/2) evaluates in double precision is a small positive quantity, so the returned value will be very close to (or exactly) pi/2 (rounded to a double). Even if cos(pi/2) is zero, we would have sinh(asinh(inf)) -> inf, and atan(inf) = pi/2.
The only issue is if the precision is changed. With long double (64 bit precision), cos(pi/2) is negative and then we would get -pi/2 for the result. I suppose we could guard against this with fabs(cphi); but so many other things would need to be tweaked in PROJ, I don't think this is warranted. N.B. with quad precision (113 bits), cos(pi/2) is positve once again.
mwtoews
changed the title
Remove possibly confusing aliases in AuxLat (CHI for CONFORMAL etc.) Remove the stutter in the naming: AuxLat::AUXNUMBER -> AuxLat::NUMBER, AuxLat::AUXORDER -> AuxLat::ORDER. Empasize the NUMBER and ORDER are distinct. Align documentation for pj_auxlat.* routines with paper on auxiliary latitudes. Detail how to regenerate the array of coefficients. In pj_auxlat_convert routines, make the order, K, an optional argument (defaulting to AuxLat::ORDER). In tmerc, give a guide to the terminology used bu Engsager and Poder.
4 participants
This follows on the discussion in #4441. In particular...
latitudes.cpp now offers methods for converting between auxiliary laitudes using trigonometric series whose coefficients are power series in n. The series are evaluated using Clenshaw summation and Horner's method. For small flattening, these often offer better accuracy than the analytical expressions.
Now mlfn.cpp and tmerc.cpp use these conversions, simplifying the implementations in these files. In particular tmerc's functions gatg and clens (which puzzlingly seemed to do exactly the same thing!) have been removed. tmerc's clenS (with a capital "S") stays -- it does complex Clenshaw summation
Along the way, I've simplified pj_conformal_lat and pj_conform_lat_inverse in latitudes.cpp. For the former, I've substituted a better conditioned formula; for the latter, I used the existing (and more efficent) approach implemented in pj_sinhpsi2tanphi.
The computation of authalic latitudes in latitudes.cpp now uses the series method if the flattening is small. However some of the test cases of the healpix projection use large flattening (is this really necessary?); so I've had to retain the poorly conditioned analytical formulas as well.
If we really wanted to support large flattening in this case, it's necessary to evaluate the authalic latitude using Eqs. (52), (53), (54) of "On auxiliary latitudes" (https://arxiv.org/abs/2212.05818) and using Newton's method on the tangents of the angles for the inverse, using Eq. (57). I don't see the need for these at present.
I've simplified the interface somewhat: passing P instead of P->e, P->onees, etc. separately; dropping "_coeffs" off the q function name and the tagging of the inverse function name with "_exact". I had hoped to retire the external use of Snyder's q(sinphi). However, some equal area projections, e.g., the cylindrical version, really do need this function instead of the authalic latitude directly.
docs/source/*.rstfor new API