diff --git a/src/layers/awareness/rangeRing.test.mjs b/src/layers/awareness/rangeRing.test.mjs index 074d449fe..b17651090 100644 --- a/src/layers/awareness/rangeRing.test.mjs +++ b/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, diff --git a/src/layers/awareness/rendering.js b/src/layers/awareness/rendering.js index 8cee87d80..36d920850 100644 --- a/src/layers/awareness/rendering.js +++ b/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;