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 - {} (or any object with both .lat, .lon properties)
* p2 - {} (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;
}
