[PROJ] spilhaus regresssion test failures

[Javier Jimenez Shaw]

Sorry, I should me more explicit
What I mean is:
 - take 9.6.0
 - replace the function pj_conformal_lat_inverse in latitudes.cpp with the
new implementation
 - replace its declaration in proj_internal.h
 - replace its usage in spilhaus.cpp

... that is this patch

------- begin patch --------------
diff --git a/src/latitudes.cpp b/src/latitudes.cpp
index 7d49f667..c5e277f8 100644
--- a/src/latitudes.cpp
+++ b/src/latitudes.cpp
@@ -36,33 +36,17 @@ double pj_conformal_lat(double phi, double e) {
 }

 /*****************************************************************************/
-double pj_conformal_lat_inverse(double chi, double e, double threshold) {
+
+double pj_conformal_lat_inverse(double chi, const PJ *P) {
     /***********************************
-     * The inverse formula of the conformal latitude
-     * for phi (geodetic latitude) in terms of chi (conformal latitude),
-     * Snyder, A working manual, formula 3-4
+     * The geodetic latitude, phi, in terms of the conformal latitude, chi.
      *
      * chi: conformal latitude, in radians
-     * e: ellipsoid eccentricity
-     * threshold: between two consecutive iteratinons to break the loop,
radians
-     * returns: geodetic latitude, in radians
+     * returns: geodetic latitude, phi, in radians
+     * Copied from merc.cpp
      ***********************************/
-    if (e == 0.0)
+    if (P->e == 0.0)
         return chi;

-    const double taninit = tan(M_PI / 4 + chi / 2);
-    double phi = chi;
-    for (int i = 0; i < 10; i++) {
-        const double esinphi = e * sin(phi);
-        const double new_phi =
-            2 * atan(taninit *
-                     std::pow(((1 + esinphi) / (1 - esinphi)), (e / 2))) -
-            0.5 * M_PI;
-        if (abs(new_phi - phi) < threshold) {
-            phi = new_phi;
-            break;
-        }
-        phi = new_phi;
-    }
-    return phi;
+    return atan(pj_sinhpsi2tanphi(P->ctx, tan(chi), P->e));
 }
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 922c7074..54b64341 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -923,7 +923,7 @@ double *pj_authset(double);
 double pj_authlat(double, double *);

 double pj_conformal_lat(double phi, double e);
-double pj_conformal_lat_inverse(double chi, double e, double threshold);
+double pj_conformal_lat_inverse(double chi, const PJ *P);

 COMPLEX pj_zpoly1(COMPLEX, const COMPLEX *, int);
 COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *);
diff --git a/src/projections/spilhaus.cpp b/src/projections/spilhaus.cpp
index 01b86800..ba9a6c57 100644
--- a/src/projections/spilhaus.cpp
+++ b/src/projections/spilhaus.cpp
@@ -98,8 +98,7 @@ static PJ_LP spilhaus_inverse(PJ_XY xy, PJ *P) {
              aatan2(cosphi_s * sinlam_s,
                     Q->sinalpha * cosphi_s * coslam_s - Q->cosalpha *
sinphi_s);

-    const double THRESHOLD = 1e-10;
-    lp.phi = pj_conformal_lat_inverse(lp.phi, P->e, THRESHOLD);
+    lp.phi = pj_conformal_lat_inverse(lp.phi, P);

     return lp;
 }
------- end patch --------------


On Sun, 4 May 2025 at 20:27, Greg Troxel <gdt at lexort.com> wrote:

> Javier Jimenez Shaw <j1 at jimenezshaw.com> writes:
>
> > The PR was https://github.com/OSGeo/PROJ/pull/4446
> >
> > Maybe the function that impacts the tests is "pj_conformal_lat_inverse".
> > You can try with the new implementation.
> > The changes should be small. (See that the signature of the function is a
> > bit different). It is used only in one place.
>
>
> Do you mean:
>
>   take proj 9.6.0
>
>   read the diff to pj_conformal_lat_inverse from the PR, and hand patch
>   it in, except not changing the signature
>
>   build and rerun tests
>
> I guess I could further keep the code and compare the two results.
>
>
> I looked at that code, and it jumped out at me that are abs() was used,
> which is a C function, not C++, and is integer only.
>
> With this patch (to 9.6.0), I get a pass:
>
> --- src/latitudes.cpp.~1~       2025-03-13 10:31:27.000000000 -0400
> +++ src/latitudes.cpp   2025-05-04 14:23:04.656697860 -0400
> @@ -58,7 +58,7 @@
>              2 * atan(taninit *
>                       std::pow(((1 + esinphi) / (1 - esinphi)), (e / 2))) -
>              0.5 * M_PI;
> -        if (abs(new_phi - phi) < threshold) {
> +        if (std::abs(new_phi - phi) < threshold) {
>              phi = new_phi;
>              break;
>          }
>
> Interesting that this only causes trouble with spilhaus, at least as far
> as tests catch it.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20250504/334932df/attachment.htm>