java 两个点球面距离_计算球面两点间距离实现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 - {} (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;

}