// Thanks to Movable Type Scripts
google.maps.LatLng.prototype.distanceFrom = function(latLng) {
if (!latLng) {
return NaN; // not a valid LatLng
}
// WGS-84 ellipsoid params
var a = 6378137; // m
var b = 6356752.314245; m
var f = 1/298.257223563;
var lat1 = this.lat();
var lon1 = this.lng();
var lat2 = latLng.lat();
var lon2 = latLng.lng();
var L = (lon2-lon1).toRad();
var U1 = Math.atan((1-f) * Math.tan(lat1.toRad()));
var U2 = Math.atan((1-f) * Math.tan(lat2.toRad()));
var sinU1 = Math.sin(U1);
var cosU1 = Math.cos(U1);
var sinU2 = Math.sin(U2);
var cosU2 = Math.cos(U2);
var lambda = L;
var lambdaP;
var iterLimit = 100;
do {
var sinLambda = Math.sin(lambda);
var 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.000; // 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
return s;
}