计算球面两点间距离实现Vincenty+Haversine

vincenty公式  精度很高能达到0.5毫米，但是很慢。

Haversine公式半正矢公式，比vincenty快，精度没有vincenty高，也长使用。

-------------------------------------------openlayers中实现的Vincenty----------------------------------------------------------

角度转弧度

/**
* Function: rad 
*
* Parameters:
* x - {Float}
*
* Returns:
* {Float}
*/
OpenLayers.Util.rad = function(x) {return x*Math.PI/180;};

 

弧度转角度

/**
* Function: deg
*
* Parameters:
* x - {Float}
*
* Returns:
* {Float}
*/
OpenLayers.Util.deg = function(x) {return x*180/Math.PI;};

 

 

a 长半轴

b短半轴

c 扁率

/**
* Property: VincentyConstants
* {Object} Constants for Vincenty functions.
*/
OpenLayers.Util.VincentyConstants = {
    a: 6378137,
    b: 6356752.3142,
    f: 1/298.257223563
};

 

	WGS-84	a = 6 378 137 m (±2 m)	b ≈ 6 356 752.314245 m	f ≈ 1 / 298.257223563
	GRS-80	a = 6 378 137 m	b ≈ 6 356 752.314140 m	f = 1 / 298.257222101
	Airy 1830	a = 6 377 563.396 m	b = 6 356 256.910 m	f ≈ 1 / 299.3249646
	Internat’l 1924	a = 6 378 388 m	b ≈ 6 356 911.946 m	f = 1 / 297
	Clarke mod.1880	a = 6 378 249.145 m	b ≈ 6 356 514.86955 m	f = 1 / 293.465
	GRS-67	a = 6 378 160 m	b ≈ 6 356 774.719 m	f = 1 / 298.247167

给定两个地理坐标（经纬度）返回km距离

/**
* APIFunction: distVincenty
* Given two objects representing points with geographic coordinates, this
*     calculates the distance between those points on the surface of an
*     ellipsoid.
*
* Parameters:
* p1 - {<OpenLayers.LonLat>} (or any object with both .lat, .lon properties)
* p2 - {<OpenLayers.LonLat>} (or any object with both .lat, .lon properties)
*
* Returns:
* {Float} The distance (in km) between the two input points as measured on an
*     ellipsoid.  Note that the input point objects must be in geographic
*     coordinates (decimal degrees) and the return distance is in kilometers.
*/
OpenLayers.Util.distVincenty = function(p1, p2) {
    var ct = OpenLayers.Util.VincentyConstants;
    var a = ct.a, b = ct.b, f = ct.f;

    var L = OpenLayers.Util.rad(p2.lon - p1.lon);
    var U1 = Math.atan((1-f) * Math.tan(OpenLayers.Util.rad(p1.lat)));
    var U2 = Math.atan((1-f) * Math.tan(OpenLayers.Util.rad(p2.lat)));
    var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
    var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
    var lambda = L, lambdaP = 2*Math.PI;
    var iterLimit = 20;
    while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0) {
        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 alpha = Math.asin(cosU1 * cosU2 * sinLambda / sinSigma);
        var cosSqAlpha = Math.cos(alpha) * Math.cos(alpha);
        var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
        var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
        lambdaP = lambda;
        lambda = L + (1-C) * f * Math.sin(alpha) *
        (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
    }
    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);
    var d = s.toFixed(3)/1000; // round to 1mm precision
    return d;
};

-------------------------------------------end----------------------------------------------------------

Haversine公式实现

 

function toRadians(degree) {

return degree * Math.PI / 180;

}

function distance(latitude1, longitude1, latitude2, longitude2) {

// R is the radius of the earth in kilometers

var R = 6371;

var deltaLatitude = toRadians(latitude2-latitude1);

var deltaLongitude = toRadians(longitude2-longitude1);

latitude1 =toRadians(latitude1);

latitude2 =toRadians(latitude2);

var a = Math.sin(deltaLatitude/2) *

Math.sin(deltaLatitude/2) +

Math.cos(latitude1) *

Math.cos(latitude2) *

Math.sin(deltaLongitude/2) *

Math.sin(deltaLongitude/2);

var c = 2 * Math.atan2(Math.sqrt(a),

Math.sqrt(1-a));

var d = R * c;

return d;

}