// 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 }; }