This link might be helpful to you, as it details the use of the Haversine formula to calculate the distance.
Excerpt:
This script [in Javascript] calculates great-circle distances between the two points –
that is, the shortest distance over the earth’s surface – using the
‘Haversine’ formula.
function getDistanceFromLatLonInKm(lat1,lon1,lat2,lon2) {
var R = 6371; // Radius of the earth in km
var dLat = deg2rad(lat2-lat1); // deg2rad below
var dLon = deg2rad(lon2-lon1);
var a =
Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(deg2rad(lat1)) * Math.cos(deg2rad(lat2)) *
Math.sin(dLon/2) * Math.sin(dLon/2)
;
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c; // Distance in km
return d;
}
function deg2rad(deg) {
return deg * (Math.PI/180)
}
I implemented a WGS84 distance function using the average of the start and end altitude as the constant altitude. If you are certain that there will be relatively little altitude variation along your path this works acceptably well (error is relative to the altitude difference of your two LLA points).
Here's my code (C#):
/// <summary>
/// Gets the geodesic distance between two pathpoints in the current mode's coordinate system
/// </summary>
/// <param name="point1">First point</param>
/// <param name="point2">Second point</param>
/// <param name="mode">Coordinate mode that both points are in</param>
/// <returns>Distance between the two points in the current coordinate mode</returns>
public static double GetGeodesicDistance(PathPoint point1, PathPoint point2, CoordMode mode) {
// calculate proper geodesics for LLA paths
if (mode == CoordMode.LLA) {
// meeus approximation
double f = (point1.Y + point2.Y) / 2 * LatLonAltTransformer.DEGTORAD;
double g = (point1.Y - point2.Y) / 2 * LatLonAltTransformer.DEGTORAD;
double l = (point1.X - point2.X) / 2 * LatLonAltTransformer.DEGTORAD;
double sinG = Math.Sin(g);
double sinL = Math.Sin(l);
double sinF = Math.Sin(f);
double s, c, w, r, d, h1, h2;
// not perfect but use the average altitude
double a = (LatLonAltTransformer.A + point1.Z + LatLonAltTransformer.A + point2.Z) / 2.0;
sinG *= sinG;
sinL *= sinL;
sinF *= sinF;
s = sinG * (1 - sinL) + (1 - sinF) * sinL;
c = (1 - sinG) * (1 - sinL) + sinF * sinL;
w = Math.Atan(Math.Sqrt(s / c));
r = Math.Sqrt(s * c) / w;
d = 2 * w * a;
h1 = (3 * r - 1) / 2 / c;
h2 = (3 * r + 1) / 2 / s;
return d * (1 + (1 / LatLonAltTransformer.RF) * (h1 * sinF * (1 - sinG) - h2 * (1 - sinF) * sinG));
}
PathPoint diff = new PathPoint(point2.X - point1.X, point2.Y - point1.Y, point2.Z - point1.Z, 0);
return Math.Sqrt(diff.X * diff.X + diff.Y * diff.Y + diff.Z * diff.Z);
}
In practice we've found that the altitude difference rarely makes a large difference, our paths are typically 1-2km long with altitude varying on the order of 100m and we see about ~5m change on average versus using the WGS84 ellipsoid unmodified.
Edit:
To add to this, if you do expect large altitude changes, you can convert your WGS84 coordinates to ECEF (earth centered earth fixed) and evaluate straight-line paths as shown at the bottom of my function. Converting a point to ECEF is simple to do:
/// <summary>
/// Converts a point in the format (Lon, Lat, Alt) to ECEF
/// </summary>
/// <param name="point">Point as (Lon, Lat, Alt)</param>
/// <returns>Point in ECEF</returns>
public static PathPoint WGS84ToECEF(PathPoint point) {
PathPoint outPoint = new PathPoint(0);
double lat = point.Y * DEGTORAD;
double lon = point.X * DEGTORAD;
double e2 = 1.0 / RF * (2.0 - 1.0 / RF);
double sinLat = Math.Sin(lat), cosLat = Math.Cos(lat);
double chi = A / Math.Sqrt(1 - e2 * sinLat * sinLat);
outPoint.X = (chi + point.Z) * cosLat * Math.Cos(lon);
outPoint.Y = (chi + point.Z) * cosLat * Math.Sin(lon);
outPoint.Z = (chi * (1 - e2) + point.Z) * sinLat;
return outPoint;
}
Edit 2:
I was asked about some of the other variables in my code:
// RF is the eccentricity of the WGS84 ellipsoid
public const double RF = 298.257223563;
// A is the radius of the earth in meters
public const double A = 6378137.0;
LatLonAltTransformer
is a class I used to convert from LatLonAlt coordinates to ECEF coordinates and defines the constants above.
Best Answer
Haversine and Vincenty are two algorithms for solving different problems. Haversine computes the great circle distance on a sphere while Vincenty computes the shortest (geodesic) distance on the surface of an ellipsoid of revolution. So the answer to your question can be broken into 2 parts:
For terrestrial applications, an ellipsoid of revolution is a reasonable approximation to "mean sea level"; the error is ± 100 m. The flattening of this ellipsoid is small, about 1/300, and so can be approximated by a sphere (of equal volume, for example).
Great circle distances differ from geodesic distances by up to 0.5%. In some applications, e.g., what's the distance from the Cape to Cairo?, this error can be neglected. In other applications, e.g., determining maritime boundaries, it is far too large (it's 5 m over a distance of 1 km). In general, you're safer using the geodesic distance.
If you're interested is distance traveled (by car, boat, or plane), there are lots of constraints on the path taken and neither the great circle or geodesic distance, which measure the length of shortest paths on an ideal surface, would be appropriate.
On the question of whether the algorithms are accurate:
Haversine is accurate to round-off unless the points are nearly antipodal. Better formulas are given in the Wikipedia article on great-circle distances.
Vincenty is usually accurate to about 0.1 mm. However if the points are nearly antipodal, the algorithm fails to converge and the error is much larger. I give a better algorithm for solving the geodesic problem in Algorithms for geodesics. See also the Wikipedia article on geodesics on an ellipsoid.
Solving the geodesic problem is slower than solving for the great-circle. But it's still very fast (about 1 μs per calculation), so this shouldn't be a reason to prefer great circle distances.
ADENDUM
Here is the Java package which implements my algorithm for finding geodesic distances. Unlike Vincenty's method, this is accurate to round-off and converges everywhere.