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.

Author: IeuanWalker

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