// @wuys/map · math/geodesy.js
// Earth-sphere geometry helpers — used by data/centroids, layer/arcs, etc.

const DEG = Math.PI / 180;

/** Great-circle angular distance (radians) between two (lon, lat) in degrees. */
export function angularDistance(a, b) {
  const lat1 = a.lat * DEG, lat2 = b.lat * DEG;
  const dLon = (b.lon - a.lon) * DEG;
  const cos = Math.sin(lat1) * Math.sin(lat2) + Math.cos(lat1) * Math.cos(lat2) * Math.cos(dLon);
  return Math.acos(Math.max(-1, Math.min(1, cos)));
}

/** Bearing (radians) from a to b. */
export function bearing(a, b) {
  const lat1 = a.lat * DEG, lat2 = b.lat * DEG;
  const dLon = (b.lon - a.lon) * DEG;
  const y = Math.sin(dLon) * Math.cos(lat2);
  const x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(dLon);
  return Math.atan2(y, x);
}

/**
 * Polygon shoelace area + centroid in the polygon's plane (here used on lon/lat
 * which is fine for our purposes — small relative-area errors at high latitudes,
 * acceptable since we use this only to pick the "biggest" sub-polygon).
 */
export function shoelaceCentroid(points) {
  let twoArea = 0, cx = 0, cy = 0;
  const n = points.length;
  for (let i = 0; i < n; i++) {
    const [x0, y0] = points[i];
    const [x1, y1] = points[(i + 1) % n];
    const cross = x0 * y1 - x1 * y0;
    twoArea += cross;
    cx += (x0 + x1) * cross;
    cy += (y0 + y1) * cross;
  }
  if (Math.abs(twoArea) < 1e-9) return null;
  return {
    x: cx / (3 * twoArea),
    y: cy / (3 * twoArea),
    area: Math.abs(twoArea) / 2,
  };
}
