worldmap-viewer: initial public history

1 files changed, 100 insertions, 0 deletions

diff --git a/viewer/js/geo.js b/viewer/js/geo.js
new file mode 100644
index 0000000..c6fc3e9
--- /dev/null
+++ b/viewer/js/geo.js
@@ -0,0 +1,100 @@
+// Sphere maths for World Maps. Scene convention (y up): x = cos φ cos λ, y = sin φ, z = −cos φ sin λ.
+import { PLANET } from './planet.js';
+export const R_KM = PLANET.radius_km;
+const D = Math.PI / 180;
+
+export function toVec(lat, lon) {
+  const la = lat * D, lo = lon * D, c = Math.cos(la);
+  return [c * Math.cos(lo), Math.sin(la), -c * Math.sin(lo)];
+}
+
+export function toLatLon(v) {
+  const n = Math.hypot(v[0], v[1], v[2]) || 1;
+  return { lat: Math.asin(Math.max(-1, Math.min(1, v[1] / n))) / D, lon: Math.atan2(-v[2], v[0]) / D };
+}
+
+const dot = (a, b) => a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+const cross = (a, b) => [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];
+
+export function angleRad(p, q) {
+  const a = toVec(p.lat, p.lon), b = toVec(q.lat, q.lon);
+  return Math.atan2(Math.hypot(...cross(a, b)), dot(a, b));
+}
+
+export const distanceKm = (p, q) => R_KM * angleRad(p, q);
+
+export function interpolate(p, q, stepDeg = 1) {
+  const a = toVec(p.lat, p.lon), b = toVec(q.lat, q.lon), w = angleRad(p, q);
+  const n = Math.max(1, Math.ceil(w / (stepDeg * D) - 1e-9));   // tolerate rounding: 3° at 1° steps = 3 pieces
+  const out = [];
+  for (let i = 0; i <= n; i++) {
+    if (w < 1e-12) { out.push({ lat: p.lat, lon: p.lon }); continue; }
+    const t = i / n, s1 = Math.sin((1 - t) * w) / Math.sin(w), s2 = Math.sin(t * w) / Math.sin(w);
+    out.push(toLatLon([a[0] * s1 + b[0] * s2, a[1] * s1 + b[1] * s2, a[2] * s1 + b[2] * s2]));
+  }
+  return out;
+}
+
+export function polygonAreaKm2(pts) {
+  if (pts.length < 3) return 0;
+  const v = pts.map(p => toVec(p.lat, p.lon));
+  let e = 0;
+  for (let i = 1; i + 1 < v.length; i++) {       // signed spherical excess of a triangle fan (Van Oosterom–Strackee)
+    const a = v[0], b = v[i], c = v[i + 1];
+    e += 2 * Math.atan2(dot(a, cross(b, c)), 1 + dot(a, b) + dot(b, c) + dot(c, a));
+  }
+  return Math.abs(e) * R_KM * R_KM;
+}
+
+// Equal Earth (Šavrič et al. 2018): same constants as map/mapgen/projections.py
+const A1 = 1.340264, A2 = -0.081106, A3 = 0.000893, A4 = 0.003796, M = Math.sqrt(3) / 2, P3 = Math.PI / 3;
+export const EE_XMAX = (2 * Math.sqrt(3) * Math.PI) / (3 * A1);
+export const EE_YMAX = A1 * P3 + A2 * P3 ** 3 + A3 * P3 ** 7 + A4 * P3 ** 9;
+export const PC_XMAX = Math.PI;
+export const PC_YMAX = Math.PI / 2;
+const eeD = t => A1 + 3 * A2 * t ** 2 + 7 * A3 * t ** 6 + 9 * A4 * t ** 8;
+
+export function project(view, lat, lon) {
+  if (view === 'plate_carree') return { x: lon * D, y: lat * D };
+  const t = Math.asin(M * Math.sin(lat * D));
+  return { x: (2 * Math.sqrt(3) * lon * D * Math.cos(t)) / (3 * eeD(t)), y: A1 * t + A2 * t ** 3 + A3 * t ** 7 + A4 * t ** 9 };
+}
+
+export function unproject(view, x, y) {
+  if (view === 'plate_carree') {
+    if (Math.abs(x) > PC_XMAX || Math.abs(y) > PC_YMAX) return null;
+    return { lat: y / D, lon: x / D };
+  }
+  if (Math.abs(y) > EE_YMAX) return null;
+  let t = y / A1;
+  for (let i = 0; i < 12; i++) t -= (A1 * t + A2 * t ** 3 + A3 * t ** 7 + A4 * t ** 9 - y) / eeD(t);
+  const s = Math.sin(t) / M;
+  if (Math.abs(s) > 1 + 1e-12) return null;
+  const lon = (3 * x * eeD(t)) / (2 * Math.sqrt(3) * Math.cos(t)) / D;
+  if (!Number.isFinite(lon) || Math.abs(lon) > 180 + 1e-9) return null;
+  return { lat: Math.asin(Math.max(-1, Math.min(1, s))) / D, lon: Math.max(-180, Math.min(180, lon)) };
+}
+
+export function sphereMesh(latSegs, lonSegs) {
+  const nv = (latSegs + 1) * (lonSegs + 1);
+  const positions = new Float32Array(nv * 3), uvs = new Float32Array(nv * 2);
+  for (let i = 0; i <= latSegs; i++) {
+    const lat = 90 - (180 * i) / latSegs;
+    for (let j = 0; j <= lonSegs; j++) {
+      const k = i * (lonSegs + 1) + j;
+      positions.set(toVec(lat, -180 + (360 * j) / lonSegs), 3 * k);
+      uvs[2 * k] = j / lonSegs;
+      uvs[2 * k + 1] = 1 - i / latSegs;
+    }
+  }
+  const indices = new Uint32Array(latSegs * lonSegs * 6);
+  let n = 0;
+  for (let i = 0; i < latSegs; i++) {
+    for (let j = 0; j < lonSegs; j++) {
+      const a = i * (lonSegs + 1) + j, b = a + lonSegs + 1, c = b + 1, d = a + 1;
+      indices.set([a, b, c, a, c, d], n);
+      n += 6;
+    }
+  }
+  return { positions, uvs, indices };
+}