From f2eba2e008414010eb7ce0d7a8993d5f4e01754f Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Fri, 16 Dec 2022 18:27:42 -0500 Subject: [PATCH] Improved series for the meridian length and its inverse Replace the series for the meridian length in the functions pj_enfn pj_mlfn pj_inv_mlfn in mlfn.cpp with improved versions. The old version of pj_mlfn used Snyder's Eq. (3-21) extended from 3rd to 4th order in e^2 and pj_inv_mlfn inverted this by an iterative method. The new version uses 6th order expansions in the third flattening, n, for both pj_mlfn and pj_inv_mlfn. Expansions in n have a smaller truncation error that the corresponding expansion in e^2. Increasing the order to 6th (from 4th) means that the routines give full double precision accuracy for |f| <= 1/150. Because the meridian length is used in several projections, it's important that it be accurate as possible. For example, now the azimuthal equidistant projection is equally accurate in the polar (uses pj_mlfn) and oblique (uses the geodesic routines) aspects. The increased accuracy tripped one of the tests: an instance of the projection +proj=utm +approx. The fix is to use the non-approx results for the expected values. Presumably this change results in pj_mlfn being a little slower while pj_inv_mlfn is considerably faster (since there's no iteration involved). Incidentally the new method for pj_inv_mlfn uses Snyder's recommended series Eq. (3-23) extended from 4th to 6th order (Snyder writes n as e_1). More discussion of the series used here is given in my recent paper On auxiliary latitudes, https://arxiv.org/abs/2212.05818 --- src/mlfn.cpp | 103 ++++++++++++++++++++++++------------ src/mlfn.hpp | 73 ------------------------- src/proj_internal.h | 2 +- src/projections/aeqd.cpp | 8 +-- src/projections/bonne.cpp | 2 +- src/projections/cass.cpp | 2 +- src/projections/eqdc.cpp | 2 +- src/projections/gn_sinu.cpp | 2 +- src/projections/lcca.cpp | 2 +- src/projections/tmerc.cpp | 7 ++- test/gie/builtins.gie | 2 +- 11 files changed, 81 insertions(+), 124 deletions(-) delete mode 100644 src/mlfn.hpp diff --git a/src/mlfn.cpp b/src/mlfn.cpp index 763d83ee78..32023fe1cc 100644 --- a/src/mlfn.cpp +++ b/src/mlfn.cpp @@ -1,52 +1,87 @@ #include - -#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); } diff --git a/src/mlfn.hpp b/src/mlfn.hpp deleted file mode 100644 index 426dcb0ed2..0000000000 --- a/src/mlfn.hpp +++ /dev/null @@ -1,73 +0,0 @@ -#ifndef MLFN_HPP -#define 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. -*/ - -inline static double inline_pj_mlfn(double phi, double sphi, double cphi, const double *en) { - cphi *= sphi; - sphi *= sphi; - return(en[0] * phi - cphi * (en[1] + sphi*(en[2] - + sphi*(en[3] + sphi*en[4])))); -} - -inline static double -inline_pj_inv_mlfn(PJ_CONTEXT *ctx, double arg, double es, const double *en, - double* sinphi, double* cosphi) { - const double k = 1./(1.-es); - constexpr double INV_MLFN_EPS = 1e-11; - constexpr int INV_MFN_MAX_ITER = 10; - double phi = arg; - double s = sin(phi); - double c = cos(phi); - for (int i = INV_MFN_MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */ - double t = 1. - es * s * s; - t = (inline_pj_mlfn(phi, s, c, en) - arg) * (t * sqrt(t)) * k; - phi -= t; - if (fabs(t) < INV_MLFN_EPS) - { - // Instead of recomputing sin(phi), cos(phi) from scratch, - // use sin(phi-t) and cos(phi-t) approximate formulas with - // 1-term approximation of sin(t) and cos(t) - *sinphi = s - c * t; - *cosphi = c + s * t; - return phi; - } - if (fabs(t) < 1e-3) - { - // 2-term approximation of sin(t) and cos(t) - // Max relative error is 4e-14 on cos(t), and 8e-15 on sin(t) - const double t2 = t * t; - const double cos_t = 1 - 0.5 * t2; - const double sin_t = t * (1 - (1. / 6) * t2); - const double s_new = s * cos_t - c * sin_t; - c = c * cos_t + s * sin_t; - s = s_new; - } - else if (fabs(t) < 1e-2) - { - // 3-term approximation of sin(t) and cos(t) - // Max relative error is 2e-15 on cos(t), and 2e-16 on sin(t) - const double t2 = t * t; - const double cos_t = 1 - 0.5 * t2 * (1 - (1. / 12) * t2); - const double sin_t = t * (1 - (1. / 6) * t2 * (1 - (1. / 20) * t2)); - const double s_new = s * cos_t - c * sin_t; - c = c * cos_t + s * sin_t; - s = s_new; - } - else - { - s = sin(phi); - c = cos(phi); - } - } - *sinphi = s; - *cosphi = c; - proj_context_errno_set( ctx, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN ); - return phi; -} - -#endif // MLFN_HPP diff --git a/src/proj_internal.h b/src/proj_internal.h index bf5bb797da..a4c6315f09 100644 --- a/src/proj_internal.h +++ b/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); diff --git a/src/projections/aeqd.cpp b/src/projections/aeqd.cpp index f01edff940..a601785045 100644 --- a/src/projections/aeqd.cpp +++ b/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; } - - diff --git a/src/projections/bonne.cpp b/src/projections/bonne.cpp index 99a0cc8f2e..3a779bd711 100644 --- a/src/projections/bonne.cpp +++ b/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) * diff --git a/src/projections/cass.cpp b/src/projections/cass.cpp index 40baa8f0ca..3a76f128a7 100644 --- a/src/projections/cass.cpp +++ b/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(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); diff --git a/src/projections/eqdc.cpp b/src/projections/eqdc.cpp index e408b7774b..85f40bcc52 100644 --- a/src/projections/eqdc.cpp +++ b/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.; diff --git a/src/projections/gn_sinu.cpp b/src/projections/gn_sinu.cpp index 7115b42d02..a2f5622959 100644 --- a/src/projections/gn_sinu.cpp +++ b/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(P->opaque)->en); + lp.phi = pj_inv_mlfn(xy.y, static_cast(P->opaque)->en); s = fabs(lp.phi); if (s < M_HALFPI) { s = sin(lp.phi); diff --git a/src/projections/lcca.cpp b/src/projections/lcca.cpp index 7b5eec7a37..18f1dbe015 100644 --- a/src/projections/lcca.cpp +++ b/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; } diff --git a/src/projections/tmerc.cpp b/src/projections/tmerc.cpp index 3385ffaac3..6ad829d002 100644 --- a/src/projections/tmerc.cpp +++ b/src/projections/tmerc.cpp @@ -18,7 +18,6 @@ #include "proj.h" #include "proj_internal.h" #include -#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(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; diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 567abea0f3..61f236f42a 100644 --- a/test/gie/builtins.gie +++ b/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 -------------------------------------------------------------------------------