Point to point distances are needed all throughout applications we've written, most notably, Smache. So getting accurate distance values was obviously very important. Unfortunately, the pythagorean theorem doesn't cut it when we're looking for distances across the surface of a planet, so a slightly more complicated function had to be used.
What we wanted was the great circle distance. Great circles are lines that encircle a sphere leaving equal-sized halves on either side. Any two points on a sphere can be connected with one and only one great circle, except for the cases where the points are exactly opposite each other or exactly the same in which infinite great circles can pass through the two points. To come up with the distance between two points along their respective great circle, we need a little trig and the haversine formula.
NOTE: Smache was written in .NET with C# so my examples are also, but I included the javascript representation too because, well, javascript rocks.
The input parameters to the distance calculation function are the latitude and longitude for point A and the latitude and longitude for point B. We also needed a radius constant which was set at 6372 kilometers, widely agreed upon as the average radius of the earth. Note that because of the irregular shape of the earth, the resulting distance may be off by as much as half a percent.
const double EARTH_RADIUS = 6372;
double GreatCircleDistance(double latA, double longA, double latB, double longB) { }The haversine formula expects values to be in radians rather than degrees, so the first task is to make that conversion, simply dividing by 180 over pi.
double latARad = latA / (180 / Math.PI);
double longARad = longA / (180 / Math.PI);
double latBRad = latB / (180 / Math.PI);
double longBRad = longB / (180 / Math.PI);The formula will also use delta values, so in the interest of comprehension we'll create those variables as well.
double latDelta = latBRad - latARad;
double longDelta = longBRad - longARad;
Next is where the haversine comes in. In order to eliminate duplicate equations I've split the formula into two statements.
double haversineA = Math.Pow(Math.Sin(latDelta / 2), 2) + Math.Cos(latARad) * Math.Cos(latBRad) * Math.Pow(Math.Sin(longDelta / 2), 2);
double haversineB = 2 * Math.Atan2(Math.Sqrt(haversineA), Math.Sqrt(1 - haversineA));
Finally we simply need to multiply by the radius to get the final distance.
double finalDistance = haversineB * EARTH_RADIUS;So the function ends up like so...
const double EARTH_RADIUS = 6372;
double GreatCircleDistance(double latA, double longA, double latB, double longB) {
double latARad = latA / (180 / Math.PI);
double longARad = longA / (180 / Math.PI);
double latBRad = latB / (180 / Math.PI);
double longBRad = longB / (180 / Math.PI);
double latDelta = latBRad - latARad;
double longDelta = longBRad - longARad;
double haversineA = Math.Pow(Math.Sin(latDelta / 2), 2) + Math.Cos(latARad) * Math.Cos(latBRad) * Math.Pow(Math.Sin(longDelta / 2), 2);
double haversineB = 2 * Math.Atan2(Math.Sqrt(haversineA), Math.Sqrt(1 - haversineA));
double finalDistance = haversineB * EARTH_RADIUS;
return finalDistance;
}or in javascript...
var EARTH_RADIUS = 6372;
function GreatCircleDistance(latA, longA, latB, longB) {
var latARad = latA / (180 / Math.PI);
var longARad = longA / (180 / Math.PI);
var latBRad = latB / (180 / Math.PI);
var longBRad = longB / (180 / Math.PI);
var latDelta = latBRad - latARad;
var longDelta = longBRad - longARad;
var haversineA = Math.pow(Math.sin(latDelta / 2), 2) + Math.cos(latARad) * Math.cos(latBRad) * Math.pow(Math.sin(longDelta / 2), 2);
var haversineB = 2 * Math.atan2(Math.sqrt(haversineA), Math.sqrt(1 - haversineA));
var finalDistance = haversineB * EARTH_RADIUS;
return finalDistance;
}And if you're into one-liners...
double GreatCircleDistance(double latA, double longA, double latB, double longB) {
return 6372 * 2 * Math.Atan2(Math.Sqrt(Math.Pow(Math.Sin(((latB / (180 / Math.PI)) - (latA / (180 / Math.PI))) / 2), 2) + Math.Cos((latA / (180 / Math.PI))) * Math.Cos((latB / (180 / Math.PI))) * Math.Pow(Math.Sin(((longB / (180 / Math.PI)) - (longA / (180 / Math.PI))) / 2), 2)), Math.Sqrt(1 - (Math.Pow(Math.Sin(((latB / (180 / Math.PI)) - (latA / (180 / Math.PI))) / 2), 2) + Math.Cos((latA / (180 / Math.PI))) * Math.Cos((latB / (180 / Math.PI))) * Math.Pow(Math.Sin(((longB / (180 / Math.PI)) - (longA / (180 / Math.PI))) / 2), 2))));
}
function GreatCircleDistance(latA, longA, latB, longB) {
return 6372 * 2 * Math.atan2(Math.sqrt(Math.pow(Math.sin(((latB / (180 / Math.PI)) - (latA / (180 / Math.PI))) / 2), 2) + Math.cos((latA / (180 / Math.PI))) * Math.cos((latB / (180 / Math.PI))) * Math.pow(Math.sin(((longB / (180 / Math.PI)) - (longA / (180 / Math.PI))) / 2), 2)), Math.sqrt(1 - (Math.pow(Math.sin(((latB / (180 / Math.PI)) - (latA / (180 / Math.PI))) / 2), 2) + Math.cos((latA / (180 / Math.PI))) * Math.cos((latB / (180 / Math.PI))) * Math.pow(Math.sin(((longB / (180 / Math.PI)) - (longA / (180 / Math.PI))) / 2), 2))));
}