From 3d000ce54bbea0da092df0c66f175065fde04dcc Mon Sep 17 00:00:00 2001
From: ClemensFischer <fischer.clemens@gmail.com>
Date: Fri, 23 Jan 2026 20:37:17 +0100
Subject: [PATCH] Updated TransverseMercatorProjection projections

---
 .../Shared/TransverseMercatorProjection.cs    | 70 ++++++++++++++-----
 .../TransverseMercatorProjectionSnyder.cs     | 10 +--
 2 files changed, 57 insertions(+), 23 deletions(-)

diff --git a/MapControl/Shared/TransverseMercatorProjection.cs b/MapControl/Shared/TransverseMercatorProjection.cs
index 365ad75c..d2bf8241 100644
--- a/MapControl/Shared/TransverseMercatorProjection.cs
+++ b/MapControl/Shared/TransverseMercatorProjection.cs
@@ -63,13 +63,6 @@ namespace MapControl
 
         public override Point? LocationToMap(double latitude, double longitude)
         {
-#if NETFRAMEWORK
-            static double Atanh(double x) => Math.Log((1d + x) / (1d - x)) / 2d;
-#else
-            static double Atanh(double x) => Math.Atanh(x);
-#endif
-            // k0 * A
-            var k0A = ScaleFactor * EquatorialRadius * f1;
             // φ
             var phi = latitude * Math.PI / 180d;
             var sinPhi = Math.Sin(phi);
@@ -81,20 +74,20 @@ namespace MapControl
             var xi_ = Math.Atan2(t, Math.Cos(dLambda));
             // η'
             var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t));
-            // ξ
-            var xi = xi_ +
-                a1 * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_) +
-                a2 * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_) +
-                a3 * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_);
-            // η
-            var eta = eta_ +
+            // k0 * A
+            var k0A = ScaleFactor * EquatorialRadius * f1;
+            // E
+            var x = FalseEasting + k0A * (eta_ +
                 a1 * Math.Cos(2d * xi_) * Math.Sinh(2d * eta_) +
                 a2 * Math.Cos(4d * xi_) * Math.Sinh(4d * eta_) +
-                a3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_);
+                a3 * Math.Cos(6d * xi_) * Math.Sinh(6d * eta_));
+            // N
+            var y = FalseNorthing + k0A * (xi_ +
+                a1 * Math.Sin(2d * xi_) * Math.Cosh(2d * eta_) +
+                a2 * Math.Sin(4d * xi_) * Math.Cosh(4d * eta_) +
+                a3 * Math.Sin(6d * xi_) * Math.Cosh(6d * eta_));
 
-            return new Point(
-                k0A * eta + FalseEasting,
-                k0A * xi + FalseNorthing);
+            return new Point(x, y);
         }
 
         public override Location MapToLocation(double x, double y)
@@ -129,5 +122,46 @@ namespace MapControl
                 phi * 180d / Math.PI,
                 dLambda * 180d / Math.PI + CentralMeridian);
         }
+
+#if NETFRAMEWORK
+        private static double Atanh(double x) => Math.Log((1d + x) / (1d - x)) / 2d;
+#else
+        private static double Atanh(double x) => Math.Atanh(x);
+#endif
+        /*
+         *  Relative scale is usually < 1.001 and hence neglectable.
+         * 
+        public override Matrix RelativeScale(double latitude, double longitude)
+        {
+            // φ
+            var phi = latitude * Math.PI / 180d;
+            var sinPhi = Math.Sin(phi);
+            // t
+            var t = Math.Sinh(Atanh(sinPhi) - f2 * Atanh(f2 * sinPhi));
+            // λ - λ0
+            var dLambda = (longitude - CentralMeridian) * Math.PI / 180d;
+            var cosLambda = Math.Cos(dLambda);
+            // ξ'
+            var xi_ = Math.Atan2(t, cosLambda);
+            // η'
+            var eta_ = Atanh(Math.Sin(dLambda) / Math.Sqrt(1d + t * t));
+            // σ
+            var sigma = 1 +
+                2d * a1 * Math.Cos(2d * xi_) * Math.Cosh(2d * eta_) +
+                4d * a2 * Math.Cos(4d * xi_) * Math.Cosh(4d * eta_) +
+                6d * a3 * Math.Cos(6d * xi_) * Math.Cosh(6d * eta_);
+            // τ
+            var tau =
+                2d * a1 * Math.Sin(2d * xi_) * Math.Sinh(2d * eta_) +
+                4d * a2 * Math.Sin(4d * xi_) * Math.Sinh(4d * eta_) +
+                6d * a3 * Math.Sin(6d * xi_) * Math.Sinh(6d * eta_);
+
+            // h = k/k0 for relative scale
+            var n = Flattening / (2d - Flattening);
+            var u = (1d - n) / (1d + n) * Math.Tan(phi);
+            var h = f1 * Math.Sqrt((1d + u * u) * (sigma * sigma + tau * tau) / (t * t + cosLambda * cosLambda));
+
+            return new Matrix(h, 0d, 0d, h, 0d, 0d);
+        }*/
     }
 }
diff --git a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
index f202eee4..881e04be 100644
--- a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
+++ b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
@@ -63,15 +63,15 @@ namespace MapControl
             var e6 = e2 * e4;
 
             return EquatorialRadius *
-                ((1d - e2 / 4d - 3d * e4 / 64d - 5d * e6 / 256d) * phi -
-                (3d * e2 / 8d + 3d * e4 / 32d + 45d * e6 / 1024d) * Math.Sin(2d * phi) +
-                (15d * e4 / 256d + 45d * e6 / 1024d) * Math.Sin(4d * phi) -
-                35d * e6 / 3072d * Math.Sin(6d * phi)); // (3-21)
+                ((1d - e2 / 4d - 3d / 64d * e4 - 5d / 256d * e6) * phi -
+                (3d / 8d * e2 + 3d / 32d * e4 + 45d / 1024d * e6) * Math.Sin(2d * phi) +
+                (15d / 256d * e4 + 45d / 1024d * e6) * Math.Sin(4d * phi) -
+                35d / 3072d * e6 * Math.Sin(6d * phi)); // (3-21)
         }
 
         public override Matrix RelativeScale(double latitude, double longitude)
         {
-            var k = ScaleFactor;
+            var k = 1d; // omit k0 for relative scale
 
             if (latitude > -90d && latitude < 90d)
             {