fix(awareness): keep the range ring a full circle at the poles

At or within about 1e-7 degrees of a pole the latitude/longitude
destination formula cancels out, and the 128-point ring collapsed to
three distinct points joined by chords of about 500 km. Step along the
great circle in 3D from the centre's up, north and east unit vectors
instead. Results at ordinary latitudes and across the antimeridian are
unchanged.

src/layers/awareness/rangeRing.test.mjs

@@ -81,6 +81,41 @@ for (const [place, lat, lon] of [
   });
 }
 
+for (const [place, lat, lon] of [
+  ['the North Pole', 90, 0],
+  ['the South Pole', -90, 0],
+  ['1e-7° from the North Pole', 90 - 1e-7, 45],
+  ['1e-8° from the South Pole', -90 + 1e-8, -120],
+]) {
+  test(`the ring stays a full circle at ${place}`, () => {
+    const center = Cesium.Cartesian3.fromDegrees(lon, lat, 0);
+    const segments = AWARENESS_RANGE_RING_STYLE.segments;
+    const positions = awarenessRangeRingPositions(
+      center,
+      AWARENESS_RADIUS_M,
+      CONTEXT_RIM_HEIGHT_M,
+    );
+    assert.equal(positions.length, segments + 1);
+    assert.ok(
+      Cesium.Cartesian3.equalsEpsilon(positions[0], positions.at(-1), 0, 1e-6),
+      'the ring closes on itself',
+    );
+    const distinct = new Set(
+      positions
+        .slice(0, segments)
+        .map(({ x, y, z }) => [x, y, z].map((v) => v.toFixed(1)).join()),
+    );
+    assert.equal(distinct.size, segments, 'every vertex is its own point');
+    // Ordinary chords are 2π × 250 km / 128 ≈ 12.3 km long.
+    for (let i = 0; i < segments; i++) {
+      const chord = Cesium.Cartesian3.distance(positions[i], positions[i + 1]);
+      assert.ok(chord > 11_500 && chord < 13_000, `chord ${i}: ${chord} m`);
+    }
+    for (const distance of surfaceDistances(center, positions))
+      assert.ok(Math.abs(distance - AWARENESS_RADIUS_M) < 100, `${distance} m`);
+  });
+}
+
 test('the ring core and halo contrast on light and dark ground', () => {
   const core = Cesium.Color.fromCssColorString(
     AWARENESS_RANGE_RING_STYLE.color,

src/layers/awareness/rendering.js

@@ -40,9 +40,17 @@ export function awarenessRangeRingPositions(
   const e2 = 1 - (b * b) / (a * a);
   const sinLat = Math.sin(latitude);
   const cosLat = Math.cos(latitude);
+  const sinLon = Math.sin(longitude);
+  const cosLon = Math.cos(longitude);
   const w = 1 - e2 * sinLat * sinLat;
   const meridional = (a * (1 - e2)) / (w * Math.sqrt(w));
   const primeVertical = a / Math.sqrt(w);
+  // Step along the great circle in 3D from the centre's up, north and east
+  // unit vectors. Unlike the latitude/longitude destination formula, this
+  // stays well conditioned at and next to the poles.
+  const up = [cosLat * cosLon, cosLat * sinLon, sinLat];
+  const north = [-sinLat * cosLon, -sinLat * sinLon, cosLat];
+  const east = [-sinLon, cosLon, 0];
   const positions = [];
   for (let i = 0; i <= segments; i++) {
     const bearing = ((i % segments) / segments) * Cesium.Math.TWO_PI;
@@ -50,15 +58,13 @@ export function awarenessRangeRingPositions(
     const sinB = Math.sin(bearing);
     const angular =
       radiusM * ((cosB * cosB) / meridional + (sinB * sinB) / primeVertical);
-    const lat = Math.asin(
-      sinLat * Math.cos(angular) + cosLat * Math.sin(angular) * cosB,
+    const cosA = Math.cos(angular);
+    const sinA = Math.sin(angular);
+    const [x, y, z] = up.map(
+      (u, k) => u * cosA + (north[k] * cosB + east[k] * sinB) * sinA,
     );
-    const lon =
-      longitude +
-      Math.atan2(
-        sinB * Math.sin(angular) * cosLat,
-        Math.cos(angular) - sinLat * Math.sin(lat),
-      );
+    const lat = Math.atan2(z, Math.hypot(x, y));
+    const lon = Math.atan2(y, x);
     positions.push(Cesium.Cartesian3.fromRadians(lon, lat, heightM));
   }
   return positions;