Skip to content

Instantly share code, notes, and snippets.

@bellbind
Created Apr 26, 2022
Embed
What would you like to do?
[javascript] geodesic direct/inverse by Vincenty formulae
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);
}
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};
}
}
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment