Merge pull request #2397 from cffk/merc-update

Update Mercator projection, more accurate, faster

2 files changed, 136 insertions, 67 deletions

diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index cfce5041..def30206 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -1423,9 +1423,9 @@ expect  -0.001790143 -0.000895247
 ===============================================================================
 
 -------------------------------------------------------------------------------
-operation +proj=etmerc   +ellps=GRS80 +zone=30
+operation +proj=etmerc   +ellps=GRS80
 -------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 50 nm
 accept  2 1
 expect  222650.796797586 110642.229411933
 accept  2 -1
@@ -1434,17 +1434,28 @@ accept  -2 1
 expect  -222650.796797586 110642.229411933
 accept  -2 -1
 expect  -222650.796797586 -110642.229411933
+# near pole
+accept  30 89.9999
+expect  5.584698978 10001956.056248082
+# 3900 km from central meridian
+accept  44.69 35.37
+expect  4168136.489446198 4985511.302287407
 
 direction inverse
 accept  200 100
-expect  0.001796631 0.000904369
+expect  0.00179663056816 0.00090436947663
 accept  200 -100
-expect  0.001796631 -0.000904369
+expect  0.00179663056816 -0.00090436947663
 accept  -200 100
-expect  -0.001796631 0.000904369
+expect  -0.00179663056816 0.00090436947663
 accept  -200 -100
-expect  -0.001796631 -0.000904369
-
+expect  -0.00179663056816 -0.00090436947663
+# near pole
+accept  6 1.0001e7
+expect  0.35596960759234 89.99135362646302
+# 3900 km from central meridian
+accept  4168136.489446198 4985511.302287407
+expect  44.69 35.37
 
 ===============================================================================
 # Fahey
@@ -3355,30 +3366,54 @@ expect  -0.001953415 -0.000820580
 -------------------------------------------------------------------------------
 operation +proj=merc   +ellps=GRS80
 -------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 0 m
+accept  0 0
+expect  0 0
+tolerance 50 nm
 accept  2 1
-expect  222638.981586547 110579.965218250
+expect  222638.981586547 110579.965218249
 accept  2 -1
 expect  222638.981586547 -110579.965218249
 accept  -2 1
-expect  -222638.981586547 110579.965218250
+expect  -222638.981586547 110579.965218249
 accept  -2 -1
 expect  -222638.981586547 -110579.965218249
+# inflate tolerance by scale (k = 5.7e15)
+tolerance 3e8
+accept  0 89.99999999999999
+expect  0 235805185.015130176
+accept  0 -89.99999999999999
+expect  0 -235805185.015130176
 
 direction inverse
+tolerance 0 m
+accept  0 0
+expect  0 0
+tolerance 50 nm
 accept  200 100
-expect  0.001796631 0.000904369
+expect  0.00179663056824 0.00090436947704
 accept  200 -100
-expect  0.001796631 -0.000904369
+expect  0.00179663056824 -0.00090436947704
 accept  -200 100
-expect  -0.001796631 0.000904369
+expect  -0.00179663056824 0.00090436947704
 accept  -200 -100
-expect  -0.001796631 -0.000904369
+expect  -0.00179663056824 -0.00090436947704
+accept  0 235805185.015130176
+expect  0 89.99999999999999
+accept  0 -235805185.015130176
+expect  0 -89.99999999999999
+accept  0 1e10
+expect  0 90
+accept  0 -1e10
+expect  0 -90
 
 -------------------------------------------------------------------------------
 operation +proj=merc   +R=6400000
 -------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 0 m
+accept  0 0
+expect  0 0
+tolerance 50 nm
 accept  2 1
 expect  223402.144255274 111706.743574944
 accept  2 -1
@@ -3389,25 +3424,32 @@ accept  -2 -1
 expect  -223402.144255274 -111706.743574944
 
 direction inverse
+tolerance 0 m
+accept  0 0
+expect  0 0
+tolerance 50 nm
 accept  200 100
-expect  0.001790493 0.000895247
+expect  0.00179049310978 0.00089524655486
 accept  200 -100
-expect  0.001790493 -0.000895247
+expect  0.00179049310978 -0.00089524655486
 accept  -200 100
-expect  -0.001790493 0.000895247
+expect  -0.00179049310978 0.00089524655486
 accept  -200 -100
-expect  -0.001790493 -0.000895247
-
+expect  -0.00179049310978 -0.00089524655486
 
 -------------------------------------------------------------------------------
 operation +proj=merc +R=1
 -------------------------------------------------------------------------------
 # Test the numerical stability of the inverse spherical Mercator
 -------------------------------------------------------------------------------
-tolerance   1e-15 m
-accept  0   1e-15
+tolerance   1e-17 m
+accept  0   57.295779513e-15
 expect  0   1e-15
 
+direction inverse
+accept  0   1e-15
+expect  0   57.295779513e-15
+
 
 ===============================================================================
 # Miller Oblated Stereographic
@@ -5658,25 +5700,37 @@ expect  -0.001790143 0.511651393
 -------------------------------------------------------------------------------
 operation +proj=tmerc   +ellps=GRS80
 -------------------------------------------------------------------------------
-tolerance 0.1 mm
+tolerance 50 nm
 accept  2 1
-expect  222650.796795778 110642.229411927
+expect  222650.796797586 110642.229411933
 accept  2 -1
-expect  222650.796795778 -110642.229411927
+expect  222650.796797586 -110642.229411933
 accept  -2 1
-expect  -222650.796795778 110642.229411927
+expect  -222650.796797586 110642.229411933
 accept  -2 -1
-expect  -222650.796795778 -110642.229411927
+expect  -222650.796797586 -110642.229411933
+# near pole
+accept  30 89.9999
+expect  5.584698978 10001956.056248082
+# 3900 km from central meridian
+accept  44.69 35.37
+expect  4168136.489446198 4985511.302287407
 
 direction inverse
 accept  200 100
-expect  0.001796631 0.000904369
+expect  0.00179663056816 0.00090436947663
 accept  200 -100
-expect  0.001796631 -0.000904369
+expect  0.00179663056816 -0.00090436947663
 accept  -200 100
-expect  -0.001796631 0.000904369
+expect  -0.00179663056816 0.00090436947663
 accept  -200 -100
-expect  -0.001796631 -0.000904369
+expect  -0.00179663056816 -0.00090436947663
+# near pole
+accept  6 1.0001e7
+expect  0.35596960759234 89.99135362646302
+# 3900 km from central meridian
+accept  4168136.489446198 4985511.302287407
+expect  44.69 35.37
 
 -------------------------------------------------------------------------------
 operation +proj=tmerc   +R=6400000
diff --git a/test/unit/pj_phi2_test.cpp b/test/unit/pj_phi2_test.cpp
index c4db6e52..fdedae98 100644
--- a/test/unit/pj_phi2_test.cpp
+++ b/test/unit/pj_phi2_test.cpp
@@ -39,47 +39,62 @@ namespace {
 TEST(PjPhi2Test, Basic) {
     projCtx ctx = pj_get_default_ctx();
 
-    EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, 0.0));
-
-    EXPECT_NEAR(0.0, pj_phi2(ctx, 1.0, 0.0), 1e-16);
-    EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, 1.0));
-    EXPECT_DOUBLE_EQ(M_PI, pj_phi2(ctx, -1.0, 0.0));
-    EXPECT_DOUBLE_EQ(M_PI_2, pj_phi2(ctx, 0.0, -1.0));
-
-    EXPECT_NEAR(0.0, pj_phi2(ctx, 1.0, 1.0), 1e-16);
-    EXPECT_DOUBLE_EQ(M_PI, pj_phi2(ctx, -1.0, -1.0));
-
-    // TODO(schwehr): M_PI_4, M_PI_2, M_PI, M_E
-    // https://www.gnu.org/software/libc/manual/html_node/Mathematical-Constants.html
-
-    EXPECT_DOUBLE_EQ(-0.95445818456292697, pj_phi2(ctx, M_PI, 0.0));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, M_PI)));
-    EXPECT_DOUBLE_EQ(4.0960508381527205, pj_phi2(ctx, -M_PI, 0.0));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, -M_PI)));
-
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, M_PI, M_PI)));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -M_PI, -M_PI)));
+    // Expectation is that only sane values of e (and nan is here reckoned to
+    // be sane) are passed to pj_phi2.  Thus the return value with other values
+    // of e is "implementation dependent".
+
+    constexpr auto inf = std::numeric_limits<double>::infinity();
+    constexpr auto nan = std::numeric_limits<double>::quiet_NaN();
+
+    // Strict equality is demanded here.
+    EXPECT_EQ( M_PI_2, pj_phi2(ctx, +0.0, 0.0));
+    EXPECT_EQ(  0.0  , pj_phi2(ctx,  1.0, 0.0));
+    EXPECT_EQ(-M_PI_2, pj_phi2(ctx,  inf, 0.0));
+    // We don't expect pj_phi2 to be called with negative ts (since ts =
+    // exp(-psi)).  However, in the current implementation it is odd in ts.
+    // N.B. ts = +0.0 and ts = -0.0 return different results.
+    EXPECT_EQ(-M_PI_2, pj_phi2(ctx, -0.0, 0.0));
+    EXPECT_EQ(  0.0  , pj_phi2(ctx, -1.0, 0.0));
+    EXPECT_EQ(+M_PI_2, pj_phi2(ctx, -inf, 0.0));
+
+    constexpr double e = 0.2;
+    EXPECT_EQ( M_PI_2, pj_phi2(ctx, +0.0, e));
+    EXPECT_EQ(  0.0  , pj_phi2(ctx,  1.0, e));
+    EXPECT_EQ(-M_PI_2, pj_phi2(ctx,  inf, e));
+    EXPECT_EQ(-M_PI_2, pj_phi2(ctx, -0.0, e));
+    EXPECT_EQ(  0.0  , pj_phi2(ctx, -1.0, e));
+    EXPECT_EQ(+M_PI_2, pj_phi2(ctx, -inf, e));
+
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx,  nan, 0.0)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx,  nan, e  )));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, +0.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx,  1.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx,  inf, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -0.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -1.0, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -inf, nan)));
+    EXPECT_TRUE(std::isnan(pj_phi2(ctx,  nan, nan)));
+
+    EXPECT_DOUBLE_EQ( M_PI/3, pj_phi2(ctx, 1/(sqrt(3.0)+2), 0.0));
+    EXPECT_DOUBLE_EQ( M_PI/4, pj_phi2(ctx, 1/(sqrt(2.0)+1), 0.0));
+    EXPECT_DOUBLE_EQ( M_PI/6, pj_phi2(ctx, 1/ sqrt(3.0)   , 0.0));
+    EXPECT_DOUBLE_EQ(-M_PI/3, pj_phi2(ctx,    sqrt(3.0)+2 , 0.0));
+    EXPECT_DOUBLE_EQ(-M_PI/4, pj_phi2(ctx,    sqrt(2.0)+1 , 0.0));
+    EXPECT_DOUBLE_EQ(-M_PI/6, pj_phi2(ctx,    sqrt(3.0)   , 0.0));
+
+    // Generated with exp(e * atanh(e * sin(phi))) / (tan(phi) + sec(phi))
+    EXPECT_DOUBLE_EQ( M_PI/3, pj_phi2(ctx, 0.27749174377027023413, e));
+    EXPECT_DOUBLE_EQ( M_PI/4, pj_phi2(ctx, 0.42617788119104192995, e));
+    EXPECT_DOUBLE_EQ( M_PI/6, pj_phi2(ctx, 0.58905302448626726064, e));
+    EXPECT_DOUBLE_EQ(-M_PI/3, pj_phi2(ctx,  3.6037108218537833089, e));
+    EXPECT_DOUBLE_EQ(-M_PI/4, pj_phi2(ctx,  2.3464380582241712935, e));
+    EXPECT_DOUBLE_EQ(-M_PI/6, pj_phi2(ctx,  1.6976400399134411849, e));
 }
 
-TEST(PjPhi2Test, AvoidUndefinedBehavior) {
-    auto ctx = pj_get_default_ctx();
+} // namespace
 
-    const auto nan = std::numeric_limits<double>::quiet_NaN();
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, 0.0)));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, nan)));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, nan, nan)));
 
-    // We do not really care about the values that follow.
-    const auto inf = std::numeric_limits<double>::infinity();
 
-    EXPECT_DOUBLE_EQ(-M_PI_2, pj_phi2(ctx, inf, 0.0));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, inf)));
 
-    EXPECT_DOUBLE_EQ(4.7123889803846897, pj_phi2(ctx, -inf, 0.0));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, 0.0, -inf)));
 
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, inf, inf)));
-    EXPECT_TRUE(std::isnan(pj_phi2(ctx, -inf, -inf)));
-}
 
-} // namespace