How to find distance from the latitude and longitu

问题:

I have a set of latitudes and longitudes of locations.

• How to find distance from one location in the set to another?
• Is there a formula ?

回答1:

The Haversine formula assumes a spherical earth. However, the shape of the earh is more complex. An oblate spheroid model will give better results.

If such accuracy is needed, you should better use Vincenty inverse formula. See http://en.wikipedia.org/wiki/Vincenty's_formulae for details. Using it, you can get a 0.5mm accuracy for the spheroid model.

There is no perfect formula, since the real shape of the earth is too complex to be expressed by a formula. Moreover, the shape of earth changes due to climate events (see http://www.nasa.gov/centers/goddard/earthandsun/earthshape.html), and also changes over time due to the rotation of the earth.

You should also note that the method above does not take altitudes into account, and assumes a sea-level oblate spheroid.

Edit 10-Jul-2010: I found out that there are rare situations for which Vincenty inverse formula does not converge to the declared accuracy. A better idea is to use GeographicLib (see http://sourceforge.net/projects/geographiclib/) which is also more accurate.

回答2:

Here's one: http://www.movable-type.co.uk/scripts/latlong.html

Using Haversine formula:

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 

回答3:

Apply the Haversine formula to find the distance. See the C# code below to find the distance between 2 coordinates. Better still if you want to say find a list of stores within a certain radius, you could apply a WHERE clause in SQL or a LINQ filter in C# to it.

The formula here is in kilometres, you will have to change the relevant numbers and it will work for miles.

E.g: Convert 6371.392896 to miles.

    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

If you would like to perform the Haversine formula in 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 is a function I custom created, its is a simple 1 liner of"Math.PI * angle / 180.0

My blog entry - SQL Haversine

回答4:

Are you looking for

Haversine formula

The haversine formula is an equation important in navigation, giving great-circle distances between two points on a sphere from their longitudes and latitudes. It is a special case of a more general formula in spherical trigonometry, the law of haversines, relating the sides and angles of spherical "triangles".

回答5:

Have a look at this.. has a javascript example as well.

Find Distance

回答6:

Use the Great Circle Distance Formula.

回答7:

This link has all the info you need, either on it or linked.

回答8:

here is a fiddle with finding locations / near locations to long/lat by given IP:

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

And here is the function I use to calculate the distance in straight line:

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;
    }

It returns the distance in Kilometers

回答9:

If you are measuring distances less than (perhaps) 1 degree lat/long change, are looking for a very high performance approximation, and are willing to accept more inaccuracy than Haversine formula, consider these two alternatives:

(1) "Polar Coordinate Flat-Earth Formula" from Computing Distances:

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

(2) Pythagorean theorem adjusted for latitude, as seen in Ewan Todd's SO post:

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

NOTES:
Compared to Ewan's post, I've substituted average(lat0, lat1) for lat0 inside of cos( lat0 ).

#2 is vague on whether values are degrees, radians, or kilometers; you will need some conversion code as well. See my complete code at bottom of this post.

#1 is designed to work well even near the poles, though if you are measuring a distance whose endpoints are on "opposite" sides of the pole (longitudes differ by more than 90 degrees?), Haversine is recommended instead, even for small distances.

I haven't thoroughly measured errors of these approaches, so you should take representative points for your application, and compare results to some high-quality library, to decide if the accuracies are acceptable. For distances less than a few kilometers my gut sense is that these are within 1% of correct measurement.

---

An alternative way to gain high performance (when applicable):

If you have a large set of static points, within one or two degrees of longitude/latitude, that you will then be calculating distances from a small number of dynamic (moving) points, consider converting your static points ONCE to the containing UTM zone (or to any other local Cartesian coordinate system), and then doing all your math in that Cartesian coordinate system.
Cartesian = flat earth = Pythagorean theorem applies, so distance = sqrt(dx^2 + dy^2).

Then the cost of accurately converting the few moving points to UTM is easily afforded.

---

CAVEAT for #1 (Polar): May be very wrong for distances less than 0.1 (?) meter. Even with double precision math, the following coordinates, whose true distance is about 0.005 meters, was given as "zero" by my implementation of Polar algorithm:

inputs:

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

results:

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

this was due to the two factors u and v exactly canceling each other:

    u   0.632619944868587   double
    v   -0.632619944868587  double

In another case, it gave a distance of 0.067129 m when the oblate spheroid answer was 0.002887 m. The problem was that cos(lon2 - lon1) was too close to 1, so cos function returned exactly 1.

Other than measuring sub-meter distances, the max errors (compared to an oblate spheroid formula) I found for the limited small-distance data I've fed in so far:

    maxHaversineErrorRatio  0.00350976281908381 double
    maxPolarErrorRatio  0.0510789996931342  double

where "1" would represent a 100% error in the answer; e.g. when it returned "0", that was an error of "1" (excluded from above "maxPolar"). So "0.01" would be an error of "1 part in 100" or 1%.

Comparing Polar error with Haversine error over distances less than 2000 meters to see how much worse this simpler formula is. So far, the worst I've seen is 51 parts per 1000 for Polar vs 4 parts per 1000 for Haversine. At about 58 degrees latitude.

---

Now implemented "Pythagorean with Latitude Adjustment".

It is MUCH more consistent than Polar for distances < 2000 m.
I originally thought the Polar problems were only when < 1 m,
but the result shown immediately below is quite troubling.

As distances approach zero, pythagorean/latitude approaches haversine. For example this measurement ~ 217 meters:

    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 has a much worse error with these inputs; either there is some mistake in my code, or in Cos function I am running on, or I have to recommend not using Polar, even though most Polar measurements were much closer than this.

OTOH, Pythagorean, even with * cos(latitude) adjustment, has error that increases more rapidly than distance (ratio of max_error/distance increases for larger distances), so you need to carefully consider the maximum distance you will measure, and the acceptable error. In addition, it is not advisable to COMPARE two nearly-equal distances using Pythagorean, to decide which is shorter, as the error is different in different DIRECTIONS (evidence not shown).

Worst case measurements, errorRatio = Abs(error) / distance (Sweden; up to 2000 m):

    t_maxHaversineErrorRatio    0.00351012021578681 double
    t_maxPolarErrorRatio        66.0825360597085    double
    t_maxPythagoreanErrorRatio  0.00350976281416454 double

As mentioned before, the extreme polar errors are for sub-meter distances, where it could report zero instead of 6 cm, or report over 0.5 m for a distance of 1 cm (hence the "66 x" worst case shown in t_maxPolarErrorRatio), but there are also some poor results at larger distances. [Needs to be tested again with a Cosine function that is known to be highly accurate.]

Measurements taken in C# code in Xamarin.Android running on a Moto E4.

---

C# code:

    // 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;
    }

回答10:

I am done using SQL query

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 

回答11:

Following is the module (coded in f90) containing three formulas discussed in the previous answers. You can either put this module at the top of your program (before PROGRAM MAIN) or compile it separately and include the module directory during compilation. The following module contains three formulas. First two are great-circle distances based on the assumption that earth is spherical.

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

回答12:

On this page you can see the whole code and formulas how distances of locations are calculated in Android Location class

android/location/Location.java

EDIT: According the hint from @Richard I put the code of the linked function into my answer, to avoid invalidated link:

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;
}

回答13:

just use the distance formula Sqrt( (x2-x1)^2 + (y2-y1)^2 )