From 5d22d137b1be282bdc14c1d12ab8f669f58d41a6 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Mon, 26 Oct 2020 12:44:18 -0400 Subject: Update Mercator projection Introduction ------------ The existing formulation for the Mercator projection is "satisfactory"; it is reasonably accurate. However for a core projection like Mercator, I think we should strive for full double precision accuracy. This commit uses cleaner, more accurate, and faster methods for computing the forward and inverse projections. These use the formulation in terms of hyperbolic functions that are manifestly odd in latitude psi = asinh(tan(phi)) - e * atanh(e * sin(phi)) (phi = latitude; psi = isometric latitude = Mercator y coordinate). Contrast this with the existing formulation psi = log(tan(pi/4 - phi/2)) - e/2 * log((1 + e * sin(phi)) / (1 - e * sin(phi))) where psi(-phi) isn't exactly equal to -psi(phi) and psi(0) isn't guaranteed to be 0. Implementation -------------- There's no particular issue implementing the forward projection, just apply the formulas above. The inverse projection is tricky because there's no closed form solution for the inverse. The existing code for the inverse uses an iterative method from Snyder. This is the usual hokey function iteration, and, as usual, the convergence rate is linear (error reduced by a constant factor on each iteration). This is OK (just) for low accuracy work. But nowadays, something with quadratic convergence (e.g., Newton's method, number of correct digits doubles on each iteration) is preferred (and used here). More on this later. The solution for phi(psi) I use is described in my TM paper and I lifted the specific formulation from GeographicLib's Math::tauf, which uses the same underlying machinery for all conformal projections. It solves for tan(phi) in terms of sinh(psi) which as a near identity mapping is ideal for Newton's method. For comparison I also look at the approach adopted by Poder + Engsager in their TM paper and implemented in etmerc. This uses trigonometric series (accurate to n^6) to convert phi <-> chi. psi is then given by psi = asinh(tan(chi)) Accuracy -------- I tested just the routines for transforming phi <-> psi from merc.cpp and measured the errors (converted to true nm = nanometers) for the forward and inverse mapping. I also included in my analysis the method used by etmerc. This uses a trigonometric series to convert phi <-> chi = atan(sinh(psi)), the conformal latitude. forward inverse max rms max rms old merc 3.60 0.85 2189.47 264.81 etmerc 1.82 0.38 1.42 0.37 new merc 1.83 0.30 2.12 0.31 1 nm is pretty much the absolute limit for accuracy in double precision (1 nm = 10e6 m / 2^53, approximately), and 5 nm is probably the limit on what you should routinely expect. So the old merc inverse is considerably less accurate that it could be. The old merc forward is OK on accuracy -- except that if does not preserve the parity of the projection. The accuracy of etmerc is fine (the truncation error of the 6th order series is small compared with the round-off error). However, situation reverses as the flattening is increased. E.g., at f = 1/150, the max error for the inverse projection is 8 nm. etmerc is OK for terrestrial applications, but couldn't be used for Mars. Timing ------ Here's what I get with g++ -O3 on various Linux machines with recent versions of g++. As always, you should take these with a grain of salt. You might expect the relative timings to vary by 20% or so when switching between compilers/machines. Times per call in ns = nanoseconds. forward inverse old merc 121 360 etmerc 4e-6 1.4 new merc 20 346 The new merc method is 6 times faster at the forward projection and modestly faster at the inverse projection (despite being more accurate). The latter result is because it only take 2 iterations of Newton's method to get full accuracy compared with an average of 5 iterations for the old method to get only um accuracy. A shocking aspect of these timings is how fast etmerc is. Another is that forward etmerc is streaks faster that inverse etmerc (it made be doubt my timing code). Evidently, asinh(tan(chi)) is a lot faster to compute than atan(sinh(psi)). The hesitation about adopting etmerc then comes down to: * the likelihood that Mercator may be used for non-terrestrial bodies; * the question of whether the timing benefits for the etmerc method would be noticeable in a realistic application; * need to duplicate the machinery for evaluating the coefficients for the series and for Clenshaw summation in the current code layout. Ripple effects ============== The Mercator routines used the the Snyder method, pj_tsfn and pj_phi2, are used in other projections. These relate phi to t = exp(-psi) (a rather bizarre choice in my book). I've retrofitted these to use the more accurate methods. These do the "right thing" for phi in [-pi/2, pi/2] , t in [0, inf], and e in [0, 1). NANs are properly handled. Of course, phi = pi/2 in double precision is actually less than pi/2, so cos(pi/2) > 0. So no special handling is needed for pi/2. Even if angles were handled in such a way that 90deg were exactly represented, these routines would still "work", with, e.g., tan(pi/2) -> inf. (A caution: with long doubles = a 64-bit fraction, we have cos(pi/2) < 0; and now we would need to be careful.) As a consequence, there no need for error handling in pj_tsfn; the HUGE_VAL return has gone and, of course, HUGE_VAL is a perfectly legal input to tsfn's inverse, phi2, which would return -pi/2. This "error handling" was only needed for e = 1, a case which is filtered out upstream. I will note that bad argument handling is much more natural using NAN instead of HUGE_VAL. See issue #2376 I've renamed the error condition for non-convergence of the inverse projection from "non-convergent inverse phi2" to "non-convergent sinh(psi) to tan(phi)". Now that pj_tsfn and pj_phi2 now return "better" results, there were some malfunctions in the projections that called them, specifically gstmerc, lcc, and tobmerc. * gstmerc invoked pj_tsfn(phi, sinphi, e) with a value of sinphi that wasn't equal to sin(phi). Disaster followed. I fixed this. I also replaced numerous occurrences of "-1.0 * x" by "-x". (Defining a function with arguments phi and sinphi is asking for trouble.) * lcc incorrectly thinks that the projection isn't defined for standard latitude = +/- 90d. This happens to be false (it reduces to polar stereographic in this limit). The check was whether tsfn(phi) = 0 (which only tested for the north pole not the south pole). However since tsfn(pi/2) now (correctly) returns a nonzero result, this test fails. I now just test for |phi| = pi/2. This is clearer and catches both poles (I'm assuming that the current implementation will probably fail in these cases). * tobmerc similarly thinks that phi close to +/- pi/2 can't be transformed even though psi(pi/2) is only 38. I'm disincline to fight this. However I did tighten up the failure condition (strict equality of |phi| == pi/2). OTHER STUFF =========== Testing ------- builtins.gei: I tightened up the tests for merc (and while I was about it etmerc and tmerc) to reflect full double precision accuracy. My test values are generated with MPFR enabled code and so should be accurate to all digits given. For the record, for GRS80 I use f = 1/298.2572221008827112431628366 in these calculations. pj_phi2_test: many of the tests were bogus testing irrelevant input parameters, like negative values of exp(-psi), and freezing in the arbitrary behavior of phi2. I've reworked most for the tests to be semi-useful. @schwehr can you review. Documentation ------------- I've updated merc.rst to outline the calculation of the inverse projection. phi2.cpp includes detailed notes about applying Newton's method to find tan(phi) in terms of sinh(psi). Future work ----------- lcc needs some tender loving care. It can easily (and should) be modified to allow stdlat = +/- 90 (reduces to polar stereographic), stdlat = 0 and stdlat_1 + stdlat_2 = 0 (reduces to Mercator). A little more elbow grease will allow the treatment of stdlat_1 close to stdlat_2 using divided differences. (See my implementation of the LambertConformalConic class in GeographicLib.) All the places where pj_tsfn and pj_phi2 are called need to be reworked to cut out the use of Snyder's t = exp(-psi() variable and instead use sinh(psi). Maybe include the machinery for series conversions between all auxiliary latitudes as "support functions". Then etmerc could use this (as could mlfn for computing meridional distance). merc could offer the etmerc style projection via chi as an option when the flattening is sufficiently small. --- src/phi2.cpp | 175 ++++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 120 insertions(+), 55 deletions(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index b81456b0..eb6d5c82 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -1,6 +1,7 @@ /* Determine latitude angle phi-2. */ #include +#include #include "proj.h" #include "proj_internal.h" @@ -8,61 +9,125 @@ static const double TOL = 1.0e-10; static const int N_ITER = 15; +double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { + /**************************************************************************** + * Convert tau' = sinh(psi) = tan(chi) to tau = tan(phi). The code is taken + * from GeographicLib::Math::tauf(taup, e). + * + * Here + * phi = geographic latitude (radians) + * psi is the isometric latitude + * psi = asinh(tan(phi)) - e * atanh(e * sin(phi)) + * = asinh(tan(chi)) + * chi is the conformal latitude + * + * The representation of latitudes via their tangents, tan(phi) and tan(chi), + * maintains full *relative* accuracy close to latitude = 0 and +/- pi/2. + * This is sometimes important, e.g., to compute the scale of the transverse + * Mercator projection which involves cos(phi)/cos(chi) tan(phi) + * + * From Karney (2011), Eq. 7, + * + * tau' = sinh(psi) = sinh(asinh(tan(phi)) - e * atanh(e * sin(phi))) + * = tan(phi) * cosh(e * atanh(e * sin(phi))) - + * sec(phi) * sinh(e * atanh(e * sin(phi))) + * = tau * sqrt(1 + sigma^2) - sqrt(1 + tau^2) * sigma + * where + * sigma = sinh(e * atanh( e * tau / sqrt(1 + tau^2) )) + * + * For e small, + * + * tau' = (1 - e^2) * tau + * + * The relation tau'(tau) can therefore by reliably inverted by Newton's + * method with + * + * tau = tau' / (1 - e^2) + * + * as an initial guess. Newton's method requires dtau'/dtau. Noting that + * + * dsigma/dtau = e^2 * sqrt(1 + sigma^2) / + * (sqrt(1 + tau^2) * (1 + (1 - e^2) * tau^2)) + * d(sqrt(1 + tau^2))/dtau = tau / sqrt(1 + tau^2) + * + * we have + * + * dtau'/dtau = (1 - e^2) * sqrt(1 + tau'^2) * sqrt(1 + tau^2) / + * (1 + (1 - e^2) * tau^2) + * + * This works fine unless tau^2 and tau'^2 overflows. This may be partially + * cured by writing, e.g., sqrt(1 + tau^2) as hypot(1, tau). However, nan + * will still be generated with tau' = inf, since (inf - inf) will appear in + * the Newton iteration. + * + * If we note that for sufficiently large |tau|, i.e., |tau| >= 2/sqrt(eps), + * sqrt(1 + tau^2) = |tau| and + * + * tau' = exp(- e * atanh(e)) * tau + * + * So + * + * tau = exp(e * atanh(e)) * tau' + * + * can be returned unless |tau| >= 2/sqrt(eps); this then avoids overflow + * problems for large tau' and returns the correct result for tau' = +/-inf + * and nan. + * + * Newton's method usually take 2 iterations to converge to double precision + * accuracy (for WGS84 flattening). However only 1 iteration is needed for + * |chi| < 3.35 deg. In addition, only 1 iteration is needed for |chi| > + * 89.18 deg (tau' > 70), if tau = exp(e * atanh(e)) * tau' is used as the + * starting guess. + ****************************************************************************/ + + constexpr int numit = 5; + // min iterations = 1, max iterations = 2; mean = 1.954 + constexpr double tol = sqrt(std::numeric_limits::epsilon()) / 10; + constexpr double tmax = 2 / sqrt(std::numeric_limits::epsilon()); + double + e2m = 1 - e * e, + tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m, + stol = tol * std::max(1.0, fabs(taup)); + if (!(fabs(tau) < tmax)) return tau; // handles +/-inf and nan and e = 1 + // if (e2m < 0) return std::numeric_limits::quiet_NaN(); + int i = numit; + for (; i; --i) { + double tau1 = sqrt(1 + tau * tau), + sig = sinh( e * atanh(e * tau / tau1) ), + taupa = sqrt(1 + sig * sig) * tau - sig * tau1, + dtau = (taup - taupa) * (1 + e2m * tau * tau) / + ( e2m * sqrt(1 + tau * tau) * sqrt(1 + taupa * taupa) ); + tau += dtau; + if (!(fabs(dtau) >= stol)) + break; + } + if (i == 0) + pj_ctx_set_errno(ctx, PJD_ERR_NON_CONV_SINHPSI2TANPHI); + return tau; +} + /*****************************************************************************/ double pj_phi2(projCtx ctx, const double ts0, const double e) { -/****************************************************************************** -Determine latitude angle phi-2. -Inputs: - ts = exp(-psi) where psi is the isometric latitude (dimensionless) - e = eccentricity of the ellipsoid (dimensionless) -Output: - phi = geographic latitude (radians) -Here isometric latitude is defined by - psi = log( tan(pi/4 + phi/2) * - ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) ) - = asinh(tan(phi)) - e * atanh(e * sin(phi)) -This routine inverts this relation using the iterative scheme given -by Snyder (1987), Eqs. (7-9) - (7-11) -*******************************************************************************/ - const double eccnth = .5 * e; - double ts = ts0; -#ifdef no_longer_used_original_convergence_on_exact_dphi - double Phi = M_HALFPI - 2. * atan(ts); -#endif - int i = N_ITER; - - for(;;) { - /* - * sin(Phi) = sin(PI/2 - 2* atan(ts)) - * = cos(2*atan(ts)) - * = 2*cos^2(atan(ts)) - 1 - * = 2 / (1 + ts^2) - 1 - * = (1 - ts^2) / (1 + ts^2) - */ - const double sinPhi = (1 - ts * ts) / (1 + ts * ts); - const double con = e * sinPhi; - double old_ts = ts; - ts = ts0 * pow((1. - con) / (1. + con), eccnth); -#ifdef no_longer_used_original_convergence_on_exact_dphi - /* The convergence criterion is nominally on exact dphi */ - const double newPhi = M_HALFPI - 2. * atan(ts); - const double dphi = newPhi - Phi; - Phi = newPhi; -#else - /* If we don't immediately update Phi from this, we can - * change the conversion criterion to save us computing atan() at each step. - * Particularly we can observe that: - * |atan(ts) - atan(old_ts)| <= |ts - old_ts| - * So if |ts - old_ts| matches our convergence criterion, we're good. - */ - const double dphi = 2 * (ts - old_ts); -#endif - if (fabs(dphi) > TOL && --i) { - continue; - } - break; - } - if (i <= 0) - pj_ctx_set_errno(ctx, PJD_ERR_NON_CON_INV_PHI2); - return M_HALFPI - 2. * atan(ts); + /**************************************************************************** + * Determine latitude angle phi-2. + * Inputs: + * ts = exp(-psi) where psi is the isometric latitude (dimensionless) + * e = eccentricity of the ellipsoid (dimensionless) + * Output: + * phi = geographic latitude (radians) + * Here isometric latitude is defined by + * psi = log( tan(pi/4 + phi/2) * + * ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) ) + * = asinh(tan(phi)) - e * atanh(e * sin(phi)) + * + * OLD: This routine inverts this relation using the iterative scheme given + * by Snyder (1987), Eqs. (7-9) - (7-11). + * + * NEW: This routine writes converts t = exp(-psi) to + * + * tau' = sinh(psi) = (1/t - t)/2 + * + * returns atan(sinpsi2tanphi(tau')) + ***************************************************************************/ + return atan(pj_sinhpsi2tanphi(ctx, (1/ts0 - ts0) / 2, e)); } -- cgit v1.2.3 From b3a20c5b9c6efbeb0d6f528aad1d4e6bb4332bfa Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Mon, 26 Oct 2020 13:58:42 -0400 Subject: Try to fix compiler complaints for max and constexpr sqrt --- src/phi2.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index eb6d5c82..6e810240 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -2,6 +2,7 @@ #include #include +#include #include "proj.h" #include "proj_internal.h" @@ -82,8 +83,8 @@ double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { constexpr int numit = 5; // min iterations = 1, max iterations = 2; mean = 1.954 - constexpr double tol = sqrt(std::numeric_limits::epsilon()) / 10; - constexpr double tmax = 2 / sqrt(std::numeric_limits::epsilon()); + static const double tol = sqrt(std::numeric_limits::epsilon()) / 10; + static const double tmax = 2 / sqrt(std::numeric_limits::epsilon()); double e2m = 1 - e * e, tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m, -- cgit v1.2.3 From 6c2363e88f8eb86e1331b3dc5f05897436ad9c43 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Mon, 26 Oct 2020 14:57:52 -0400 Subject: phi2.cpp: remove unused static consts + minor code tweak --- src/phi2.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index 6e810240..c60ca055 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -7,9 +7,6 @@ #include "proj.h" #include "proj_internal.h" -static const double TOL = 1.0e-10; -static const int N_ITER = 15; - double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { /**************************************************************************** * Convert tau' = sinh(psi) = tan(chi) to tau = tan(phi). The code is taken @@ -97,7 +94,7 @@ double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { sig = sinh( e * atanh(e * tau / tau1) ), taupa = sqrt(1 + sig * sig) * tau - sig * tau1, dtau = (taup - taupa) * (1 + e2m * tau * tau) / - ( e2m * sqrt(1 + tau * tau) * sqrt(1 + taupa * taupa) ); + ( e2m * tau1 * sqrt(1 + taupa * taupa) ); tau += dtau; if (!(fabs(dtau) >= stol)) break; -- cgit v1.2.3 From ee35d15db597801a86314dc8a47da7de10c9d8f2 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Mon, 26 Oct 2020 17:50:51 -0400 Subject: phi2.cpp: Slight cosmetic changes to sinpsi2tanphi. --- src/phi2.cpp | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index c60ca055..7bbb5679 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -80,23 +80,29 @@ double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { constexpr int numit = 5; // min iterations = 1, max iterations = 2; mean = 1.954 - static const double tol = sqrt(std::numeric_limits::epsilon()) / 10; - static const double tmax = 2 / sqrt(std::numeric_limits::epsilon()); - double + static const double + rooteps = sqrt(std::numeric_limits::epsilon()), + tol = rooteps / 10, // the convergence criterion for Newton's method + tmax = 2 / rooteps; // the large arg limit is exact for tau > tmax + const double e2m = 1 - e * e, - tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m, stol = tol * std::max(1.0, fabs(taup)); - if (!(fabs(tau) < tmax)) return tau; // handles +/-inf and nan and e = 1 + // The initial guess. 70 corresponds to chi = 89.18 deg (see above) + double tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m; + if (!(fabs(tau) < tmax)) // handles +/-inf and nan and e = 1 + return tau; + // If we need to deal with e > 1, then we could include: // if (e2m < 0) return std::numeric_limits::quiet_NaN(); int i = numit; for (; i; --i) { - double tau1 = sqrt(1 + tau * tau), + double + tau1 = sqrt(1 + tau * tau), sig = sinh( e * atanh(e * tau / tau1) ), taupa = sqrt(1 + sig * sig) * tau - sig * tau1, - dtau = (taup - taupa) * (1 + e2m * tau * tau) / - ( e2m * tau1 * sqrt(1 + taupa * taupa) ); + dtau = ( (taup - taupa) * (1 + e2m * (tau * tau)) / + (e2m * tau1 * sqrt(1 + taupa * taupa)) ); tau += dtau; - if (!(fabs(dtau) >= stol)) + if (!(fabs(dtau) >= stol)) // backwards test to allow nans to succeed. break; } if (i == 0) -- cgit v1.2.3 From 6bd7c777f8e789f8ea34a6aa68104ab44a31beee Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Mon, 26 Oct 2020 19:03:26 -0400 Subject: Fix/add some comments. --- src/phi2.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index 7bbb5679..b9b37765 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -127,7 +127,7 @@ double pj_phi2(projCtx ctx, const double ts0, const double e) { * OLD: This routine inverts this relation using the iterative scheme given * by Snyder (1987), Eqs. (7-9) - (7-11). * - * NEW: This routine writes converts t = exp(-psi) to + * NEW: This routine converts t = exp(-psi) to * * tau' = sinh(psi) = (1/t - t)/2 * -- cgit v1.2.3 From 94e36270ca393bd7b107bf690f09fd8ec1cd046b Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Tue, 27 Oct 2020 10:02:27 -0400 Subject: Use nm units in builtins.gie. Remove backward looking comments in code. --- src/phi2.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index b9b37765..0fdca47c 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -116,6 +116,7 @@ double pj_phi2(projCtx ctx, const double ts0, const double e) { * Determine latitude angle phi-2. * Inputs: * ts = exp(-psi) where psi is the isometric latitude (dimensionless) + * this variable is defined in Snyder (1987), Eq. (7-10) * e = eccentricity of the ellipsoid (dimensionless) * Output: * phi = geographic latitude (radians) @@ -123,13 +124,12 @@ double pj_phi2(projCtx ctx, const double ts0, const double e) { * psi = log( tan(pi/4 + phi/2) * * ( (1 - e*sin(phi)) / (1 + e*sin(phi)) )^(e/2) ) * = asinh(tan(phi)) - e * atanh(e * sin(phi)) + * = asinh(tan(chi)) + * chi = conformal latitude * - * OLD: This routine inverts this relation using the iterative scheme given - * by Snyder (1987), Eqs. (7-9) - (7-11). + * This routine converts t = exp(-psi) to * - * NEW: This routine converts t = exp(-psi) to - * - * tau' = sinh(psi) = (1/t - t)/2 + * tau' = tan(chi) = sinh(psi) = (1/t - t)/2 * * returns atan(sinpsi2tanphi(tau')) ***************************************************************************/ -- cgit v1.2.3 From 692fc26b6d494aeaa85658314bc020a5cd6da7a1 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Tue, 27 Oct 2020 16:44:43 -0400 Subject: merc.cpp + phi2.cpp: Avoid declaring multiple variables in 1 statement. --- src/phi2.cpp | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index 0fdca47c..1c48d67f 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -80,13 +80,11 @@ double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { constexpr int numit = 5; // min iterations = 1, max iterations = 2; mean = 1.954 - static const double - rooteps = sqrt(std::numeric_limits::epsilon()), - tol = rooteps / 10, // the convergence criterion for Newton's method - tmax = 2 / rooteps; // the large arg limit is exact for tau > tmax - const double - e2m = 1 - e * e, - stol = tol * std::max(1.0, fabs(taup)); + static const double rooteps = sqrt(std::numeric_limits::epsilon()); + static const double tol = rooteps / 10; // the criterion for Newton's method + static const double tmax = 2 / rooteps; // threshold for large arg limit exact + const double e2m = 1 - e * e; + const double stol = tol * std::max(1.0, fabs(taup)); // The initial guess. 70 corresponds to chi = 89.18 deg (see above) double tau = fabs(taup) > 70 ? taup * exp(e * atanh(e)) : taup / e2m; if (!(fabs(tau) < tmax)) // handles +/-inf and nan and e = 1 @@ -95,12 +93,11 @@ double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { // if (e2m < 0) return std::numeric_limits::quiet_NaN(); int i = numit; for (; i; --i) { - double - tau1 = sqrt(1 + tau * tau), - sig = sinh( e * atanh(e * tau / tau1) ), - taupa = sqrt(1 + sig * sig) * tau - sig * tau1, - dtau = ( (taup - taupa) * (1 + e2m * (tau * tau)) / - (e2m * tau1 * sqrt(1 + taupa * taupa)) ); + double tau1 = sqrt(1 + tau * tau); + double sig = sinh( e * atanh(e * tau / tau1) ); + double taupa = sqrt(1 + sig * sig) * tau - sig * tau1; + double dtau = ( (taup - taupa) * (1 + e2m * (tau * tau)) / + (e2m * tau1 * sqrt(1 + taupa * taupa)) ); tau += dtau; if (!(fabs(dtau) >= stol)) // backwards test to allow nans to succeed. break; -- cgit v1.2.3 From aa46197d66ce70ece382bf955326c46b13f35864 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sat, 14 Nov 2020 22:55:31 +0100 Subject: Weed out proj_api.h datatypes and replace them with their proj.h counterparts --- src/phi2.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index 1c48d67f..214d1058 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -7,7 +7,7 @@ #include "proj.h" #include "proj_internal.h" -double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { +double pj_sinhpsi2tanphi(PJ_CONTEXT *ctx, const double taup, const double e) { /**************************************************************************** * Convert tau' = sinh(psi) = tan(chi) to tau = tan(phi). The code is taken * from GeographicLib::Math::tauf(taup, e). @@ -108,7 +108,7 @@ double pj_sinhpsi2tanphi(projCtx ctx, const double taup, const double e) { } /*****************************************************************************/ -double pj_phi2(projCtx ctx, const double ts0, const double e) { +double pj_phi2(PJ_CONTEXT *ctx, const double ts0, const double e) { /**************************************************************************** * Determine latitude angle phi-2. * Inputs: -- cgit v1.2.3 From a74b985b5006c2d279353a245cfcb850cf7fcc94 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Tue, 17 Nov 2020 12:54:24 +0100 Subject: Remove pj_ctx_* functions and use their proj_context counterparts --- src/phi2.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/phi2.cpp') diff --git a/src/phi2.cpp b/src/phi2.cpp index 214d1058..2f258e47 100644 --- a/src/phi2.cpp +++ b/src/phi2.cpp @@ -103,7 +103,7 @@ double pj_sinhpsi2tanphi(PJ_CONTEXT *ctx, const double taup, const double e) { break; } if (i == 0) - pj_ctx_set_errno(ctx, PJD_ERR_NON_CONV_SINHPSI2TANPHI); + proj_context_errno_set(ctx, PJD_ERR_NON_CONV_SINHPSI2TANPHI); return tau; } -- cgit v1.2.3