Improved series for the meridian length and its inverse

src/mlfn.cpp

@@ -1,52 +1,87 @@
 #include <math.h>
-
-#include "proj.h"
 #include "proj_internal.h"
-#include "mlfn.hpp"
 
-/* meridional distance for ellipsoid and inverse
-**	8th degree - accurate to < 1e-5 meters when used in conjunction
-**		with typical major axis values.
-**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
+/* meridional distance for ellipsoid and inverse using 6th-order expansion in
+** the third flattening n.  This gives full double precision accuracy for |f|
+** <= 1/150.
 */
-#define C00 1.
-#define C02 .25
-#define C04 .046875
-#define C06 .01953125
-#define C08 .01068115234375
-#define C22 .75
-#define C44 .46875
-#define C46 .01302083333333333333
-#define C48 .00712076822916666666
-#define C66 .36458333333333333333
-#define C68 .00569661458333333333
-#define C88 .3076171875
-#define EN_SIZE 5
+
+#define Lmax 6
+
+// Evaluation sum(p[i] * x^i, i, 0, N) via Horner's method.  N.B. p is of
+// length N+1.
+static double polyval(double x, const double p[], int N) {
+    double y = N < 0 ? 0 : p[N];
+    for (; N > 0;) y = y * x + p[--N];
+    return y;
+}
+
+// Evaluate y = sum(c[k] * sin((2*k+2) * zeta), i, 0, K-1)
+static double clenshaw(double szeta, double czeta,
+                       const double c[], int K) {
+    // Approx operation count = (K + 5) mult and (2 * K + 2) add
+    double u0 = 0, u1 = 0,                         // accumulators for sum
+        X = 2 * (czeta - szeta) * (czeta + szeta); // 2 * cos(2*zeta)
+    for (; K > 0;) {
+        double t = X * u0 - u1 + c[--K];
+        u1 = u0; u0 = t;
+    }
+    return 2 * szeta * czeta * u0; // sin(2*zeta) * u0
+}
 
 double *pj_enfn(double es) {
-    double t, *en;
 
-    en = (double *) malloc(EN_SIZE * sizeof (double));
-    if (nullptr==en)
-        return nullptr;
+    // Expansion of (quarter meridian) / ((a+b)/2 * pi/2) as series in n^2; these
+    // coefficients are ( (2*k - 3)!! / (2*k)!! )^2 for k = 0..3
+    static const double coeff_rad[] = {1, 1.0/4, 1.0/64, 1.0/256};
+
+    // Coefficients to convert phi to mu, Eq. A5 in arXiv:2212.05818
+    // with 0 terms dropped
+    static const double coeff_mu_phi[] = {
+        -3.0/2, 9.0/16, -3.0/32,
+        15.0/16, -15.0/32, 135.0/2048,
+        -35.0/48, 105.0/256,
+        315.0/512, -189.0/512,
+        -693.0/1280,
+        1001.0/2048,
+    };
+    // Coefficients to convert mu to phi, Eq. A6 in arXiv:2212.05818
+    // with 0 terms dropped
+    static const double coeff_phi_mu[] = {
+        3.0/2, -27.0/32, 269.0/512,
+        21.0/16, -55.0/32, 6759.0/4096,
+        151.0/96, -417.0/128,
+        1097.0/512, -15543.0/2560,
+        8011.0/2560,
+        293393.0/61440,
+    };
 
-    en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08)));
-    en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08)));
-    en[2] = (t = es * es) * (C44 - es * (C46 + es * C48));
-    en[3] = (t *= es) * (C66 - es * C68);
-    en[4] = t * es * C88;
+    double n = 1 + sqrt(1 - es); n = es / (n * n);
+    double n2 = n * n, d = n, *en;
 
+    // 2*Lmax for the Fourier coeffs for each direction of conversion + 1 for
+    // overall multiplier.
+    en = (double *) malloc((2*Lmax + 1) * sizeof (double));
+    if (nullptr==en)
+        return nullptr;
+    en[0] = polyval(n2, coeff_rad, Lmax/2) / (1 + n);
+    for (int l = 0, o = 0; l < Lmax; ++l) {
+        int m = (Lmax - l - 1) / 2;
+        en[l + 1       ] = d * polyval(n2, coeff_mu_phi + o, m);
+        en[l + 1 + Lmax] = d * polyval(n2, coeff_phi_mu + o, m);
+        d *= n;
+        o += m + 1;
+    }
     return en;
 }
 
 double
 pj_mlfn(double phi, double sphi, double cphi, const double *en) {
-    return inline_pj_mlfn(phi, sphi, cphi, en);
+    return en[0] * (phi + clenshaw(sphi, cphi, en + 1, Lmax));
 }
 
 double
-pj_inv_mlfn(PJ_CONTEXT *ctx, double arg, double es, const double *en) {
-    double sinphi_ignored;
-    double cosphi_ignored;
-    return inline_pj_inv_mlfn(ctx, arg, es, en, &sinphi_ignored, &cosphi_ignored);
+pj_inv_mlfn(double mu, const double *en) {
+    mu /= en[0];
+    return mu + clenshaw(sin(mu), cos(mu), en + 1 + Lmax, Lmax);
 }

src/proj_internal.h

@@ -815,7 +815,7 @@ void     *free_params (PJ_CONTEXT *ctx, paralist *start, int errlev);
 
 double *pj_enfn(double);
 double  pj_mlfn(double, double, double, const double *);
-double  pj_inv_mlfn(PJ_CONTEXT *, double, double, const double *);
+double  pj_inv_mlfn(double, const double *);
 double  pj_qsfn(double, double, double);
 double  pj_tsfn(double, double, double);
 double  pj_msfn(double, double, double);

src/projections/aeqd.cpp

@@ -191,8 +191,7 @@ static PJ_LP e_guam_inv(PJ_XY xy, PJ *P) { /* Guam elliptical */
     for (i = 0; i < 3; ++i) {
         t = P->e * sin(lp.phi);
         t = sqrt(1. - t * t);
-        lp.phi = pj_inv_mlfn(P->ctx, Q->M1 + xy.y -
-            x2 * tan(lp.phi) * t, P->es, Q->en);
+        lp.phi = pj_inv_mlfn(Q->M1 + xy.y - x2 * tan(lp.phi) * t, Q->en);
     }
     lp.lam = xy.x * t / cos(lp.phi);
     return lp;
@@ -223,8 +222,7 @@ static PJ_LP aeqd_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse
         lp.lam = lon2 * DEG_TO_RAD;
         lp.lam -= P->lam0;
     } else { /* Polar */
-        lp.phi = pj_inv_mlfn(P->ctx, Q->mode == N_POLE ? Q->Mp - c : Q->Mp + c,
-            P->es, Q->en);
+        lp.phi = pj_inv_mlfn(Q->mode == N_POLE ? Q->Mp - c : Q->Mp + c, Q->en);
         lp.lam = atan2(xy.x, Q->mode == N_POLE ? -xy.y : xy.y);
     }
     return lp;
@@ -328,5 +326,3 @@ PJ *PROJECTION(aeqd) {
 
     return P;
 }
-
-

src/projections/bonne.cpp

@@ -83,7 +83,7 @@ static PJ_LP bonne_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, invers
 
     xy.y = Q->am1 - xy.y;
     rh = hypot(xy.x, xy.y);
-    lp.phi = pj_inv_mlfn(P->ctx, Q->am1 + Q->m1 - rh, P->es, Q->en);
+    lp.phi = pj_inv_mlfn(Q->am1 + Q->m1 - rh, Q->en);
     if ((s = fabs(lp.phi)) < M_HALFPI) {
         s = sin(lp.phi);
         lp.lam = rh * atan2(xy.x, xy.y) *

src/projections/cass.cpp

@@ -68,7 +68,7 @@ static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse
     PJ_LP lp = {0.0, 0.0};
     struct cass_data *Q = static_cast<struct cass_data*>(P->opaque);
 
-    const double phi1 = pj_inv_mlfn (P->ctx, Q->m0 + xy.y, P->es, Q->en);
+    const double phi1 = pj_inv_mlfn (Q->m0 + xy.y, Q->en);
     const double tanphi1 = tan (phi1);
     const double T1 = tanphi1*tanphi1;
     const double sinphi1 = sin (phi1);

src/projections/eqdc.cpp

@@ -51,7 +51,7 @@ static PJ_LP eqdc_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse
         }
         lp.phi = Q->c - Q->rho;
         if (Q->ellips)
-            lp.phi = pj_inv_mlfn(P->ctx, lp.phi, P->es, Q->en);
+            lp.phi = pj_inv_mlfn(lp.phi, Q->en);
         lp.lam = atan2(xy.x, xy.y) / Q->n;
     } else {
         lp.lam = 0.;

src/projections/gn_sinu.cpp

@@ -38,7 +38,7 @@ static PJ_LP gn_sinu_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inve
     PJ_LP lp = {0.0,0.0};
     double s;
 
-    lp.phi = pj_inv_mlfn(P->ctx, xy.y, P->es, static_cast<struct pj_opaque*>(P->opaque)->en);
+    lp.phi = pj_inv_mlfn(xy.y, static_cast<struct pj_opaque*>(P->opaque)->en);
     s = fabs(lp.phi);
     if (s < M_HALFPI) {
         s = sin(lp.phi);

src/projections/lcca.cpp

@@ -115,7 +115,7 @@ static PJ_LP lcca_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse
         proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
         return lp;
     }
-    lp.phi = pj_inv_mlfn(P->ctx, S + Q->M0, P->es, Q->en);
+    lp.phi = pj_inv_mlfn(S + Q->M0, Q->en);
 
     return lp;
 }

src/projections/tmerc.cpp

@@ -18,7 +18,6 @@
 #include "proj.h"
 #include "proj_internal.h"
 #include <math.h>
-#include "mlfn.hpp"
 
 PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell\n\tapprox";
 PROJ_HEAD(etmerc, "Extended Transverse Mercator") "\n\tCyl, Sph";
@@ -106,7 +105,7 @@ static PJ_XY approx_e_fwd (PJ_LP lp, PJ *P)
         FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t)
         + FC7 * als * (61. + t * ( t * (179. - t) - 479. ) )
         )));
-    xy.y = P->k0 * (inline_pj_mlfn(lp.phi, sinphi, cosphi, Q->en) - Q->ml0 +
+    xy.y = P->k0 * (pj_mlfn(lp.phi, sinphi, cosphi, Q->en) - Q->ml0 +
         sinphi * al * lp.lam * FC2 * ( 1. +
         FC4 * als * (5. - t + n * (9. + 4. * n) +
         FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t)
@@ -155,12 +154,12 @@ static PJ_LP approx_e_inv (PJ_XY xy, PJ *P) {
     PJ_LP lp = {0.0,0.0};
     const auto *Q = &(static_cast<struct tmerc_data*>(P->opaque)->approx);
 
-    double sinphi, cosphi;
-    lp.phi = inline_pj_inv_mlfn(P->ctx, Q->ml0 + xy.y / P->k0, P->es, Q->en, &sinphi, &cosphi);
+    lp.phi = pj_inv_mlfn(Q->ml0 + xy.y / P->k0, Q->en);
     if (fabs(lp.phi) >= M_HALFPI) {
         lp.phi = xy.y < 0. ? -M_HALFPI : M_HALFPI;
         lp.lam = 0.;
     } else {
+        double sinphi = sin(lp.phi), cosphi = cos(lp.phi);
         double t = fabs (cosphi) > 1e-10 ? sinphi/cosphi : 0.;
         const double n = Q->esp * cosphi * cosphi;
         double con = 1. - P->es * sinphi * sinphi;

test/gie/builtins.gie

@@ -6926,7 +6926,7 @@ expect  687071.43910944  6210141.32674801    0.00000000     2000.0000
 operation +proj=utm +zone=32 +approx
 tolerance 0.001 mm
 accept	12 56 0 2000
-expect  687071.43911000  6210141.32675053    0.00000000     2000.0000
+expect  687071.43910944  6210141.32674801    0.00000000     2000.0000
 
 
 -------------------------------------------------------------------------------