The haversine formula is useful in calculating the distance between two locations on the Earth, when the latitudes & longitudes of the two locations are provided.
An implementation of calcDistance function based on haversine formula in JavaScript is as follows: (the degToRad function is solely a unit converter)
function degToRad(degrees) {
return degrees * Math.PI / 180;
}
function calcDistance(lat1, lon1, lat2, lon2) {
var earthRadiusKm = 6371;
var dLat = degToRad(lat2-lat1);
var dLon = degToRad(lon2-lon1);
lat1 = degToRad(lat1);
lat2 = degToRad(lat2);
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
return earthRadiusKm * c;
}The algorithm in the calcDistance function and the meaning of intermediate values (a and c) can be explained using the following geometrical derivation.
The following is a isosceles triangle, with the length of the segment AO and BO equal to 1.
(please ignore the notation and
as the image is shared in explaining another phase)
We can tell that: defining the angle ∠AOB as θ, then the length of the segment AB that is oppsite to the angle will be
A unit sphere (with radius as 1) with two points on the surface of the sphere A(lat1, lon1) and B(lat2, lon2), as well as some auxillary points and lines are drawn as follows. Point N is the north pole and point O is the center of the sphere.
We can tell that:
- in the triangle formed by point A, O and C, the angle ∠AOC shows the difference between latitude values of A and B. Therefore the angle will be
(or
provided that
lat1andlat2have been converted to radian) - the radius for a unit sphere is 1. Therefore, the length of chord AC will be
according to the conclusion in phase1
- as the difference of latitude between AC and BD are the same, the length of chord BD will be
- similarly, in the triangle EOF, the angle ∠EOF shows the difference between longitude values of A and B. Therefore, the length of chord EF will be
- AG is a perpendicular line to OE in the triangle AOE, the length of OG will be
, which is
- the ratio between the length of OG and OE is the length of OG since OE is the radius 1
- OG is the radius of the circular plane parallel to the ecliptic plane at lat1. As lat1 goes near the north pole, OG goes to 0
- the chord EF is at the ecliptic plane and the chord AD is at the circular plane at lat1. As lat1 goes near the north pole, AD goes to 0
- we can tell that OG and AD will decrese to 0 at the same rate, which is the value of OG. Therefore, the length of AD will be
, which is
- similarly, the length of BC will be
At the end of derivation phase 2, we get the length of all 4 sides of a quadrilateral ▱ABCD
The 2-dimensional shape ▱ABCD derived from the previous phase is as follows. We need to calculate the length of segment AB (which is the chord of AB from the sphere).
We can tell that:
- AH is a perpendicular line to BC
as the shape is a isosceles trapezoid
by the Pythagorean Theorem
by the Pythagorean Theorem
- By substituting results from the previous phase, eventually we will get that:
- For convenience in the next phase, we define an intermediate value a to be half of the chord length AB squared. That is:
In the sphere in phase 2, the triangle formed by point A, B and O is as follows:
We can tell that:
- OM is a perpendicular line to AB. As the triangle is isosceles, M will be the middle point of AB
- from the previous phase,
. Therefore,
by the Pythagorean Thereom
, and
The angle ∠AOB in radian tells about the spherical distance between A and B. In the case of a unit sphere, the distance between A and B will be ∠AOB in radian. In the case of the Earth, the answer will be