计算球面两点间距离实现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;

}