feat(foundry): add map camera spiral controls

2026-07-16 12:47:48 +03:00
parent 5e8c1cc3fc
commit e96401358b
11 changed files with 802 additions and 43 deletions
@@ -0,0 +1,180 @@
+const TWO_PI = Math.PI * 2;
+const WGS84_SEMI_MAJOR_METERS = 6_378_137;
+const WGS84_FLATTENING = 1 / 298.257_223_563;
+const WGS84_SEMI_MINOR_METERS = WGS84_SEMI_MAJOR_METERS * (1 - WGS84_FLATTENING);
+
+// This survey mode is deliberately regional. Radially projecting an
+// Archimedean spiral onto WGS84 stays within a negligible speed/tangent error
+// at this radius; a future continental mode needs an ellipsoid-arc
+// reparameterization instead of silently pretending the flat formula is exact.
+export const MAX_SPIRAL_RADIUS_METERS = 250_000;
+
+export type GeodeticRadians = {
+  longitude: number;
+  latitude: number;
+};
+
+export type SpiralSurfaceFrame = GeodeticRadians & {
+  angle: number;
+  radiusMeters: number;
+  radialBearing: number;
+  tangentHeading: number;
+};
+
+export function normalizeRadians(value: number) {
+  const normalized = ((value + Math.PI) % TWO_PI + TWO_PI) % TWO_PI - Math.PI;
+  return normalized === -Math.PI ? Math.PI : normalized;
+}
+
+export function spiralRadiusAtAngle(angle: number, pitchMetersPerTurn: number) {
+  const safeAngle = Math.max(0, finiteNumber(angle, "spiral_angle_invalid"));
+  const pitch = positiveNumber(pitchMetersPerTurn, "spiral_pitch_invalid");
+  return (pitch / TWO_PI) * safeAngle;
+}
+
+export function spiralArcLengthAtAngle(angle: number, pitchMetersPerTurn: number) {
+  const safeAngle = Math.max(0, finiteNumber(angle, "spiral_angle_invalid"));
+  const b = positiveNumber(pitchMetersPerTurn, "spiral_pitch_invalid") / TWO_PI;
+  return (b / 2) * (
+    safeAngle * Math.sqrt(1 + safeAngle * safeAngle)
+    + Math.asinh(safeAngle)
+  );
+}
+
+/**
+ * Inverts the Archimedean spiral arc-length function. This keeps the camera's
+ * configured surface speed constant instead of accelerating as the radius
+ * grows. Newton is monotone here; the bounded fallback protects the first
+ * samples near the origin from numerical overshoot.
+ */
+export function spiralAngleAtArcLength(surfaceDistanceMeters: number, pitchMetersPerTurn: number) {
+  const distance = Math.max(0, finiteNumber(surfaceDistanceMeters, "spiral_distance_invalid"));
+  const pitch = positiveNumber(pitchMetersPerTurn, "spiral_pitch_invalid");
+  if (distance === 0) return 0;
+  const b = pitch / TWO_PI;
+  let angle = distance <= b ? distance / b : Math.sqrt((2 * distance) / b);
+  for (let iteration = 0; iteration < 10; iteration += 1) {
+    const error = spiralArcLengthAtAngle(angle, pitch) - distance;
+    if (Math.abs(error) <= Math.max(1e-7, distance * 1e-12)) break;
+    const derivative = b * Math.sqrt(1 + angle * angle);
+    const next = angle - error / derivative;
+    angle = Number.isFinite(next) && next >= 0 ? next : angle / 2;
+  }
+  return angle;
+}
+
+/**
+ * Vincenty's WGS84 direct solution: move from a geodetic point along an
+ * initial bearing by an exact ellipsoid surface distance. Longitude/latitude
+ * are radians. The animation therefore remains stable at datelines and high
+ * latitudes and never treats degrees as a flat Cartesian plane.
+ */
+export function directGeodesicDestination(
+  origin: GeodeticRadians,
+  initialBearing: number,
+  distanceMeters: number,
+): GeodeticRadians {
+  const latitude1 = finiteNumber(origin.latitude, "spiral_origin_invalid");
+  const longitude1 = finiteNumber(origin.longitude, "spiral_origin_invalid");
+  const bearing = finiteNumber(initialBearing, "spiral_bearing_invalid");
+  const distance = Math.max(0, finiteNumber(distanceMeters, "spiral_distance_invalid"));
+  if (Math.abs(latitude1) > Math.PI / 2 + 1e-12) throw new Error("spiral_origin_invalid");
+  if (distance === 0) return { longitude: normalizeRadians(longitude1), latitude: latitude1 };
+
+  const sinBearing = Math.sin(bearing);
+  const cosBearing = Math.cos(bearing);
+  const tanReducedLatitude1 = (1 - WGS84_FLATTENING) * Math.tan(latitude1);
+  const cosReducedLatitude1 = 1 / Math.sqrt(1 + tanReducedLatitude1 * tanReducedLatitude1);
+  const sinReducedLatitude1 = tanReducedLatitude1 * cosReducedLatitude1;
+  const sigma1 = Math.atan2(tanReducedLatitude1, cosBearing);
+  const sinAlpha = cosReducedLatitude1 * sinBearing;
+  const cosSquaredAlpha = 1 - sinAlpha * sinAlpha;
+  const uSquared = cosSquaredAlpha
+    * (WGS84_SEMI_MAJOR_METERS ** 2 - WGS84_SEMI_MINOR_METERS ** 2)
+    / (WGS84_SEMI_MINOR_METERS ** 2);
+  const coefficientA = 1 + (uSquared / 16_384)
+    * (4096 + uSquared * (-768 + uSquared * (320 - 175 * uSquared)));
+  const coefficientB = (uSquared / 1024)
+    * (256 + uSquared * (-128 + uSquared * (74 - 47 * uSquared)));
+
+  let sigma = distance / (WGS84_SEMI_MINOR_METERS * coefficientA);
+  let previousSigma = Number.POSITIVE_INFINITY;
+  let sinSigma = 0;
+  let cosSigma = 1;
+  let cosTwoSigmaMiddle = 0;
+  for (let iteration = 0; iteration < 24 && Math.abs(sigma - previousSigma) > 1e-12; iteration += 1) {
+    cosTwoSigmaMiddle = Math.cos(2 * sigma1 + sigma);
+    sinSigma = Math.sin(sigma);
+    cosSigma = Math.cos(sigma);
+    const deltaSigma = coefficientB * sinSigma * (
+      cosTwoSigmaMiddle
+      + (coefficientB / 4) * (
+        cosSigma * (-1 + 2 * cosTwoSigmaMiddle ** 2)
+        - (coefficientB / 6) * cosTwoSigmaMiddle
+          * (-3 + 4 * sinSigma ** 2)
+          * (-3 + 4 * cosTwoSigmaMiddle ** 2)
+      )
+    );
+    previousSigma = sigma;
+    sigma = distance / (WGS84_SEMI_MINOR_METERS * coefficientA) + deltaSigma;
+  }
+
+  sinSigma = Math.sin(sigma);
+  cosSigma = Math.cos(sigma);
+  cosTwoSigmaMiddle = Math.cos(2 * sigma1 + sigma);
+  const temporary = sinReducedLatitude1 * sinSigma
+    - cosReducedLatitude1 * cosSigma * cosBearing;
+  const latitude2 = Math.atan2(
+    sinReducedLatitude1 * cosSigma + cosReducedLatitude1 * sinSigma * cosBearing,
+    (1 - WGS84_FLATTENING) * Math.sqrt(sinAlpha ** 2 + temporary ** 2),
+  );
+  const lambda = Math.atan2(
+    sinSigma * sinBearing,
+    cosReducedLatitude1 * cosSigma - sinReducedLatitude1 * sinSigma * cosBearing,
+  );
+  const coefficientC = (WGS84_FLATTENING / 16) * cosSquaredAlpha
+    * (4 + WGS84_FLATTENING * (4 - 3 * cosSquaredAlpha));
+  const longitudeDelta = lambda - (1 - coefficientC) * WGS84_FLATTENING * sinAlpha * (
+    sigma + coefficientC * sinSigma * (
+      cosTwoSigmaMiddle + coefficientC * cosSigma * (-1 + 2 * cosTwoSigmaMiddle ** 2)
+    )
+  );
+  return {
+    longitude: normalizeRadians(longitude1 + longitudeDelta),
+    latitude: latitude2,
+  };
+}
+
+export function spiralSurfaceFrame(
+  origin: GeodeticRadians,
+  initialHeading: number,
+  surfaceDistanceMeters: number,
+  pitchMetersPerTurn: number,
+): SpiralSurfaceFrame {
+  const angle = spiralAngleAtArcLength(surfaceDistanceMeters, pitchMetersPerTurn);
+  const radiusMeters = spiralRadiusAtAngle(angle, pitchMetersPerTurn);
+  if (radiusMeters > MAX_SPIRAL_RADIUS_METERS) throw new Error("spiral_extent_limit");
+  const radialBearing = normalizeRadians(initialHeading + angle);
+  const destination = directGeodesicDestination(origin, radialBearing, radiusMeters);
+  // In polar coordinates dr/dθ=b and r=bθ. The path tangent is therefore
+  // rotated atan2(r, dr/dθ)=atan(θ) from the outward radial direction.
+  const tangentHeading = normalizeRadians(radialBearing + Math.atan(angle));
+  return {
+    ...destination,
+    angle,
+    radiusMeters,
+    radialBearing,
+    tangentHeading,
+  };
+}
+
+function finiteNumber(value: number, errorCode: string) {
+  if (!Number.isFinite(value)) throw new Error(errorCode);
+  return value;
+}
+
+function positiveNumber(value: number, errorCode: string) {
+  const finite = finiteNumber(value, errorCode);
+  if (finite <= 0) throw new Error(errorCode);
+  return finite;
+}