Ytukyg

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 

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

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".

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

Find Distance

Use the Great Circle Distance Formula.

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

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

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;

    }

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 

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

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;

}