-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathVincentyFormula.js
66 lines (58 loc) · 3.14 KB
/
VincentyFormula.js
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2012 */
/* www.movable-type.co.uk/scripts/latlong-vincenty.html */
/* from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the */
/* Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975 */
/* http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
function toRad(value) {
/** Converts numeric degrees to radians */
return value * Math.PI / 180;
}
/**
* Calculates geodetic distance between two points specified by latitude/longitude using
* Vincenty inverse formula for ellipsoids
*
* @param {number} lat1 - latitude first point, in decimal degrees
* @param {number} lon1 - longitude first point, in decimal degrees
* @param {number} lat2 - latitude second point, in decimal degrees
* @param {number} lon2 - longitude second point, in decimal degrees
* @returns {number} - distance in metres between points
*/
function distVincenty(lat1, lon1, lat2, lon2) {
var a = 6378137, b = 6356752.314245, f = 1/298.257223563; // WGS-84 ellipsoid params
var L = toRad(lon2-lon1);
var U1 = Math.atan((1-f) * Math.tan(toRad(lat1)));
var U2 = Math.atan((1-f) * Math.tan(toRad(lat2)));
var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
var lambda = L, lambdaP, iterLimit = 100;
do {
var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
var sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
(cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
if (sinSigma==0) return 0; // co-incident points
var cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
var sigma = Math.atan2(sinSigma, cosSigma);
var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
var cosSqAlpha = 1 - sinAlpha*sinAlpha;
var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
if (isNaN(cos2SigmaM)) cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6)
var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
lambdaP = lambda;
lambda = L + (1-C) * f * sinAlpha *
(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
} while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);
if (iterLimit==0) {
return NaN; // formula failed to converge
}
var uSq = cosSqAlpha * (a*a - b*b) / (b*b);
var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
var s = b*A*(sigma-deltaSigma);
//s = s.toFixed(3); // round to 1mm precision
s = Math.round( s * 1000 ) / 1000; // round to 1mm precision
return s;
}