Polygon - Generate Geo polygon around a specific point (#53)

* Add Polygon class with geospatial methods and tests

Implemented a new static `Polygon` class in the `GeoUK` namespace for generating polygons around geographic points. The class includes input validation and methods for degree-radian conversion.

Also added a suite of unit tests in `PolygonTests.cs` using Xunit to cover various scenarios, ensuring the correctness and performance of the `GeneratePolygonAroundPoint` method.

* Fix longitude normalization in Polygon.cs

Updated longitude normalization to the range of [-180, 180]
instead of [-180, 180). Adjusted logic to correctly handle
values exceeding the bounds by adding or subtracting 360.0
as necessary, ensuring proper wrapping of longitude values.

* Add GeoUK.OSTN.XUnit project and update README.md

- Included a new project `GeoUK.OSTN.XUnit` in the solution.
- Added `pr.yml` to the Solution Items section.
- Enhanced `README.md` with a new section on generating polygons around points, including code examples and a link to GeoJson.io for testing.

GeoUK.OSTN.sln

@@ -18,6 +18,7 @@ EndProject
 Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution Items", "{8EC462FD-D22E-90A8-E5CE-7E832BA40C5D}"
 	ProjectSection(SolutionItems) = preProject
 		.github\workflows\pr.yml = .github\workflows\pr.yml
+		README.md = README.md
 		.github\workflows\release.yml = .github\workflows\release.yml
 	EndProjectSection
 EndProject

GeoUK/Polygon.cs

@@ -0,0 +1,123 @@
+using System;
+using System.Collections.Generic;
+
+namespace GeoUK
+{
+    public static class Polygon
+    {
+        /// <summary>
+        /// WGS84 equatorial radius in meters
+        /// </summary>
+        const double wgs84_equatorial_radius = 6378137.0;
+
+        /// <summary>
+        /// Maximum radius to prevent numerical instability (quarter of Earth's circumference)
+        /// </summary>
+        const double max_radius_meters = Math.PI * wgs84_equatorial_radius / 2.0;
+
+        /// <summary>
+        /// Generates a polygon around a given point.
+        /// Each point is 'radius' from the center
+        /// </summary>
+        /// <param name="longitude">Longitude of the center point</param>
+        /// <param name="latitude">Latitude of the center point</param>
+        /// <param name="radius">Radius of the polygon in meters</param>
+        /// <param name="numberOfPoints">Number of points to make up the polygon (minimum 3)</param>
+        /// <returns>List of [longitude, latitude] coordinate arrays forming the polygon</returns>
+        public static List<double[]> GeneratePolygonAroundPoint(double longitude, double latitude, double radius, int numberOfPoints)
+        {
+            if(latitude < -90 || latitude > 90)
+            {
+                throw new ArgumentOutOfRangeException(nameof(latitude), "Latitude must be between -90 and 90 degrees");
+            }
+
+            if(longitude < -180 || longitude > 180)
+            {
+                throw new ArgumentOutOfRangeException(nameof(longitude), "Longitude must be between -180 and 180 degrees");
+            }
+
+            if(radius <= 0)
+            {
+                throw new ArgumentOutOfRangeException(nameof(radius), "Radius must be positive");
+            }
+
+            if(radius > max_radius_meters)
+            {
+                throw new ArgumentOutOfRangeException(nameof(radius), $"Radius must not exceed {max_radius_meters:F0} meters");
+            }
+
+            // Handle polar regions - near poles, create a simple circle in projected coordinates
+            if(Math.Abs(latitude) > 89.0)
+            {
+                throw new ArgumentOutOfRangeException(nameof(latitude), "Geodesic calculations are unstable near the poles (latitude > ±89°)");
+            }
+
+            if(numberOfPoints < 3)
+            {
+                numberOfPoints = 3;
+            }
+
+            // Pre-allocate with exact capacity (including closing point)
+            List<double[]> coordinates = new List<double[]>(numberOfPoints + 1);
+
+            // Pre-calculate common values to avoid repeated calculations
+            double angleIncrement = 2.0 * Math.PI / numberOfPoints;
+            double latRad = DegreesToRadians(latitude);
+            double lonRad = DegreesToRadians(longitude);
+            double angularDistance = radius / wgs84_equatorial_radius;
+
+            // Pre-calculate trigonometric values for the center point
+            double sinLatRad = Math.Sin(latRad);
+            double cosLatRad = Math.Cos(latRad);
+            double sinAngularDistance = Math.Sin(angularDistance);
+            double cosAngularDistance = Math.Cos(angularDistance);
+
+            for(int i = 0; i < numberOfPoints; i++)
+            {
+                double angle = angleIncrement * i;
+
+                // Calculate destination latitude 
+                double asinArg = sinLatRad * cosAngularDistance + cosLatRad * sinAngularDistance * Math.Cos(angle);
+                asinArg = Math.Max(-1.0, Math.Min(1.0, asinArg));
+                double lat2Rad = Math.Asin(asinArg);
+
+                // Calculate destination longitude
+                double lon2Rad = lonRad + Math.Atan2(
+                    Math.Sin(angle) * sinAngularDistance * cosLatRad,
+                    cosAngularDistance - sinLatRad * Math.Sin(lat2Rad)
+                );
+
+                // Convert to degrees and normalize longitude properly
+                double lat2 = RadiansToDegrees(lat2Rad);
+                double lon2 = RadiansToDegrees(lon2Rad);
+
+                // Normalize longitude to [-180, 180)
+                if(lon2 > 180.0)
+                {
+                    lon2 -= 360.0;
+                }
+                else if(lon2 <= -180.0)
+                {
+                    lon2 += 360.0;
+                }
+
+                coordinates.Add(new double[] { lon2, lat2 });
+            }
+
+            // Close the polygon by repeating the first point
+            coordinates.Add(new double[] { coordinates[0][0], coordinates[0][1] });
+
+            return coordinates;
+        }
+
+        /// <summary>
+        /// Converts degrees to radians.
+        /// </summary>
+        static double DegreesToRadians(double degrees) => degrees * Math.PI / 180.0;
+
+        /// <summary>
+        /// Converts radians to degrees.
+        /// </summary>
+        static double RadiansToDegrees(double radians) => radians * 180.0 / Math.PI;
+    }
+}

README.md

@@ -51,6 +51,30 @@ Cartesian bngCartesian = Transform.Etrs89ToOsgb36(cartesian);
 EastingNorthing bngEN = Convert.ToEastingNorthing(new Airy1830(), new BritishNationalGrid(), bngCartesian); 
 ```
 
+## Generate Polygon Around Point
+The `Polygon.GeneratePolygonAroundPoint` method creates a geodesically accurate polygon around a given latitude/longitude point.
+
+```csharp
+// Generate an 8-sided polygon with 1km radius around London
+var polygon = Polygon.GeneratePolygonAroundPoint(
+    longitude: -0.1276,    // London longitude
+    latitude: 51.5074,     // London latitude  
+    radius: 1000,          // 1000 meters (1km)
+    numberOfPoints: 8      // 8-sided polygon
+);
+
+// Result is List<double[]> where each array is [longitude, latitude]
+// The polygon is automatically closed (first point repeated at end)
+foreach (var coordinate in polygon)
+{
+    double longitude = coordinate[0];
+    double latitude = coordinate[1];
+    Console.WriteLine($"Point: {longitude:F6}, {latitude:F6}");
+}
+```
+> GeoJson is a useful website to test the generated polygon - [GeoJson.io](http://geojson.io/)
+
+
 ## Get OS Map reference
 The map references (Easting/Northing) used in Ordnance Survey maps are divided into 500km squares which are sub-divided into 100km squares. These squares are given a two letter code. The first letter represents the 500km square and the second represents the 100km square within it. A six digit map reference would look something like TL123456 where the first two characters represents the 100km square as indicated on the map with the first three digits of the six representing the easting and the last three digits representing the northing. Using this system means that a map reference is quoted as an easting/northing (in metres) from the square's origin. An EastingNorthing coordinate object, as returned from the transformation described above, can be converted to an OS map reference by using the Osgb36 class as follows:
 ```csharp

Tests/GeoUK.OSTN.XUnit/PolygonTests.cs

@@ -0,0 +1,322 @@
+using System;
+using System.Linq;
+using Xunit;
+using GeoUK;
+
+namespace GeoUK.OSTN.XUnit
+{
+    public class PolygonTests
+    {
+        #region Exception Tests - Invalid Input Parameters
+
+        [Theory]
+        [InlineData(-90.1)]
+        [InlineData(-91)]
+        [InlineData(-180)]
+        [InlineData(90.1)]
+        [InlineData(91)]
+        [InlineData(180)]
+        public void GeneratePolygonAroundPoint_InvalidLatitude_ThrowsArgumentOutOfRangeException(double invalidLatitude)
+        {
+            var exception = Assert.Throws<ArgumentOutOfRangeException>(() => 
+                Polygon.GeneratePolygonAroundPoint(0, invalidLatitude, 1000, 8));
+            Assert.Equal("latitude", exception.ParamName);
+            Assert.Contains("Latitude must be between -90 and 90 degrees", exception.Message);
+        }
+
+        [Theory]
+        [InlineData(-180.1)]
+        [InlineData(-181)]
+        [InlineData(180.1)]
+        [InlineData(181)]
+        [InlineData(360)]
+        public void GeneratePolygonAroundPoint_InvalidLongitude_ThrowsArgumentOutOfRangeException(double invalidLongitude)
+        {
+            var exception = Assert.Throws<ArgumentOutOfRangeException>(() => 
+                Polygon.GeneratePolygonAroundPoint(invalidLongitude, 0, 1000, 8));
+            Assert.Equal("longitude", exception.ParamName);
+            Assert.Contains("Longitude must be between -180 and 180 degrees", exception.Message);
+        }
+
+        [Theory]
+        [InlineData(0)]
+        [InlineData(-1)]
+        [InlineData(-1000)]
+        public void GeneratePolygonAroundPoint_InvalidRadius_ThrowsArgumentOutOfRangeException(double invalidRadius)
+        {
+            var exception = Assert.Throws<ArgumentOutOfRangeException>(() => 
+                Polygon.GeneratePolygonAroundPoint(0, 0, invalidRadius, 8));
+            Assert.Equal("radius", exception.ParamName);
+            Assert.Contains("Radius must be positive", exception.Message);
+        }
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_RadiusExceedsMaximum_ThrowsArgumentOutOfRangeException()
+        {
+            // max_radius_meters = Math.PI * 6378137.0 / 2.0 ≈ 10,018,754 meters
+            double maxRadius = Math.PI * 6378137.0 / 2.0;
+            double excessiveRadius = maxRadius + 1;
+
+            var exception = Assert.Throws<ArgumentOutOfRangeException>(() => 
+                Polygon.GeneratePolygonAroundPoint(0, 0, excessiveRadius, 8));
+            Assert.Equal("radius", exception.ParamName);
+            Assert.Contains("Radius must not exceed", exception.Message);
+        }
+
+        [Theory]
+        [InlineData(89.1)]
+        [InlineData(89.5)]
+        [InlineData(90)]
+        [InlineData(-89.1)]
+        [InlineData(-89.5)]
+        [InlineData(-90)]
+        public void GeneratePolygonAroundPoint_PolarRegions_ThrowsArgumentOutOfRangeException(double polarLatitude)
+        {
+            var exception = Assert.Throws<ArgumentOutOfRangeException>(() => 
+                Polygon.GeneratePolygonAroundPoint(0, polarLatitude, 1000, 8));
+            Assert.Equal("latitude", exception.ParamName);
+            Assert.Contains("Geodesic calculations are unstable near the poles", exception.Message);
+        }
+
+        #endregion
+
+        #region Boundary Value Tests
+
+        [Theory]
+        [InlineData(-180, -89)]
+        [InlineData(-180, 89)]
+        [InlineData(180, -89)]
+        [InlineData(180, 89)]
+        [InlineData(0, -89)]
+        [InlineData(0, 89)]
+        public void GeneratePolygonAroundPoint_BoundaryCoordinates_ReturnsValidPolygon(double longitude, double latitude)
+        {
+            var result = Polygon.GeneratePolygonAroundPoint(longitude, latitude, 1000, 8);
+            
+            Assert.NotNull(result);
+            Assert.Equal(9, result.Count); // 8 points + 1 closing point
+            Assert.True(result.All(coord => !double.IsNaN(coord[0]) && !double.IsNaN(coord[1])));
+        }
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_MinimumValidRadius_ReturnsValidPolygon()
+        {
+            double minRadius = double.Epsilon;
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, minRadius, 8);
+            
+            Assert.NotNull(result);
+            Assert.Equal(9, result.Count);
+            Assert.True(result.All(coord => !double.IsNaN(coord[0]) && !double.IsNaN(coord[1])));
+        }
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_MaximumValidRadius_ReturnsValidPolygon()
+        {
+            double maxRadius = Math.PI * 6378137.0 / 2.0; // Exactly at the limit
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, maxRadius, 8);
+            
+            Assert.NotNull(result);
+            Assert.Equal(9, result.Count);
+            Assert.True(result.All(coord => !double.IsNaN(coord[0]) && !double.IsNaN(coord[1])));
+        }
+
+        #endregion
+
+        #region NumberOfPoints Parameter Tests
+
+        [Theory]
+        [InlineData(0)]
+        [InlineData(1)]
+        [InlineData(2)]
+        public void GeneratePolygonAroundPoint_NumberOfPointsLessThanThree_DefaultsToThree(int numberOfPoints)
+        {
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, 1000, numberOfPoints);
+            
+            Assert.NotNull(result);
+            Assert.Equal(4, result.Count); // 3 points + 1 closing point
+        }
+
+        [Theory]
+        [InlineData(3, 4)]   // 3 points + 1 closing
+        [InlineData(4, 5)]   // 4 points + 1 closing
+        [InlineData(8, 9)]   // 8 points + 1 closing
+        [InlineData(12, 13)] // 12 points + 1 closing
+        [InlineData(360, 361)] // Large number of points
+        public void GeneratePolygonAroundPoint_VariousNumberOfPoints_ReturnsCorrectCount(int numberOfPoints, int expectedCount)
+        {
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, 1000, numberOfPoints);
+            
+            Assert.Equal(expectedCount, result.Count);
+        }
+
+        #endregion
+
+        #region Polygon Closure Tests
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_PolygonIsClosed_FirstAndLastPointsAreIdentical()
+        {
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, 1000, 8);
+            
+            var firstPoint = result.First();
+            var lastPoint = result.Last();
+            
+            Assert.Equal(firstPoint[0], lastPoint[0], 10); // longitude
+            Assert.Equal(firstPoint[1], lastPoint[1], 10); // latitude
+        }
+
+        #endregion
+
+        #region Coordinate Validation Tests
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_AllCoordinatesWithinValidRange_Success()
+        {
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, 1000, 8);
+            
+            foreach (var coord in result)
+            {
+                Assert.True(coord[0] >= -180.0 && coord[0] <= 180.0, $"Longitude {coord[0]} is out of range");
+                Assert.True(coord[1] >= -90.0 && coord[1] <= 90.0, $"Latitude {coord[1]} is out of range");
+            }
+        }
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_FloatingPointPrecision_DoesNotReturnNaN()
+        {
+            // Test edge cases that might cause floating-point precision issues
+            var result = Polygon.GeneratePolygonAroundPoint(179.999999, 89.0, 1000, 8);
+            Assert.True(result.All(coord => !double.IsNaN(coord[0]) && !double.IsNaN(coord[1])));
+        }
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_AntimeridianCrossing_HandlesCorrectly()
+        {
+            // Test longitude normalization across ±180°
+            var result = Polygon.GeneratePolygonAroundPoint(179.5, 0, 100000, 8);
+            Assert.True(result.All(coord => coord[0] >= -180.0 && coord[0] <= 180.0));
+        }
+
+        #endregion
+
+        #region Longitude Normalization Tests
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_LongitudeNormalization_HandlesAntimeridianCrossing()
+        {
+            // Test near antimeridian where polygon might cross ±180°
+            var result = Polygon.GeneratePolygonAroundPoint(179, 0, 200000, 8);
+            
+            Assert.True(result.All(coord => coord[0] >= -180.0 && coord[0] <= 180.0));
+            
+            // Should have points on both sides of antimeridian
+            bool hasPositiveLongitudes = result.Any(coord => coord[0] > 0);
+            bool hasNegativeLongitudes = result.Any(coord => coord[0] < 0);
+            Assert.True(hasPositiveLongitudes || hasNegativeLongitudes);
+        }
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_LongitudeNormalization_WesternHemisphere()
+        {
+            // Test near antimeridian from western side
+            var result = Polygon.GeneratePolygonAroundPoint(-179, 0, 200000, 8);
+            
+            Assert.True(result.All(coord => coord[0] >= -180.0 && coord[0] <= 180.0));
+        }
+
+        #endregion
+
+        #region Geometric Properties Tests
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_EquatorialRegion_ReturnsReasonableCoordinates()
+        {
+            double centerLon = 0;
+            double centerLat = 0;
+            double radius = 1000; // 1km
+            
+            var result = Polygon.GeneratePolygonAroundPoint(centerLon, centerLat, radius, 8);
+            
+            // All points should be roughly within expected distance from center
+            // For small distances at equator, degrees ≈ meters/111320
+            double expectedDegreeDelta = radius / 111320.0; // Rough conversion
+            
+            foreach (var coord in result.Take(result.Count - 1)) // Exclude closing point
+            {
+                double lonDiff = Math.Abs(coord[0] - centerLon);
+                double latDiff = Math.Abs(coord[1] - centerLat);
+                
+                Assert.True(lonDiff <= expectedDegreeDelta * 2, $"Longitude difference {lonDiff} too large");
+                Assert.True(latDiff <= expectedDegreeDelta * 2, $"Latitude difference {latDiff} too large");
+            }
+        }
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_LargeRadius_ReturnsValidCoordinates()
+        {
+            // Test with a large but valid radius
+            double largeRadius = 5000000; // 5000km, well within max limit
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, largeRadius, 8);
+            
+            Assert.NotNull(result);
+            Assert.Equal(9, result.Count);
+            Assert.True(result.All(coord => !double.IsNaN(coord[0]) && !double.IsNaN(coord[1])));
+        }
+
+        #endregion
+
+        #region Different Geographic Regions Tests
+
+        [Theory]
+        [InlineData(0, 0)]      // Equator, Prime Meridian
+        [InlineData(-120, 45)]  // North America
+        [InlineData(120, -30)]  // Australia
+        [InlineData(-60, -20)]  // South America
+        [InlineData(30, 60)]    // Northern Europe/Asia
+        [InlineData(100, 10)]   // Southeast Asia
+        public void GeneratePolygonAroundPoint_DifferentGeographicRegions_ReturnsValidPolygons(double longitude, double latitude)
+        {
+            var result = Polygon.GeneratePolygonAroundPoint(longitude, latitude, 10000, 8);
+            
+            Assert.NotNull(result);
+            Assert.Equal(9, result.Count);
+            Assert.True(result.All(coord => !double.IsNaN(coord[0]) && !double.IsNaN(coord[1])));
+            Assert.True(result.All(coord => coord[0] >= -180.0 && coord[0] <= 180.0));
+            Assert.True(result.All(coord => coord[1] >= -90.0 && coord[1] <= 90.0));
+        }
+
+        #endregion
+
+        #region Performance and Stress Tests
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_HighNumberOfPoints_PerformsReasonably()
+        {
+            // Test with a high number of points
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, 1000, 1000);
+            
+            Assert.Equal(1001, result.Count); // 1000 points + 1 closing
+            Assert.True(result.All(coord => !double.IsNaN(coord[0]) && !double.IsNaN(coord[1])));
+        }
+
+        #endregion
+
+        #region Array Structure Tests
+
+        [Fact]
+        public void GeneratePolygonAroundPoint_CoordinateArrayStructure_CorrectFormat()
+        {
+            var result = Polygon.GeneratePolygonAroundPoint(0, 0, 1000, 8);
+            
+            foreach (var coord in result)
+            {
+                Assert.NotNull(coord);
+                Assert.Equal(2, coord.Length);
+                Assert.True(double.IsFinite(coord[0]), "Longitude must be finite");
+                Assert.True(double.IsFinite(coord[1]), "Latitude must be finite");
+            }
+        }
+
+        #endregion
+    }
+}