geo.bounding-box

import { type BoundingBox } from "./geo_bounding_box_types.ts";

/**
 * The latitude/longitude box that contains every point within a distance of a
 * centre, on a sphere of the IUGG mean Earth radius, by the method of
 * J. P. Matuschek, "Finding Points Within a Distance of a Latitude/Longitude
 * Using Bounding Coordinates" (http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates).
 *
 * The trigonometry is written out below instead of calling Math.sin and
 * friends, because the platform maths library may differ in the last bit
 * between JavaScript, Python and Rust. Only correctly-rounded IEEE 754
 * operations (+ - * / sqrt floor) are used, in the same order in all three,
 * so every language computes the same doubles before rounding.
 */

// IUGG mean radius R1 of the GRS 80 ellipsoid (Moritz, Journal of Geodesy 74 (2000)).
const EARTH_RADIUS_METRES = 6371008.8;
// Ten millionths of a degree, about a centimetre.
const SCALE = 10000000;
// Values within a millionth of a grid step of a 7-decimal value are taken to
// be that value, so a centre of 51.5074 stays 51.5074 instead of becoming
// 51.5073999 because its double sits a hair below. That is about 1e-13
// degrees, far below anything the box is used for.
const SNAP = 0.000001;

const PI = 3.141592653589793;
const HALF_PI = PI / 2;
const DEGREES = PI / 180;

// Taylor coefficients, each one correctly-rounded division.
const S3 = -1 / 6, S5 = 1 / 120, S7 = -1 / 5040, S9 = 1 / 362880, S11 = -1 / 39916800;
const S13 = 1 / 6227020800, S15 = -1 / 1307674368000, S17 = 1 / 355687428096000;
const C2 = -1 / 2, C4 = 1 / 24, C6 = -1 / 720, C8 = 1 / 40320, C10 = -1 / 3628800;
const C12 = 1 / 479001600, C14 = -1 / 87178291200, C16 = 1 / 20922789888000, C18 = -1 / 6402373705728000;
const A3 = -1 / 3, A5 = 1 / 5, A7 = -1 / 7, A9 = 1 / 9, A11 = -1 / 11;
const A13 = 1 / 13, A15 = -1 / 15, A17 = 1 / 17, A19 = -1 / 19, A21 = 1 / 21;

// sin and cos on |r| <= pi/4, where the series converge to well below an ulp.
function sinSmall(r: number): number {
  const s = r * r;
  let p = S17;
  p = S15 + s * p;
  p = S13 + s * p;
  p = S11 + s * p;
  p = S9 + s * p;
  p = S7 + s * p;
  p = S5 + s * p;
  p = S3 + s * p;
  return r + r * s * p;
}

function cosSmall(r: number): number {
  const s = r * r;
  let p = C18;
  p = C16 + s * p;
  p = C14 + s * p;
  p = C12 + s * p;
  p = C10 + s * p;
  p = C8 + s * p;
  p = C6 + s * p;
  p = C4 + s * p;
  p = C2 + s * p;
  return 1 + s * p;
}

// Reduce to a quarter turn k and a remainder |r| <= pi/4.
function sinCos(x: number): [number, number] {
  const k = Math.floor(x / HALF_PI + 0.5);
  const r = x - k * HALF_PI;
  const sr = sinSmall(r);
  const cr = cosSmall(r);
  const quadrant = ((k % 4) + 4) % 4;
  if (quadrant === 0) return [sr, cr];
  if (quadrant === 1) return [cr, -sr];
  if (quadrant === 2) return [-sr, -cr];
  return [-cr, sr];
}

// atan for 0 <= u <= 1: halve the angle three times (tan(t/2) = u / (1 + sqrt(1 + u^2)))
// so the series runs on |u| <= tan(pi/32), then scale back by 8.
function atanUnit(u: number): number {
  let v = u;
  v = v / (1 + Math.sqrt(1 + v * v));
  v = v / (1 + Math.sqrt(1 + v * v));
  v = v / (1 + Math.sqrt(1 + v * v));
  const s = v * v;
  let p = A21;
  p = A19 + s * p;
  p = A17 + s * p;
  p = A15 + s * p;
  p = A13 + s * p;
  p = A11 + s * p;
  p = A9 + s * p;
  p = A7 + s * p;
  p = A5 + s * p;
  p = A3 + s * p;
  return 8 * (v + v * s * p);
}

// asin for 0 <= x, via atan; (1 - x)(1 + x) keeps precision near 1.
function asinPositive(x: number): number {
  if (x >= 1) return HALF_PI;
  const t = x / Math.sqrt((1 - x) * (1 + x));
  return t > 1 ? HALF_PI - atanUnit(1 / t) : atanUnit(t);
}

// Outward rounding: a lower bound goes down and an upper bound goes up, so
// rounding can only grow the box, never cut the circle.
function roundDown(x: number): number {
  return Math.floor(x * SCALE + SNAP) / SCALE + 0;
}

function roundUp(x: number): number {
  return Math.ceil(x * SCALE - SNAP) / SCALE + 0;
}

function checkNumber(value: number, what: string): void {
  if (typeof value !== "number" || !Number.isFinite(value)) {
    throw new TypeError(`${what} must be a finite number`);
  }
}

export function boundingBox(lat: number, lng: number, distanceMetres: number): BoundingBox {
  checkNumber(lat, "latitude");
  checkNumber(lng, "longitude");
  checkNumber(distanceMetres, "distance");
  if (lat < -90 || lat > 90) throw new RangeError("latitude must be between -90 and 90 degrees");
  if (lng < -180 || lng > 180) throw new RangeError("longitude must be between -180 and 180 degrees");
  if (distanceMetres < 0) throw new RangeError("distance must be 0 or more metres");

  const r = distanceMetres / EARTH_RADIUS_METRES;
  const rDegrees = r / DEGREES;
  let minLat = lat - rDegrees;
  let maxLat = lat + rDegrees;
  let minLng: number;
  let maxLng: number;

  if (minLat > -90 && maxLat < 90) {
    const [sinR] = sinCos(r);
    const [, cosLat] = sinCos(lat * DEGREES);
    const dLng = asinPositive(sinR / cosLat) / DEGREES;
    minLng = lng - dLng;
    maxLng = lng + dLng;
    // Past the antimeridian the bound comes round the other side, leaving
    // minLng > maxLng: the box is the two strips either side of 180.
    if (minLng < -180) minLng += 360;
    if (maxLng > 180) maxLng -= 360;
  } else {
    // The circle reaches a pole, so it covers every longitude.
    if (minLat < -90) minLat = -90;
    if (maxLat > 90) maxLat = 90;
    minLng = -180;
    maxLng = 180;
  }

  return { minLat: roundDown(minLat), minLng: roundDown(minLng), maxLat: roundUp(maxLat), maxLng: roundUp(maxLng) };
}