Run results for: Geo distance 2

Run details:

User agent: Mozilla/5.0 (Linux; Android 13) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/113.0.0.0 Mobile Safari/537.36

Browser: Chrome Mobile 113

Operating system: Android

Device Platform: Mobile

Date tested: one year ago

Test name	Executions per second
Vincenty reference	47457.5 Ops/sec
Vincenty	46601.5 Ops/sec
planar	2948912.8 Ops/sec
haversine	2301007.5 Ops/sec

Script Preparation code:

Benchmark.prototype.setup = function() {

    /** Converts numeric degrees to radians */
    if (typeof(Number.prototype.toRad) === "undefined") {
        Number.prototype.toRad = function() {
            return this * Math.PI / 180;
        }
    }

    /**
     * Calculates geodetic distance between two points specified by latitude/longitude using 
     * Vincenty inverse formula for ellipsoids
     *
     * @param   {Number} lat1, lon1: first point in decimal degrees
     * @param   {Number} lat2, lon2: second point in decimal degrees
     * @returns (Number} distance in metres between points
     */

    function distVincentyRef(lat1, lon1, lat2, lon2) {
        var a = 6378137,
            b = 6356752.314245,
            f = 1 / 298.257223563; // WGS-84 ellipsoid params
        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),
            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
        console.log(iterLimit);
        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;
    }

    function distVincenty(lat1, lon1, lat2, lon2) {
        const
        	a = 6378137,
            b = 6356752.314245,
            f = 1 / 298.257223563, // WGS-84 ellipsoid params
            degToRad = Math.PI / 180,
            f2 = f * f,
            f_4_4_16 = (f + f2) / 4,
            f_3_16 = f2 * 3 / 16,
            f1 = 1 - f,
            a2 = a * a,
            b2 = b * b,
            a2b2b2 = (a2 - b2) / b2;
        var L = (lon2 - lon1) * degToRad;
        var U1 = Math.atan(f1 * Math.tan(lat1 * degToRad));
        var U2 = Math.atan(f1 * Math.tan(lat2 * degToRad));
        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 cosU2sinLambda = cosU2 * sinLambda;
            var sinU1sinU2 = sinU1 * sinU2;
            var cosU1sinU2_sinU1cosU2cosLambda = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
            var sinSigma = Math.sqrt(cosU2sinLambda * cosU2sinLambda + cosU1sinU2_sinU1cosU2cosLambda * cosU1sinU2_sinU1cosU2cosLambda);
            if (sinSigma == 0) return 0; // co-incident points
            var cosSigma = sinU1sinU2 + cosU1 * cosU2 * cosLambda;
            var sigma = Math.atan2(sinSigma, cosSigma);
            var sinAlpha = cosU1 * cosU2sinLambda / sinSigma;
            var cosSqAlpha = 1 - sinAlpha * sinAlpha;
            var cos2SigmaM = cosSigma - 2 * sinU1sinU2 / cosSqAlpha;
            if (isNaN(cos2SigmaM)) cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6)
            var C = cosSqAlpha * (f_4_4_16 - f_3_16 * 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
        console.log(iterLimit);
        var uSq = cosSqAlpha * a2b2b2;
        var A = 1 + uSq * (0.25 + uSq * (-0.046875 + uSq * (0.01953125 - 0.01068115234 * uSq)));
        var B = uSq * (0.25 + uSq * (-0.125 + uSq * (0.072265625 - 0.0458984375 * 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;
    }
  
    function distPlanar(latDeg1, lonDeg1, latDeg2, lonDeg2) {
      const
        a = 6378137,
        b = 6356752.314245,
        a2 = a * a,
        b2 = b * b,
        ab4 = a2 * a2 * b2 * b2;

      var latRad1 = latDeg1 / 180.0 * Math.PI;
      var lonRad1 = lonDeg1 / 180.0 * Math.PI;
      var latRad2 = latDeg2 / 180.0 * Math.PI;
      var lonRad2 = lonDeg2 / 180.0 * Math.PI;

      var latRad = (latRad1 + latRad2) * 0.5;
      var lonRad = (lonRad1 + lonRad2) * 0.5;

      var sinLat = Math.sin(latRad);
      var sinLat2 = sinLat * sinLat;
      var a2cos2 = a2 * (1 - sinLat2);
      var b2sin2 = b2 * sinLat2;

      var s = a2cos2 + b2sin2;
      var latScale2 = ab4 / (s * s * s);
      var lonScale2 = a2 / (1 + b2sin2 / a2cos2);

      var dLat = latRad1 - latRad2;
      var dLon = lonRad1 - lonRad2;
      return Math.sqrt(latScale2 * dLat * dLat + lonScale2 * dLon * dLon).toFixed(3);
    }

    // Distance in kilometers between two points using the Haversine algo.
    function haversine(lat1, lon1, lat2, lon2) {
        var R = 6371;
        var dLat = (lat2 - lat1).toRad();
        var dLong = (lon2 - lon1).toRad();

        var a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(lat1.toRad()) * Math.cos(lat2.toRad()) * Math.sin(dLong / 2) * Math.sin(dLong / 2);
        var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
        var d = R * c;

        return Math.round(d);
    }
};

Tests:

Vincenty reference
distVincentyRef(-31.965379, 115.822077, -33.85789, 151.214806);
Vincenty
distVincenty(-31.965379, 115.822077, -33.85789, 151.214806);
planar
distPlanar(-31.965379, 115.822077, -33.85789, 151.214806);
haversine
haversine(-31.965379, 115.822077, -33.85789, 151.214806);

Run results for: Geo distance 2

Vincenty reference

Vincenty

planar

haversine