---coding by yuangu(lifulinghan@aol.com)
--用于計算2個pgs之間的距離
function computeDistance(lat1, lon1,
lat2, lon2)
-- Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
-- using the "Inverse Formula" (section 4)
local MAXITERS = 20;
-- Convert lat/long to radians
lat1 = math.rad(lat1); -- lat1 * math.pi / 180.0;
lat2 = math.rad(lat2); --lat2 * math.pi / 180.0;
lon1 = math.rad(lon1); --lon1 * math.pi / 180.0;
lon2 = math.rad(lon2); --lon2 * math.pi / 180.0;
local a = 6378137.0; -- WGS84 major axis
local b = 6356752.3142; -- WGS84 semi-major axis
local f = (a - b) / a;
local aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);
local L = lon2 - lon1;
local A = 0.0;
local U1 = math.atan((1.0 - f) * math.tan(lat1));
local U2 = math.atan((1.0 - f) * math.tan(lat2));
local cosU1 = math.cos(U1);
local cosU2 = math.cos(U2);
local sinU1 = math.sin(U1);
local sinU2 = math.sin(U2);
local cosU1cosU2 = cosU1 * cosU2;
local sinU1sinU2 = sinU1 * sinU2;
local sigma = 0.0;
local deltaSigma = 0.0;
local cosSqAlpha = 0.0;
local cos2SM = 0.0;
local cosSigma = 0.0;
local sinSigma = 0.0;
local cosLambda = 0.0;
local sinLambda = 0.0;
local lambda = L; -- initial guess
for iter = 0, MAXITERS - 1 do
local lambdaOrig = lambda;
cosLambda = math.cos(lambda);
sinLambda = math.sin(lambda);
local t1 = cosU2 * sinLambda;
local t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
local sinSqSigma = t1 * t1 + t2 * t2;
sinSigma = math.sqrt(sinSqSigma);
cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda;
sigma = math.atan2(sinSigma, cosSigma);
local sinAlpha = 0.0
if (sinSigma ~= 0) then
sinAlpha = cosU1cosU2 * sinLambda / sinSigma
end
cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
cos2SM = 0.0
if cosSqAlpha ~= 0 then
cos2SM = cosSigma - 2.0 * sinU1sinU2 / cosSqAlpha;
end
local uSquared = cosSqAlpha * aSqMinusBSqOverBSq;
A = 1 + (uSquared / 16384.0) *
(4096.0 + uSquared *
(-768 + uSquared * (320.0 - 175.0 * uSquared)));
local B = (uSquared / 1024.0) *
(256.0 + uSquared *
(-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
local C = (f / 16.0) *
cosSqAlpha *
(4.0 + f * (4.0 - 3.0 * cosSqAlpha));
local cos2SMSq = cos2SM * cos2SM;
deltaSigma = B * sinSigma *
(cos2SM + (B / 4.0) *
(cosSigma * (-1.0 + 2.0 * cos2SMSq) -
(B / 6.0) * cos2SM *
(-3.0 + 4.0 * sinSigma * sinSigma) *
(-3.0 + 4.0 * cos2SMSq)));
lambda = L +
(1.0 - C) * f * sinAlpha *
(sigma + C * sinSigma *
(cos2SM + C * cosSigma *
(-1.0 + 2.0 * cos2SM * cos2SM)));
local delta = (lambda - lambdaOrig) / lambda;
if (math.abs(delta) < 1.0e-12) then
break;
end
end
return b * A * (sigma - deltaSigma);
end
--for test
print(computeDistance(34.8082342, 113.6125439, 34.8002478, 113.659779))
)
