如何求两个地点的经纬度距离？

Question

我有一组位置的纬度和经度。

• 如何找到从集合中的一个位置到另一个位置的距离？
• 有公式吗？

Answer 1

半正矢公式假定地球是球形的。然而，地球的形状更为复杂。扁球体模型会给出更好的结果。

如果需要这样的精度，最好使用文森特反演公式。
有关详细信息，请参阅 http://en.wikipedia.org/wiki/Vincenty's_formulae。使用它，您可以获得 0.5mm 的球体模型精度。

不存在完美的公式，因为地球的真实形状太复杂，无法用公式来表达。此外，地球的形状会因气候事件而发生变化（请参阅http://www.nasa.gov/centers/goddard/earthandsun/earthshape.html" rel="noreferrer">http://www. nasa.gov/centers/goddard/earthandsun/earthshape.html），并且由于地球自转也会随着时间而变化。

您还应该注意，上述方法没有考虑海拔高度，并假设海平面扁球体。

编辑 2010 年 7 月 10 日： 我发现在极少数情况下 Vincenty 反演公式不会收敛到声明的精度。更好的主意是使用 GeographicLib （请参阅 http://sourceforge.net/projects/geographiclib/）这也更准确。

Answer 2

这是一个： http://www.movable-type.co.uk/scripts/ latlong.html

使用半正矢公式：

R = earth’s radius (mean radius = 6,371km)
Δlat = lat2− lat1
Δlong = long2− long1
a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
d = R.c 

Answer 3

应用半正矢公式求出距离。请参阅下面的 C# 代码来查找 2 个坐标之间的距离。更好的是，如果您想要查找特定半径内的商店列表，您可以在 SQL 中应用 WHERE 子句或在 C# 中应用 LINQ 过滤器。

这里的公式以公里为单位，您必须更改相关数字，它将适用于英里。

例如：将 6371.392896 转换为英里。

    DECLARE @radiusInKm AS FLOAT
    DECLARE @lat2Compare AS FLOAT
    DECLARE @long2Compare AS FLOAT
    SET @radiusInKm = 5.000
    SET @lat2Compare = insert_your_lat_to_compare_here
    SET @long2Compare = insert_you_long_to_compare_here

    SELECT * FROM insert_your_table_here WITH(NOLOCK)
    WHERE (6371.392896*2*ATN2(SQRT((sin((radians(GeoLatitude - @lat2Compare)) / 2) * sin((radians(GeoLatitude - @lat2Compare)) / 2)) + (cos(radians(GeoLatitude)) * cos(radians(@lat2Compare)) * sin(radians(GeoLongitude - @long2Compare)/2) * sin(radians(GeoLongitude - @long2Compare)/2)))
    , SQRT(1-((sin((radians(GeoLatitude - @lat2Compare)) / 2) * sin((radians(GeoLatitude - @lat2Compare)) / 2)) + (cos(radians(GeoLatitude)) * cos(radians(@lat2Compare)) * sin(radians(GeoLongitude - @long2Compare)/2) * sin(radians(GeoLongitude - @long2Compare)/2)))
    ))) <= @radiusInKm

如果您想在 C# 中执行半正弦公式，

    double resultDistance = 0.0;
    double avgRadiusOfEarth = 6371.392896; //Radius of the earth differ, I'm taking the average.

    //Haversine formula
    //distance = R * 2 * aTan2 ( square root of A, square root of 1 - A )
    //                   where A = sinus squared (difference in latitude / 2) + (cosine of latitude 1 * cosine of latitude 2 * sinus squared (difference in longitude / 2))
    //                   and R = the circumference of the earth

    double differenceInLat = DegreeToRadian(currentLatitude - latitudeToCompare);
    double differenceInLong = DegreeToRadian(currentLongitude - longtitudeToCompare);
    double aInnerFormula = Math.Cos(DegreeToRadian(currentLatitude)) * Math.Cos(DegreeToRadian(latitudeToCompare)) * Math.Sin(differenceInLong / 2) * Math.Sin(differenceInLong / 2);
    double aFormula = (Math.Sin((differenceInLat) / 2) * Math.Sin((differenceInLat) / 2)) + (aInnerFormula);
    resultDistance = avgRadiusOfEarth * 2 * Math.Atan2(Math.Sqrt(aFormula), Math.Sqrt(1 - aFormula));

DegreesToRadian 是我自定义创建的一个函数，它是一个简单的 1 行“Math.PI * 角度” / 180.0

我的博客条目-SQL Haversine

Answer 4

您是否在寻找

Haversine 公式

半正矢公式是一个方程
在导航中很重要，给予
两点之间的大圆距离
从他们的球体上的点
经度和纬度。它是一个
更一般公式的特例
在球面三角学中，定律
半正矢，关联边和
球面“三角形”的角度。

Answer 5

看看这个..还有一个 javascript 示例。

查找距离

Answer 6

使用大圆距离公式。

Answer 7

这是一个通过给定 IP 查找长/纬度位置/附近位置的小提琴：

http://jsfiddle。 net/bassta/zrgd9qc3/2/

这是我用来计算直线距离的函数：

function distance(lat1, lng1, lat2, lng2) {
        var radlat1 = Math.PI * lat1 / 180;
        var radlat2 = Math.PI * lat2 / 180;
        var radlon1 = Math.PI * lng1 / 180;
        var radlon2 = Math.PI * lng2 / 180;
        var theta = lng1 - lng2;
        var radtheta = Math.PI * theta / 180;
        var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);
        dist = Math.acos(dist);
        dist = dist * 180 / Math.PI;
        dist = dist * 60 * 1.1515;

        //Get in in kilometers
        dist = dist * 1.609344;

        return dist;
    }

它返回以公里为单位的距离

Answer 8

如果您测量的距离小于（可能）1 度纬度/经度变化，正在寻找性能非常高的近似值，并且愿意接受比半正弦公式更多的不准确 ，考虑这两个替代方案：

（1）来自 计算距离：

a = pi/2 - lat1  
b = pi/2 - lat2  
c = sqrt( a^2 + b^2 - 2 * a * b * cos(lon2 - lon1) )   
d = R * c

(2)毕达哥拉斯定理根据纬度进行调整，如Ewan Todd 的著作中所示所以帖子：

d_ew = (long1 - long0) * cos(average(lat0, lat1))  
d_ns = (lat1 - lat0)  
d = sqrt(d_ew * d_ew + d_ns * d_ns)  

注释：
与 Ewan 的帖子相比，我在 cos( lat0 ) 中用 average(lat0, lat1) 替换了 lat0。

#2 对于值是度、弧度还是公里含糊不清；您还需要一些转换代码。请参阅本文底部的我的完整代码。

#1 的设计即使在极点附近也能正常工作，但如果您测量的距离的端点位于极点的“相反”两侧（经度相差超过 90 度？），建议使用半正弦值，即使距离很短。

我还没有彻底测量这些方法的错误，因此您应该为您的应用程序选取代表点，并将结果与​​一些高质量的库进行比较，以确定准确性是否可以接受。对于小于几公里的距离，我的直觉是这些距离与正确测量值的误差在 1% 以内。

---

获得高性能的另一种方法（如果适用）：

如果您有一大组静态点，位于经度/纬度的一到两度内，那么您将要计算少量动态（移动）点的距离，请考虑将静态点一次性转换为包含的UTM区域（或任何其他本地笛卡尔坐标系），然后在笛卡尔坐标系中进行所有数学运算。
笛卡尔 = 平坦的地球 = 毕达哥拉斯定理适用，因此距离 = sqrt(dx^2 + dy^2)。

那么，将少数移动点准确转换为UTM的成本就很容易承担。

---

#1（极地）注意事项：对于小于 0.1（？）米的距离可能会非常错误。即使使用双精度数学，以下坐标（其真实距离约为 0.005 米）在我的 Polar 算法实现中被指定为“零”：

输入：

    lon1Xdeg    16.6564465477996    double
    lat1Ydeg    57.7760262271983    double
    lon2Xdeg    16.6564466358281    double
    lat2Ydeg    57.776026248554 double

结果：

Oblate spheroid formula:  
    0.00575254911118364 double
Haversine:
    0.00573422966122257 double
Polar:
    0

这是由于两个因素 u和 v 完全相互抵消：

    u   0.632619944868587   double
    v   -0.632619944868587  double

在另一种情况下，当扁球体答案为 0.002887 m 时，它给出的距离为 0.067129 m。问题是 cos(lon2 - lon1) 与 1 太接近，因此 cos 函数返回的正是 1 。

除了测量亚米级距离之外，我还发现了迄今为止输入的有限小距离数据的最大误差（与扁球体公式相比）：

    maxHaversineErrorRatio  0.00350976281908381 double
    maxPolarErrorRatio  0.0510789996931342  double

其中“1”代表答案中 100% 的误差；例如，当它返回“0”时，即为“1”错误（从上面的“maxPolar”中排除）。因此“0.01”将是“100 分之一”或 1% 的误差。

将距离小于 2000 米的极坐标误差与半正弦误差进行比较，看看这个更简单的公式有多糟糕。到目前为止，我见过的最差情况是 Polar 为每 1000 份 51 份，而 Haversine 为每 1000 份 4 份。位于北纬58度左右。

---

现在实施“毕达哥拉斯与纬度调整”。

对于距离小于 1 的距离，它比 Polar 更加一致。 2000米
我原本以为极地问题只有在<<时才会出现。 1米，
但下面显示的结果相当令人不安。

当距离接近零时，毕达哥拉斯/纬度接近半正弦值。
例如，此测量 ~ 217 米：

    lon1Xdeg    16.6531667510102    double
    lat1Ydeg    57.7751705615804    double
    lon2Xdeg    16.6564468739869    double
    lat2Ydeg    57.7760263007586    double

    oblate      217.201200413731
    haversine   216.518428601051
    polar       226.128616011973
    pythag-cos  216.518428631907
    havErrRatio 0.00314349925958048
    polErrRatio 0.041102054598393
    pycErrRatio 0.00314349911751603

Polar 对这些输入的误差要大得多；要么是我的代码中存在一些错误，要么是我正在运行的 Cos 函数中存在错误，或者我必须建议不要使用 Polar，即使大多数 Polar 测量结果比这更接近。

OTOH，毕达哥拉斯，甚至使用 * cos(latitude) 调整时，误差的增加速度快于距离（距离越大，max_error/distance 的比率就会增加），因此您需要仔细考虑要测量的最大距离，以及可接受的错误。此外，不建议使用毕达哥拉斯比较两个几乎相等的距离来决定哪个更短，因为不同方向的误差不同（未显示证据）。

最坏情况测量，errorRatio = Abs(error) / distance（瑞典；最高 2000 m）：

    t_maxHaversineErrorRatio    0.00351012021578681 double
    t_maxPolarErrorRatio        66.0825360597085    double
    t_maxPythagoreanErrorRatio  0.00350976281416454 double

如前所述，极端极坐标误差适用于亚米距离，其中可能报告零而不是 6 cm，或者报告超过 0.5 m距离为 1 厘米（因此 t_maxPolarErrorRatio 中显示的是“66 x”最坏情况），但在较大距离下也有一些较差的结果。 [需要使用已知高度准确的余弦函数再次进行测试。]

在 Moto E4 上运行的 Xamarin.Android 中使用 C# 代码进行测量。

---

C#代码：

    // x=longitude, y= latitude. oblate spheroid formula. TODO: From where?
    public static double calculateDistanceDD_AED( double lon1Xdeg, double lat1Ydeg, double lon2Xdeg, double lat2Ydeg )
    {
        double c_dblEarthRadius = 6378.135; // km
        double c_dblFlattening = 1.0 / 298.257223563; // WGS84 inverse
                                                      // flattening
        // Q: Why "-" for longitudes??
        double p1x = -degreesToRadians( lon1Xdeg );
        double p1y = degreesToRadians( lat1Ydeg );
        double p2x = -degreesToRadians( lon2Xdeg );
        double p2y = degreesToRadians( lat2Ydeg );

        double F = (p1y + p2y) / 2;
        double G = (p1y - p2y) / 2;
        double L = (p1x - p2x) / 2;

        double sing = Math.Sin( G );
        double cosl = Math.Cos( L );
        double cosf = Math.Cos( F );
        double sinl = Math.Sin( L );
        double sinf = Math.Sin( F );
        double cosg = Math.Cos( G );

        double S = sing * sing * cosl * cosl + cosf * cosf * sinl * sinl;
        double C = cosg * cosg * cosl * cosl + sinf * sinf * sinl * sinl;
        double W = Math.Atan2( Math.Sqrt( S ), Math.Sqrt( C ) );
        if (W == 0.0)
            return 0.0;

        double R = Math.Sqrt( (S * C) ) / W;
        double H1 = (3 * R - 1.0) / (2.0 * C);
        double H2 = (3 * R + 1.0) / (2.0 * S);
        double D = 2 * W * c_dblEarthRadius;

        // Apply flattening factor
        D = D * (1.0 + c_dblFlattening * H1 * sinf * sinf * cosg * cosg - c_dblFlattening * H2 * cosf * cosf * sing * sing);

        // Transform to meters
        D = D * 1000.0;

        // tmstest
        if (true)
        {
            // Compare Haversine.
            double haversine = HaversineApproxDistanceGeo( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg );
            double error = haversine - D;
            double absError = Math.Abs( error );
            double errorRatio = absError / D;
            if (errorRatio > t_maxHaversineErrorRatio)
            {
                if (errorRatio > t_maxHaversineErrorRatio * 1.1)
                    Helper.test();
                t_maxHaversineErrorRatio = errorRatio;
            }

            // Compare Polar Coordinate Flat Earth. 
            double polarDistanceGeo = ApproxDistanceGeo_Polar( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg, D );
            double error2 = polarDistanceGeo - D;
            double absError2 = Math.Abs( error2 );
            double errorRatio2 = absError2 / D;
            if (errorRatio2 > t_maxPolarErrorRatio)
            {
                if (polarDistanceGeo > 0)
                {
                    if (errorRatio2 > t_maxPolarErrorRatio * 1.1)
                        Helper.test();
                    t_maxPolarErrorRatio = errorRatio2;
                }
                else
                    Helper.dubious();
            }

            // Compare Pythagorean Theorem with Latitude Adjustment. 
            double pythagoreanDistanceGeo = ApproxDistanceGeo_PythagoreanCosLatitude( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg, D );
            double error3 = pythagoreanDistanceGeo - D;
            double absError3 = Math.Abs( error3 );
            double errorRatio3 = absError3 / D;
            if (errorRatio3 > t_maxPythagoreanErrorRatio)
            {
                if (D < 2000)
                {
                    if (errorRatio3 > t_maxPythagoreanErrorRatio * 1.05)
                        Helper.test();
                    t_maxPythagoreanErrorRatio = errorRatio3;
                }
            }
        }


        return D;
    }

    // As a fraction of the distance.
    private static double t_maxHaversineErrorRatio, t_maxPolarErrorRatio, t_maxPythagoreanErrorRatio;


    // Average of equatorial and polar radii (meters).
    public const double EarthAvgRadius = 6371000;
    public const double EarthAvgCircumference = EarthAvgRadius * 2 * PI;
    // CAUTION: This is an average of great circles; won't be the actual distance of any longitude or latitude degree.
    public const double EarthAvgMeterPerGreatCircleDegree = EarthAvgCircumference / 360;

    // Haversine formula (assumes Earth is sphere).
    // "deg" = degrees.
    // Perhaps based on Haversine Formula in https://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
    public static double HaversineApproxDistanceGeo(double lon1Xdeg, double lat1Ydeg, double lon2Xdeg, double lat2Ydeg)
    {
        double lon1 = degreesToRadians( lon1Xdeg );
        double lat1 = degreesToRadians( lat1Ydeg );
        double lon2 = degreesToRadians( lon2Xdeg );
        double lat2 = degreesToRadians( lat2Ydeg );

        double dlon = lon2 - lon1;
        double dlat = lat2 - lat1;
        double sinDLat2 = Sin( dlat / 2 );
        double sinDLon2 = Sin( dlon / 2 );
        double a = sinDLat2 * sinDLat2 + Cos( lat1 ) * Cos( lat2 ) * sinDLon2 * sinDLon2;
        double c = 2 * Atan2( Sqrt( a ), Sqrt( 1 - a ) );
        double d = EarthAvgRadius * c;
        return d;
    }

    // From https://stackoverflow.com/a/19772119/199364
    // Based on Polar Coordinate Flat Earth in https://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
    public static double ApproxDistanceGeo_Polar( double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
    {
        double approxUnitDistSq = ApproxUnitDistSq_Polar(lon1deg, lat1deg, lon2deg, lat2deg, D);
        double c = Sqrt( approxUnitDistSq );
        return EarthAvgRadius * c;
    }

    // Might be useful to avoid taking Sqrt, when comparing to some threshold.
    // Threshold would have to be adjusted to match:  Power(threshold / EarthAvgRadius, 2)
    private static double ApproxUnitDistSq_Polar(double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
    {
        const double HalfPi = PI / 2; //1.5707963267949;

        double lon1 = degreesToRadians(lon1deg);
        double lat1 = degreesToRadians(lat1deg);
        double lon2 = degreesToRadians(lon2deg);
        double lat2 = degreesToRadians(lat2deg);

        double a = HalfPi - lat1;
        double b = HalfPi - lat2;
        double u = a * a + b * b;
        double dlon21 = lon2 - lon1;
        double cosDeltaLon = Cos( dlon21 );
        double v = -2 * a * b * cosDeltaLon;
        // TBD: Is "Abs" necessary?  That is, is "u + v" ever negative?
        //   (I think not; "v" looks like a secondary term. Though might be round-off issue near zero when a~=b.)
        double approxUnitDistSq = Abs(u + v);

        //if (approxUnitDistSq.nearlyEquals(0, 1E-16))
        //  Helper.dubious();
        //else if (D > 0)
        //{
        //  double dba = b - a;
        //  double unitD = D / EarthAvgRadius;
        //  double unitDSq = unitD * unitD;
        //  if (approxUnitDistSq > 2 * unitDSq)
        //      Helper.dubious();
        //  else if (approxUnitDistSq * 2 < unitDSq)
        //      Helper.dubious();
        //}

        return approxUnitDistSq;
    }

    // Pythagorean Theorem with Latitude Adjustment - from Ewan Todd - https://stackoverflow.com/a/1664836/199364
    // Refined by ToolmakerSteve - https://stackoverflow.com/a/53468745/199364
    public static double ApproxDistanceGeo_PythagoreanCosLatitude( double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
    {
        double approxDegreesSq = ApproxDegreesSq_PythagoreanCosLatitude( lon1deg, lat1deg, lon2deg, lat2deg );
        // approximate degrees on the great circle between the points.
        double d_degrees = Sqrt( approxDegreesSq );
        return d_degrees * EarthAvgMeterPerGreatCircleDegree;
    }

    public static double ApproxDegreesSq_PythagoreanCosLatitude( double lon1deg, double lat1deg, double lon2deg, double lat2deg )
    {
        double avgLatDeg = average( lat1deg , lat2deg );
        double avgLat = degreesToRadians( avgLatDeg );

        double d_ew = (lon2deg - lon1deg) * Cos( avgLat );
        double d_ns = (lat2deg - lat1deg);
        double approxDegreesSq = d_ew * d_ew + d_ns * d_ns;
        return approxDegreesSq;
    }

Answer 9

我已完成使用 SQL 查询

select *, (acos(sin(input_lat* 0.01745329)*sin(lattitude *0.01745329) + cos(input_lat *0.01745329)*cos(lattitude *0.01745329)*cos((input_long -longitude)*0.01745329))* 57.29577951 )* 69.16 As D  from table_name 

Answer 10

以下是包含前面答案中讨论的三个公式的模块（以 f90 编码）。您可以将此模块放在程序的顶部
（在 PROGRAM MAIN 之前）或单独编译并在编译期间包含模块目录。以下模块包含三个公式。前两个是基于地球是球形的假设的大圆距离。

module spherical_dists

contains

subroutine great_circle_distance(lon1,lat1,lon2,lat2,dist)
  !https://en.wikipedia.org/wiki/Great-circle_distance
  ! It takes lon, lats of two points on an assumed spherical earth and
  ! calculates the distance between them along the great circle connecting the two points
  implicit none
  real,intent(in)::lon1,lon2,lat1,lat2
  real,intent(out)::dist
  real,parameter::pi=3.141592,mean_earth_radius=6371.0088
  real::lonr1,lonr2,latr1,latr2
  real::delangl,dellon
  lonr1=lon1*(pi/180.);lonr2=lon2*(pi/180.)
  latr1=lat1*(pi/180.);latr2=lat2*(pi/180.)
  dellon=lonr2-lonr1
  delangl=acos(sin(latr1)*sin(latr2)+cos(latr1)*cos(latr2)*cos(dellon))
  dist=delangl*mean_earth_radius
end subroutine

subroutine haversine_formula(lon1,lat1,lon2,lat2,dist)
  ! https://en.wikipedia.org/wiki/Haversine_formula
  ! This is similar above but numerically better conditioned for small distances
  implicit none
  real,intent(in)::lon1,lon2,lat1,lat2
  !lon, lats of two points
  real,intent(out)::dist
  real,parameter::pi=3.141592,mean_earth_radius=6371.0088
  real::lonr1,lonr2,latr1,latr2
  real::delangl,dellon,dellat,a
  ! degrees are converted to radians
  lonr1=lon1*(pi/180.);lonr2=lon2*(pi/180.)
  latr1=lat1*(pi/180.);latr2=lat2*(pi/180.)
  dellon=lonr2-lonr1 ! These dels simplify the haversine formula
  dellat=latr2-latr1
  ! The actual haversine formula
  a=(sin(dellat/2))**2+cos(latr1)*cos(latr2)*(sin(dellon/2))**2
  delangl=2*asin(sqrt(a)) !2*asin(sqrt(a))
  dist=delangl*mean_earth_radius
end subroutine  

subroutine vincenty_formula(lon1,lat1,lon2,lat2,dist)
  !https://en.wikipedia.org/wiki/Vincenty%27s_formulae
  !It's a better approximation over previous two, since it considers earth to in oblate spheroid, which better approximates the shape of the earth
  implicit none
  real,intent(in)::lon1,lon2,lat1,lat2
  real,intent(out)::dist
  real,parameter::pi=3.141592,mean_earth_radius=6371.0088
  real::lonr1,lonr2,latr1,latr2
  real::delangl,dellon,nom,denom
  lonr1=lon1*(pi/180.);lonr2=lon2*(pi/180.)
  latr1=lat1*(pi/180.);latr2=lat2*(pi/180.)
  dellon=lonr2-lonr1
  nom=sqrt((cos(latr2)*sin(dellon))**2. + (cos(latr1)*sin(latr2)-sin(latr1)*cos(latr2)*cos(dellon))**2.)
  denom=sin(latr1)*sin(latr2)+cos(latr1)*cos(latr2)*cos(dellon)
  delangl=atan2(nom,denom)
  dist=delangl*mean_earth_radius
end subroutine

end module

Answer 11

在此页面上，您可以看到 Android Location 类中如何计算位置距离的完整代码和公式

android/location/Location.java

编辑：根据@Richard的提示，我把将链接函数的代码放入我的答案中，以避免链接无效：

private static void computeDistanceAndBearing(double lat1, double lon1,
    double lat2, double lon2, BearingDistanceCache results) {
    // Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
    // using the "Inverse Formula" (section 4)
    int MAXITERS = 20;
    // Convert lat/long to radians
    lat1 *= Math.PI / 180.0;
    lat2 *= Math.PI / 180.0;
    lon1 *= Math.PI / 180.0;
    lon2 *= Math.PI / 180.0;
    double a = 6378137.0; // WGS84 major axis
    double b = 6356752.3142; // WGS84 semi-major axis
    double f = (a - b) / a;
    double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);
    double L = lon2 - lon1;
    double A = 0.0;
    double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
    double U2 = Math.atan((1.0 - f) * Math.tan(lat2));
    double cosU1 = Math.cos(U1);
    double cosU2 = Math.cos(U2);
    double sinU1 = Math.sin(U1);
    double sinU2 = Math.sin(U2);
    double cosU1cosU2 = cosU1 * cosU2;
    double sinU1sinU2 = sinU1 * sinU2;
    double sigma = 0.0;
    double deltaSigma = 0.0;
    double cosSqAlpha = 0.0;
    double cos2SM = 0.0;
    double cosSigma = 0.0;
    double sinSigma = 0.0;
    double cosLambda = 0.0;
    double sinLambda = 0.0;
    double lambda = L; // initial guess
    for (int iter = 0; iter < MAXITERS; iter++) {
        double lambdaOrig = lambda;
        cosLambda = Math.cos(lambda);
        sinLambda = Math.sin(lambda);
        double t1 = cosU2 * sinLambda;
        double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
        double sinSqSigma = t1 * t1 + t2 * t2; // (14)
        sinSigma = Math.sqrt(sinSqSigma);
        cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
        sigma = Math.atan2(sinSigma, cosSigma); // (16)
        double sinAlpha = (sinSigma == 0) ? 0.0 :
            cosU1cosU2 * sinLambda / sinSigma; // (17)
        cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
        cos2SM = (cosSqAlpha == 0) ? 0.0 :
            cosSigma - 2.0 * sinU1sinU2 / cosSqAlpha; // (18)
        double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
        A = 1 + (uSquared / 16384.0) * // (3)
            (4096.0 + uSquared *
             (-768 + uSquared * (320.0 - 175.0 * uSquared)));
        double B = (uSquared / 1024.0) * // (4)
            (256.0 + uSquared *
             (-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
        double C = (f / 16.0) *
            cosSqAlpha *
            (4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
        double cos2SMSq = cos2SM * cos2SM;
        deltaSigma = B * sinSigma * // (6)
            (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))); // (11)
        double delta = (lambda - lambdaOrig) / lambda;
        if (Math.abs(delta) < 1.0e-12) {
            break;
        }
    }
    float distance = (float) (b * A * (sigma - deltaSigma));
    results.mDistance = distance;
    float initialBearing = (float) Math.atan2(cosU2 * sinLambda,
        cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);
    initialBearing *= 180.0 / Math.PI;
    results.mInitialBearing = initialBearing;
    float finalBearing = (float) Math.atan2(cosU1 * sinLambda,
            -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda);
    finalBearing *= 180.0 / Math.PI;
    results.mFinalBearing = finalBearing;
    results.mLat1 = lat1;
    results.mLat2 = lat2;
    results.mLon1 = lon1;
    results.mLon2 = lon2;
}

Answer 12

只需使用距离公式 Sqrt( (x2-x1)^2 + (y2-y1)^2 )