diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs
index 70c93925..3800c6d6 100644
--- a/MapControl/Shared/StereographicProjection.cs
+++ b/MapControl/Shared/StereographicProjection.cs
@@ -27,10 +27,25 @@ namespace MapControl
 
         public override double GridConvergence(double latitude, double longitude)
         {
-            var p0 = LocationToMap(latitude, longitude);
-            var p1 = LocationToMap(latitude + 1e-3, longitude);
+            var phi0 = LatitudeOfOrigin * Math.PI / 180d; // φ1
+            var phi1 = latitude * Math.PI / 180d;
+            var phi2 = (latitude + 1e-3) * Math.PI / 180d;
+            var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; // λ - λ0
+            var sinPhi0 = Math.Sin(phi0);
+            var cosPhi0 = Math.Cos(phi0);
+            var sinPhi1 = Math.Sin(phi1);
+            var cosPhi1 = Math.Cos(phi1);
+            var sinPhi2 = Math.Sin(phi2);
+            var cosPhi2 = Math.Cos(phi2);
+            var sinLambda = Math.Sin(dLambda);
+            var cosLambda = Math.Cos(dLambda);
+            var k1 = 2d / (1d + sinPhi0 * sinPhi1 + cosPhi0 * cosPhi1 * cosLambda);
+            var k2 = 2d / (1d + sinPhi0 * sinPhi2 + cosPhi0 * cosPhi2 * cosLambda);
+            var dCosPhi = k2 * cosPhi2 - k1 * cosPhi1;
+            var dSinPhi = k2 * sinPhi2 - k1 * sinPhi1;
 
-            return Math.Atan2(p0.X - p1.X, p1.Y - p0.Y) * 180d / Math.PI;
+            return Math.Atan2(-sinLambda * dCosPhi,
+                cosPhi0 * dSinPhi - sinPhi0 * cosLambda * dCosPhi) * 180d / Math.PI;
         }
 
         public override Matrix RelativeTransform(double latitude, double longitude)
@@ -40,20 +55,18 @@ namespace MapControl
 
         public override Point LocationToMap(double latitude, double longitude)
         {
-            var phi = latitude * Math.PI / 180d;
-            var phi1 = LatitudeOfOrigin * Math.PI / 180d;
+            var phi0 = LatitudeOfOrigin * Math.PI / 180d; // φ1
+            var phi = latitude * Math.PI / 180d; // φ
             var dLambda = (longitude - CentralMeridian) * Math.PI / 180d; // λ - λ0
-            var cosPhi = Math.Cos(phi);
+            var sinPhi0 = Math.Sin(phi0);
+            var cosPhi0 = Math.Cos(phi0);
             var sinPhi = Math.Sin(phi);
-            var cosPhi1 = Math.Cos(phi1);
-            var sinPhi1 = Math.Sin(phi1);
-            var cosLambda = Math.Cos(dLambda);
+            var cosPhi = Math.Cos(phi);
             var sinLambda = Math.Sin(dLambda);
+            var cosPhiCosLambda = cosPhi * Math.Cos(dLambda);
             var x = cosPhi * sinLambda;
-            var y = cosPhi1 * sinPhi - sinPhi1 * cosPhi * cosLambda;
-            var cosC = sinPhi1 * sinPhi + cosPhi1 * cosPhi * cosLambda; // (5-3)
-            cosC = Math.Min(Math.Max(cosC, -1d), 1d); // protect against rounding errors
-            var k = 2d / (1d + cosC); // p.157 (21-4), k0 == 1
+            var y = cosPhi0 * sinPhi - sinPhi0 * cosPhiCosLambda;
+            var k = 2d / (1d + sinPhi0 * sinPhi + cosPhi0 * cosPhiCosLambda); // p.157 (21-4), k0 == 1
 
             return new Point(
                 EquatorialRadius * k * x,
@@ -67,10 +80,10 @@ namespace MapControl
             var cosC = Math.Cos(c);
             var sinC = Math.Sin(c);
 
-            var phi1 = LatitudeOfOrigin * Math.PI / 180d;
-            var cosPhi1 = Math.Cos(phi1);
-            var sinPhi1 = Math.Sin(phi1);
-            var phi = Math.Asin(cosC * sinPhi1 + y * sinC * cosPhi1 / rho); // (20-14)
+            var phi0 = LatitudeOfOrigin * Math.PI / 180d; // φ1
+            var cosPhi0 = Math.Cos(phi0);
+            var sinPhi0 = Math.Sin(phi0);
+            var phi = Math.Asin(cosC * sinPhi0 + y * sinC * cosPhi0 / rho); // (20-14)
             double u, v;
 
             if (LatitudeOfOrigin == 90d) // (20-16)
@@ -86,7 +99,7 @@ namespace MapControl
             else // (20-15)
             {
                 u = x * sinC;
-                v = rho * cosPhi1 * cosC - y * sinPhi1 * sinC;
+                v = rho * cosPhi0 * cosC - y * sinPhi0 * sinC;
             }
 
             return new Location(
diff --git a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
index ae28a745..de62674e 100644
--- a/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
+++ b/MapControl/Shared/TransverseMercatorProjectionSnyder.cs
@@ -120,20 +120,20 @@ namespace MapControl
             var M0 = MeridianDistance(LatitudeOfOrigin * Math.PI / 180d);
             var M = M0 + (y - FalseNorthing) / ScaleFactor; // (8-20)
             var mu = M / (EquatorialRadius * (1d - e2 / 4d - e4 * 3d / 64d - e6 * 5d / 256d)); // (7-19)
-            var phi1 = mu +
+            var phi0 = mu +
                 (e1 * 3d / 2d - e13 * 27d / 32d) * Math.Sin(2d * mu) +
                 (e12 * 21d / 16d - e14 * 55d / 32d) * Math.Sin(4d * mu) +
                  e13 * 151d / 96d * Math.Sin(6d * mu) +
                  e14 * 1097d / 512d * Math.Sin(8d * mu); // (3-26)
 
-            var sinPhi1 = Math.Sin(phi1);
-            var cosPhi1 = Math.Cos(phi1);
-            var tanPhi1 = sinPhi1 / cosPhi1;
+            var sinPhi0 = Math.Sin(phi0);
+            var cosPhi0 = Math.Cos(phi0);
+            var tanPhi0 = sinPhi0 / cosPhi0;
 
             var e_2 = e2 / (1d - e2); // (8-12)
-            var C1 = e_2 * cosPhi1 * cosPhi1; // (8-21)
-            var T1 = sinPhi1 * sinPhi1 / (cosPhi1 * cosPhi1); // (8-22)
-            s = Math.Sqrt(1d - e2 * sinPhi1 * sinPhi1);
+            var C1 = e_2 * cosPhi0 * cosPhi0; // (8-21)
+            var T1 = sinPhi0 * sinPhi0 / (cosPhi0 * cosPhi0); // (8-22)
+            s = Math.Sqrt(1d - e2 * sinPhi0 * sinPhi0);
             var N1 = EquatorialRadius / s; // (8-23)
             var R1 = EquatorialRadius * (1d - e2) / (s * s * s); // (8-24)
             var D = (x - FalseEasting) / (N1 * ScaleFactor); // (8-25)
@@ -143,11 +143,11 @@ namespace MapControl
             var D5 = D * D4;
             var D6 = D * D5;
 
-            var phi = phi1 - N1 * tanPhi1 / R1 * (D2 / 2d - (5d + 3d * T1 + 10d * C1 - 4d * C1 * C1 - 9d * e_2) * D4 / 24d +
+            var phi = phi0 - N1 * tanPhi0 / R1 * (D2 / 2d - (5d + 3d * T1 + 10d * C1 - 4d * C1 * C1 - 9d * e_2) * D4 / 24d +
                 (61d + 90d * T1 + 45d * T1 * T1 + 298 * C1 - 3d * C1 * C1 - 252d * e_2) * D6 / 720d); // (8-17)
 
             var dLambda = (D - (1d + 2d * T1 + C1) * D3 / 6d +
-                (5d - 2d * C1 - 3d * C1 * C1 + 28d * T1 + 24d * T1 * T1 + 8d * e_2) * D5 / 120d) / cosPhi1; // (8-18)
+                (5d - 2d * C1 - 3d * C1 * C1 + 28d * T1 + 24d * T1 * T1 + 8d * e_2) * D5 / 120d) / cosPhi0; // (8-18)
 
             return new Location(
                 phi * 180d / Math.PI,