[javascript] geodesic direct/inverse by Vincenty formulae

example.js

import {inverse, direct} from "./vincenty.js";

{
  // example from https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2stf.html
  const latlon1 = [36.10377477777778, 140.08785502777778], latlon2 = [35.65502847222223, 139.74475044444443];
  console.log(inverse(latlon1, latlon2));
  const s =58643.80450350114, a1to2 = 211.99252761069826;
  console.log(direct(latlon1, a1to2, s));
  console.log(latlon2);
}

vincenty.js

export const wgs84 = Object.freeze({
  a: 6378137.06,
  b: 6356752.314245,
  f: 1 / 298.257223563,
});
export const grs80 = Object.freeze({
  a: 6378137,
  b: 6356582.314,
  f: 1 / 298.257222101,
});

// Vincenty formula: https://en.wikipedia.org/wiki/Vincenty%27s_formulae
export const inverse = ([lat1, lon1], [lat2, lon2], {a, b, f} = wgs84,  eps = Number.EPSILON) => {
  const {abs, atan, atan2, cos, sin, tan, PI, sqrt} = Math;
  const P1 = lat1 / 180 * PI, P2 = lat2 / 180 * PI, L = (lon2 - lon1) / 180 * PI;
  const U1 = atan((1 - f) * tan(P1)), U2 = atan((1 - f) * tan(P2));
  const sinU1 = sin(U1), cosU1 = cos(U1), sinU2 = sin(U2), cosU2 = cos(U2);
  const cosU1sinU2 = cosU1 * sinU2, sinU1cosU2 = sinU1 * cosU2, sinU1sinU2 = sinU1 * sinU2, cosU1cosU2 = cosU1 * cosU2;
  for (let l = L, i = 0; i < 1000; i++) {
    const sinL = sin(l), cosL = cos(l);
    const sinD = sqrt((cosU2 * sinL) ** 2 + (cosU1sinU2 - sinU1cosU2 * cosL) ** 2);
    const cosD = sinU1sinU2 + cosU1cosU2 * cosL;
    const D = atan2(sinD, cosD);
    const sinA = cosU1cosU2 * sinL / sinD;
    const cos2A = 1 - sinA ** 2;
    const cos2Dm = cosD - 2 * sinU1sinU2 / cos2A;
    const C = f / 16 * cos2A * (4 + f * (4 - 3 * cos2A));
    const cos22Dm = cos2Dm ** 2;
    const l_ = l;
    l = L + (1 - C) * f * sinA * (D + C * sinD * (cos2Dm + C * cosD * (-1 + 2 * cos22Dm)));
    if (abs(l - l_) < eps) {
      //console.log(i);
      const a2 = a ** 2, b2 = b ** 2, sin2D = sinD ** 2;
      const u2 = cos2A * (a2 - b2) / b2;
      const A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)));
      const B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)));
      const dD = B * sinD * (cos2Dm + B / 4 * (cosD * (-1 + 2 * cos22Dm) - B / 6 * cos2Dm * (-3 + 4 * sin2D) * (-3 + 4 * cos22Dm)));
      const s = b * A * (D - dD);
      const A1 = atan2(cosU2 * sinL, cosU1sinU2 - sinU1cosU2 * cosL);  // dir at latlon1
      const A2 = atan2(cosU1 * sinL, -sinU1cosU2 + cosU1sinU2 * cosL); // dir at latlong2, but dir from latlon1 to latlon2
      return {s, A1, A2, a1to2: (A1 / PI * 180 + 360) % 360, a2to1: (A2 / PI * 180 + 180) % 360};
    }
  }
};

export const direct = ([lat1, lon1], a1, s, {a, b, f} = wgs84, eps = Number.EPSILON) => {
  const {abs, atan, atan2, cos, sin, tan, PI, sqrt} = Math;
  const P1 = lat1 / 180 * PI, L1 = lon1 / 180 * PI, A1 = a1 / 180 * PI;
  const tanU1 = (1 - f) * tan(P1);
  const U1 = atan(tanU1);
  const cosU1 = cos(U1), sinU1 = sin(U1), cosA1 = cos(A1), sinA1 = sin(A1);
  const D1 = atan2(tanU1, cosA1);
  const sinA = cosU1 * sinA1;
  const cos2A = 1 - sinA ** 2;
  const u2 = cos2A * (a ** 2 - b ** 2) / (b ** 2);
  const A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)));
  const B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)));
  const D0 = s / b / A;
  for (let D = D0, i = 0; i < 1000; i++) {
    const sinD = sin(D), cosD = cos(D);
    const cos2Dm = cos(2 * D1 + D);
    const cos22Dm = cos2Dm ** 2, sin2D = sinD ** 2;
    const dD = B * sinD * (cos2Dm + B / 4 * (cosD * (-1 + 2 * cos22Dm) - B / 6 * cos2Dm * (-3 + 4 * sin2D) * (-3 + 4 * cos22Dm)));
    const D_ = D;
    D = D0 + dD;
    if (abs(D - D_) < eps) {
      //console.log(i);
      const P2 = atan2(sinU1 * cosD + cosU1 * sinD * cosA1, (1 - f) * sqrt(sinA ** 2 + (sinU1 * sinD - cosU1 * cosD * cosA1) ** 2));
      const l = atan2(sinD * sinA1, cosU1 * cosD - sinU1 * sinD * cosA1);
      const C = f / 16 * cos2A * (4 + f * (4 - 3 * cos2A));
      const L = l - (1 - C) * f * sinA * (D + C * sinD * (cos2Dm + C * cosD * (-1 + 2 * cos22Dm)));
      const L2 = L + L1;
      const A2 = atan2(sinA, -sinU1 * sinD + cosU1 * cosD * cosA1); // dir at latlong2, but dir from latlon1 to latlon2
      return {latlon2: [P2 / PI * 180, L2 / PI * 180], A2, a2to1: (A2 / PI * 180 + 180) % 360};
    }
  }
};