From 2e27fc740dff3becd69b452eac343cb482501a3c Mon Sep 17 00:00:00 2001
From: ClemensFischer <fischer.clemens@gmail.com>
Date: Wed, 21 Jan 2026 14:15:26 +0100
Subject: [PATCH] Updated azimuthal projections

---
 .../Shared/AzimuthalEquidistantProjection.cs   | 18 +++++++++++++-----
 MapControl/Shared/AzimuthalProjection.cs       | 18 +++++++++++++++++-
 MapControl/Shared/GnomonicProjection.cs        |  4 +---
 MapControl/Shared/OrthographicProjection.cs    |  5 +----
 MapControl/Shared/StereographicProjection.cs   |  2 +-
 5 files changed, 33 insertions(+), 14 deletions(-)

diff --git a/MapControl/Shared/AzimuthalEquidistantProjection.cs b/MapControl/Shared/AzimuthalEquidistantProjection.cs
index 279769b9..1649b54b 100644
--- a/MapControl/Shared/AzimuthalEquidistantProjection.cs
+++ b/MapControl/Shared/AzimuthalEquidistantProjection.cs
@@ -30,28 +30,36 @@ namespace MapControl
         {
             var p = GetProjectedPoint(latitude, longitude);
 
-            if (Math.Abs(p.CosC) >= 1d)
+            if (p.CosC == 1d)
             {
                 return new Matrix(1d, 0d, 0d, 1d, 0d, 0d);
             }
 
+            if (p.CosC == -1d)
+            {
+                return new Matrix(1d, 0d, 0d, double.PositiveInfinity, 0d, 0d);
+            }
+
             var c = Math.Acos(p.CosC);
             var k = c / Math.Sin(c); // p.195 (25-2)
 
-            var scale = new Matrix(1d, 0d, 0d, k, 0d, 0d); // h == 1
-            scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI);
-            return scale;
+            return p.RelativeScale(1d, k);
         }
 
         public override Point? LocationToMap(double latitude, double longitude)
         {
             var p = GetProjectedPoint(latitude, longitude);
 
-            if (Math.Abs(p.CosC) >= 1d)
+            if (p.CosC == 1d) // p.195 "If cos c = 1, ... k' = 1, and x = y = 0."
             {
                 return new Point();
             }
 
+            if (p.CosC == -1)
+            {
+                return null; // p.195 "If cos c = -1, the point ... is plotted as a circle of radius πR."
+            }
+
             var c = Math.Acos(p.CosC);
             var k = c / Math.Sin(c); // p.195 (25-2)
 
diff --git a/MapControl/Shared/AzimuthalProjection.cs b/MapControl/Shared/AzimuthalProjection.cs
index e31d8874..dc9a3bf1 100644
--- a/MapControl/Shared/AzimuthalProjection.cs
+++ b/MapControl/Shared/AzimuthalProjection.cs
@@ -1,4 +1,7 @@
 using System;
+#if WPF
+using System.Windows.Media;
+#endif
 
 namespace MapControl
 {
@@ -32,10 +35,23 @@ namespace MapControl
                 var sinPhi1 = Math.Sin(phi1);
                 var cosLambda = Math.Cos(dLambda);
                 var sinLambda = Math.Sin(dLambda);
+                var cosC = sinPhi1 * sinPhi + cosPhi1 * cosPhi * cosLambda; // (5-3)
 
                 X = cosPhi * sinLambda;
                 Y = cosPhi1 * sinPhi - sinPhi1 * cosPhi * cosLambda;
-                CosC = sinPhi1 * sinPhi + cosPhi1 * cosPhi * cosLambda; // (5-3)
+                CosC = Math.Min(Math.Max(cosC, -1d), 1d); // protect against rounding errors
+            }
+
+            public Matrix RelativeScale(double radialScale, double perpendicularScale)
+            {
+                var scale = new Matrix(radialScale, 0d, 0d, perpendicularScale, 0d, 0d);
+
+                if (radialScale != perpendicularScale)
+                {
+                    scale.Rotate(-Math.Atan2(Y, X) * 180d / Math.PI);
+                }
+
+                return scale;
             }
         }
 
diff --git a/MapControl/Shared/GnomonicProjection.cs b/MapControl/Shared/GnomonicProjection.cs
index 17917f6e..f20d3e5f 100644
--- a/MapControl/Shared/GnomonicProjection.cs
+++ b/MapControl/Shared/GnomonicProjection.cs
@@ -32,9 +32,7 @@ namespace MapControl
             var k = 1d / p.CosC; // p.165 (22-3)
             var h = k * k; // p.165 (22-2)
 
-            var scale = new Matrix(h, 0d, 0d, k, 0d, 0d);
-            scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI);
-            return scale;
+            return p.RelativeScale(h, k);
         }
 
         public override Point? LocationToMap(double latitude, double longitude)
diff --git a/MapControl/Shared/OrthographicProjection.cs b/MapControl/Shared/OrthographicProjection.cs
index a7de14e6..5a76415f 100644
--- a/MapControl/Shared/OrthographicProjection.cs
+++ b/MapControl/Shared/OrthographicProjection.cs
@@ -29,11 +29,8 @@ namespace MapControl
         public override Matrix RelativeScale(double latitude, double longitude)
         {
             var p = GetProjectedPoint(latitude, longitude);
-            var h = p.CosC; // p.149 (20-5)
 
-            var scale = new Matrix(h, 0d, 0d, 1d, 0d, 0d); // k == 1
-            scale.Rotate(-Math.Atan2(p.Y, p.X) * 180d / Math.PI);
-            return scale;
+            return p.RelativeScale(p.CosC, 1d); // p.149 (20-5), k == 1
         }
 
         public override Point? LocationToMap(double latitude, double longitude)
diff --git a/MapControl/Shared/StereographicProjection.cs b/MapControl/Shared/StereographicProjection.cs
index b166c489..66b23f6d 100644
--- a/MapControl/Shared/StereographicProjection.cs
+++ b/MapControl/Shared/StereographicProjection.cs
@@ -31,7 +31,7 @@ namespace MapControl
             var p = GetProjectedPoint(latitude, longitude);
             var k = 2d / (1d + p.CosC); // p.157 (21-4), k0 == 1
 
-            return new Matrix(k, 0d, 0d, k, 0d, 0d);
+            return p.RelativeScale(k, k);
         }
 
         public override Point? LocationToMap(double latitude, double longitude)