Calculate two points distance in Google maps V3

// 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;
}