move conformal lat funcs to latitudes.cpp and few details

src/latitudes.cpp

@@ -0,0 +1,69 @@
+/*
+ * Project:  PROJ
+ *
+ * Helper functions to compute latitudes
+ *
+ * Some from Map Projections - A Working Manual. 1987. John P. Snyder
+ * https://neacsu.net/docs/geodesy/snyder/2-general/sect_3/
+ *
+ * Copyright (c) 2025 Javier Jimenez Shaw
+ */
+
+#include <math.h>
+
+#include "proj.h"
+#include "proj_internal.h"
+
+/*****************************************************************************/
+double pj_conformal_lat(double phi, double e) {
+    /***********************************
+     * The conformal latitude chi
+     * giving a sphere which is truly conformal in accordance with the ellipsoid
+     * Snyder, A working manual, formula 3-1
+     *
+     * phi: geodetic latitude, in radians
+     * e: ellipsoid eccentricity
+     * returns: conformal latitude, in radians
+     ***********************************/
+    if (e == 0.0)
+        return phi;
+
+    const double esinphi = e * sin(phi);
+    const double chi =
+        2 * atan(tan(M_FORTPI + phi / 2) *
+                 std::pow(((1 - esinphi) / (1 + esinphi)), (e / 2))) -
+        M_HALFPI;
+    return chi;
+}
+
+/*****************************************************************************/
+double pj_conformal_lat_inverse(double chi, double e, double threshold) {
+    /***********************************
+     * The inverse formula of the conformal latitude
+     * for phi (geodetic latitude) in terms of chi (conformal latitude),
+     * Snyder, A working manual, formula 3-4
+     *
+     * chi: conformal latitude, in radians
+     * e: ellipsoid eccentricity
+     * threshold: between two consecutive iteratinons to break the loop, radians
+     * returns: geodetic latitude, in radians
+     ***********************************/
+    if (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;
+}

src/lib_proj.cmake

@@ -227,6 +227,7 @@ set(SRC_LIBPROJ_CORE
   initcache.cpp
   internal.cpp
   inv.cpp
+  latitudes.cpp
   list.cpp
   log.cpp
   malloc.cpp

src/proj_internal.h

@@ -922,6 +922,9 @@ double pj_sinhpsi2tanphi(PJ_CONTEXT *, const double, const double);
 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);
+
 COMPLEX pj_zpoly1(COMPLEX, const COMPLEX *, int);
 COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *);
 

src/projections/spilhaus.cpp

@@ -20,7 +20,6 @@
 
 #include "proj.h"
 #include "proj_internal.h"
-#include <math.h>
 
 C_NAMESPACE PJ *pj_adams_ws2(PJ *);
 
@@ -38,23 +37,14 @@ struct pj_spilhaus_data {
     PJ *adams_ws2;
 };
 
-double conformal_lat(double lat, double e) {
-    // Snyder, A working manual, formula 3-1
-    const double esinlat = e * sin(lat);
-    const double chi =
-        2 * atan(tan(M_FORTPI + lat / 2) *
-                 std::pow(((1 - esinlat) / (1 + esinlat)), (e / 2))) -
-        M_HALFPI;
-    return chi;
-}
 } // anonymous namespace
 
 static PJ_XY spilhaus_s_forward(PJ_LP lp, PJ *P) {
     PJ_XY xy = {0.0, 0.0};
     struct pj_spilhaus_data *Q =
         static_cast<struct pj_spilhaus_data *>(P->opaque);
 
-    const double phi_c = conformal_lat(lp.phi, P->e);
+    const double phi_c = pj_conformal_lat(lp.phi, P->e);
     const double cosphi_c = cos(phi_c);
     const double sinphi_c = sin(phi_c);
 
@@ -81,7 +71,7 @@ static PJ_XY spilhaus_s_forward(PJ_LP lp, PJ *P) {
     xy.x = -(xyadams.x * Q->cosrot + xyadams.y * Q->sinrot) * factor;
     xy.y = -(xyadams.x * -Q->sinrot + xyadams.y * Q->cosrot) * factor;
 
-    return (xy);
+    return xy;
 }
 
 static PJ_LP spilhaus_s_inverse(PJ_XY xy, PJ *P) {
@@ -108,17 +98,10 @@ static PJ_LP spilhaus_s_inverse(PJ_XY xy, PJ *P) {
              aatan2(cosphi_s * sinlam_s,
                     Q->sinalpha * cosphi_s * coslam_s - Q->cosalpha * sinphi_s);
 
-    // latitude (iterative formula from
-    // https://mathworld.wolfram.com/ConformalLatitude.html)
-    const double taninit = tan(M_PI / 4 + lp.phi / 2);
-    for (int i = 0; i < 10; i++) {
-        const double esinlat = P->e * sin(lp.phi);
-        lp.phi = 2 * atan(taninit * std::pow(((1 + esinlat) / (1 - esinlat)),
-                                             (P->e / 2))) -
-                 0.5 * M_PI;
-    }
-
-    return (lp);
+    const double THRESHOLD = 1e-10;
+    lp.phi = pj_conformal_lat_inverse(lp.phi, P->e, THRESHOLD);
+
+    return lp;
 }
 
 static PJ *spilhaus_destructor(PJ *P, int errlev) { /* Destructor */
@@ -161,7 +144,7 @@ PJ *PJ_PROJECTION(spilhaus) {
     Q->cosrot = cos(rotation);
     Q->sinrot = sin(rotation);
 
-    const double conformal_lat_center = conformal_lat(lat_center, P->e);
+    const double conformal_lat_center = pj_conformal_lat(lat_center, P->e);
 
     Q->sinalpha = -cos(conformal_lat_center) * cos(azimuth);
     Q->cosalpha = sqrt(1 - Q->sinalpha * Q->sinalpha);

test/gie/spilhaus.gie

@@ -3394,6 +3394,16 @@ expect   -10546875   16640625
 accept   117.75729262   34.28539857
 expect   5078125   16640625
 
+
+------------------------------------------------------------
+# Stable for default parameters
+------------------------------------------------------------
+operation  +proj=spilhaus +rot=45 +k_0=1 +lat_c=-49.56371678 +lon_c=66.94970198 +azi=40.17823482
+tolerance  1 mm
+------------------------------------------------------------
+accept   130.4   -16.2
+expect   3733410.0118   -9320.8573
+
 ------------------------------------------------------------
 # Sentitive to input parameters
 ------------------------------------------------------------
@@ -3433,4 +3443,14 @@ tolerance  1 mm
 accept   130.4   -16.2
 expect   3360069.0106   -8388.7716
 
+------------------------------------------------------------
+# Sphere
+------------------------------------------------------------
+operation  +proj=spilhaus +R=6378137
+tolerance  1 mm
+------------------------------------------------------------
+accept   130.4   -16.2
+expect   3737644.5177   -7049.7883
+
+
 </gie-strict>