Degenerate Conic

Great Circle

great-circle

If we assume that the Earth is a sphere of radius $r$, the shortest distance between two sites ($[\phi_1,\lambda_1]$ and $[\phi_2,\lambda_2]$) on the surface is the great-circle distance, which can be computed in several ways:

Law of Cosines equation:

$$ d = r \arccos\bigl(\sin\phi_1\cdot\sin\phi_2+\cos\phi_1\cdot\cos\phi_2\cdot\cos(\Delta\lambda)\bigr) $$

This equation, while mathematically correct, is not usually recommended for use, since it is ill-conditioned for sites that are very close to each other.

Using the haversine equation:

$$ d = 2 r \arcsin \sqrt{\sin^2\left(\frac{\Delta\phi}{2}\right)+\cos{\phi_1}\cdot\cos{\phi_2}\cdot\sin^2\left(\frac{\Delta\lambda}{2}\right)} $$

The haversine one is better, but can be ill-conditioned for near antipodal points.

Using the Vincenty equation:

$$ d = r \arctan \frac{\sqrt{\left(\cos\phi_2\cdot\sin(\Delta\lambda)\right)^2+\left(\cos\phi_1\cdot\sin\phi_2-\sin\phi_1\cdot\cos\phi_2\cdot\cos(\Delta\lambda)\right)^2}}{\sin\phi_1\cdot\sin\phi_2+\cos\phi_1\cdot\cos\phi_2\cdot\cos(\Delta\lambda)} $$

This one is accurate for any two points on a sphere.

Code

Fortran implementations of these algorithms are given here:

Law of Cosines:

pure function gcdist(r,lat1,long1,lat2,long2) &
              result(d)

implicit none

real(wp) :: d !! great circle distance from 1 to 2 [km]
real(wp),intent(in) :: r     !! radius of the body [km]
real(wp),intent(in) :: long1 !! longitude of first site [rad]
real(wp),intent(in) :: lat1  !! latitude of the first site [rad]
real(wp),intent(in) :: long2 !! longitude of the second site [rad]
real(wp),intent(in) :: lat2  !! latitude of the second site [rad]

real(wp) :: clat1,clat2,slat1,slat2,dlon,cdlon

clat1 = cos(lat1)
clat2 = cos(lat2)
slat1 = sin(lat1)
slat2 = sin(lat2)
dlon = long1-long2
cdlon = cos(dlon)

d = r * acos(slat1*slat2+clat1*clat2*cdlon)

end function gcdist

Haversine:

pure function gcdist_haversine(r,lat1,long1,lat2,long2) &
              result(d)

implicit none

real(wp) :: d !! great circle distance from 1 to 2 [km]
real(wp),intent(in) :: r     !! radius of the body [km]
real(wp),intent(in) :: long1 !! longitude of first site [rad]
real(wp),intent(in) :: lat1  !! latitude of the first site [rad]
real(wp),intent(in) :: long2 !! longitude of the second site [rad]
real(wp),intent(in) :: lat2  !! latitude of the second site [rad]

real(wp) :: dlat,clat1,clat2,dlon,a

dlat  = lat1-lat2
clat1 = cos(lat1)
clat2 = cos(lat2)
dlon  = long1-long2

a = sin(dlat/2.0_wp)**2 + clat1*clat2*sin(dlon/2.0_wp)**2
d = r * 2.0_wp * asin(min(1.0_wp,sqrt(a)))

end function gcdist_haversine

Vincenty:

pure function gcdist_vincenty(r,lat1,long1,lat2,long2) &
              result(d)

implicit none

real(wp) :: d !! great circle distance from 1 to 2 [km]
real(wp),intent(in) :: r     !! radius of the body [km]
real(wp),intent(in) :: long1 !! longitude of first site [rad]
real(wp),intent(in) :: lat1  !! latitude of the first site [rad]
real(wp),intent(in) :: long2 !! longitude of the second site [rad]
real(wp),intent(in) :: lat2  !! latitude of the second site [rad]

real(wp) :: c1,s1,c2,s2,dlon,clon,slon

c1   = cos(lat1)
s1   = sin(lat1)
c2   = cos(lat2)
s2   = sin(lat2)
dlon = long1-long2
clon = cos(dlon)
slon = sin(dlon)

d = r*atan2(sqrt((c2*slon)**2+(c1*s2-s1*c2*clon)**2), &
    (s1*s2+c1*c2*clon))

end function gcdist_vincenty

Test cases

Using a radius of the Earth of 6,378,137 meters, here are a few test cases (using double precision (real64) arithmetic):

Close points

The distance between the following very close points $([\phi_1,\lambda_1], [\phi_2,\lambda_2]) =([0, 0.000001], [0.0, 0.0])$ (radians) is:

Algorithm	Procedure	Result (meters)
Law of Cosines	`gcdist`	6.3784205037462689
Haversine	`gcdist_haversine`	6.3781369999999997
Vincenty	`gcdist_vincenty`	6.3781369999999997

Houston to New York

The distance between the following two points $([\phi_1,\lambda_1], [\phi_2,\lambda_2]) =([29.97, -95.35], [40.77, -73.98])$ (deg) is:

Algorithm	Procedure	Result (meters)
Law of Cosines	`gcdist`	2272779.3057236285
Haversine	`gcdist_haversine`	2272779.3057236294
Vincenty	`gcdist_vincenty`	2272779.3057236290

Antipodal points

The distance between the following two antipodal points $([\phi_1,\lambda_1], [\phi_2,\lambda_2]) = ([0, 0], [0.0, \pi/2])$ (radians) is:

Algorithm	Procedure	Result (meters)
Law of Cosines	`gcdist`	20037508.342789244
Haversine	`gcdist_haversine`	20037508.342789244
Vincenty	`gcdist_vincenty`	20037508.342789244

Nearly antipodal points

The distance between the following two nearly antipodal points $([\phi_1,\lambda_1], [\phi_2,\lambda_2]) = ([10^{-8}, 10^{-8}], [0.0, \pi/2])$ (radians) is:

Algorithm	Procedure	Result (meters)
Law of Cosines	`gcdist`	20037508.342789244
Haversine	`gcdist_haversine`	20037508.342789244
Vincenty	`gcdist_vincenty`	20037508.252588764

Other methods

For a more accurate result, using an oblate-spheroid model of the Earth, the Vincenty inverse geodetic algorithm can be used (the equation above is just a special case of the more general algorithm when the ellipsoidal major and minor axes are equal). For the Houston to New York case, the difference between the spherical and oblate-spheroid equations is 282 meters.

Algorithms • Modern Fortran Programming • Orbital Mechanics

Oct 18, 2021