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