Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Oct 18, 2021

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.

See also

Jan 31, 2016

Direct and Inverse

Thaddeus Vincenty (1920--2002) Image source: IAG Newsletter

Thaddeus Vincenty (1920--2002) Image source: IAG Newsletter

There are two standard problems in geodesy:

  • Direct geodetic problem -- Given a point on the Earth (latitude and longitude) and the direction (azimuth) and distance from that point to a second point, determine the latitude and longitude of that second point.
  • Inverse geodetic problem -- Given two points on the Earth, determine the azimuth and distance between them.

Assuming the Earth is an oblate spheroid, both of these can be solved using the iterative algorithms by Thaddeus Vincenty [1]. Modern Fortran versions of Vincenty's algorithms can be found in the geodesy module of the Fortran Astrodynamics Toolkit [2].

The direct algorithm was listed in an earlier post. If we assume a WGS84 Earth ellipsoid (a=6378.137 km, f=1/298.257223563), an example direct calculation is:

$$ \begin{array}{rcll} \phi_1 &= & 29.97^{\circ} & \mathrm{(Latitude ~of ~ first ~ point)}\\ \lambda_1 &= & -95.35^{\circ} & \mathrm{(Longitude ~of ~ first ~point)} \\ Az &= & 20^{\circ} & \mathrm{(Azimuth)}\\ s &= & 50~\mathrm{km} & \mathrm{(Distance)} \end{array} $$

The latitude and longitude of the second point is:

$$ \begin{array}{rcl} \phi_2 &= & 30.393716^{\circ} \\ \lambda_2 &= & -95.172057^{\circ} \\ \end{array} $$

The more complicated inverse algorithm is listed here (based on the USNGS code from [1]):

subroutine inverse(a, rf, b1, l1, b2, l2, faz, baz, s, it, sig, lam, kind)

    use, intrinsic :: iso_fortran_env, wp => real64

    implicit none

    real(wp), intent(in) :: a !! Equatorial semimajor axis
    real(wp), intent(in) :: rf !! reciprocal flattening (1/f)
    real(wp), intent(in) :: b1 !! latitude of point 1 (rad, positive north)
    real(wp), intent(in) :: l1 !! longitude of point 1 (rad, positive east)
    real(wp), intent(in) :: b2 !! latitude of point 2 (rad, positive north)
    real(wp), intent(in) :: l2 !! longitude of point 2 (rad, positive east)
    real(wp), intent(out) :: faz !! Forward azimuth (rad, clockwise from north)
    real(wp), intent(out) :: baz !! Back azimuth (rad, clockwise from north)
    real(wp), intent(out) :: s !! Ellipsoidal distance
    integer, intent(out) :: it !! iteration count
    real(wp), intent(out) :: sig !! spherical distance on auxiliary sphere
    real(wp), intent(out) :: lam !! longitude difference on auxiliary sphere
    integer, intent(out) :: kind !! solution flag: kind=1, long-line; kind=2, antipodal

    real(wp) :: beta1, beta2, biga, bigb, bige, bigf, boa, &
                c, cosal2, coslam, cossig, costm, costm2, &
                cosu1, cosu2, d, dsig, ep2, l, prev, sinal, &
                sinlam, sinsig, sinu1, sinu2, tem1, tem2, &
                temp, test, z

    real(wp), parameter :: pi = acos(-1.0_wp)
    real(wp), parameter :: two_pi = 2.0_wp*pi
    real(wp), parameter :: tol = 1.0e-14_wp !! convergence tolerance
    real(wp), parameter :: eps = 1.0e-15_wp !! tolerance for zero

    boa = 1.0_wp - 1.0_wp/rf

    beta1 = atan(boa*tan(b1))
    sinu1 = sin(beta1)
    cosu1 = cos(beta1)
    beta2 = atan(boa*tan(b2))
    sinu2 = sin(beta2)
    cosu2 = cos(beta2)

    l = l2 - l1
    if (l > pi) l = l - two_pi
    if (l < -pi) l = l + two_pi
    prev = l
    test = l
    it = 0
    kind = 1
    lam = l

    longline: do ! long-line loop (kind=1)

        sinlam = sin(lam)
        coslam = cos(lam)
        temp = cosu1*sinu2 - sinu1*cosu2*coslam
        sinsig = sqrt((cosu2*sinlam)**2 + temp**2)
        cossig = sinu1*sinu2 + cosu1*cosu2*coslam
        sig = atan2(sinsig, cossig)
        if (abs(sinsig) < eps) then
            sinal = cosu1*cosu2*sinlam/sign(eps, sinsig)
        else
            sinal = cosu1*cosu2*sinlam/sinsig
        end if
        cosal2 = -sinal**2 + 1.0_wp
        if (abs(cosal2) < eps) then
            costm = -2.0_wp*(sinu1*sinu2/sign(eps, &
                                              cosal2)) + cossig
        else
            costm = -2.0_wp*(sinu1*sinu2/cosal2) + &
                    cossig
        end if
        costm2 = costm*costm
        c = ((-3.0_wp*cosal2 + 4.0_wp)/rf + 4.0_wp)* &
            cosal2/rf/16.0_wp

        antipodal: do ! antipodal loop (kind=2)

            it = it + 1
            d = (((2.0_wp*costm2 - 1.0_wp)*cossig*c + &
                  costm)*sinsig*c + sig)*(1.0_wp - c)/rf
            if (kind == 1) then
                lam = l + d*sinal
                if (abs(lam - test) >= tol) then
                    if (abs(lam) > pi) then
                        kind = 2
                        lam = pi
                        if (l < 0.0_wp) lam = -lam
                        sinal = 0.0_wp
                        cosal2 = 1.0_wp
                        test = 2.0_wp
                        prev = test
                        sig = pi - abs(atan(sinu1/cosu1) + &
                                    atan(sinu2/cosu2))
                        sinsig = sin(sig)
                        cossig = cos(sig)
                        c = ((-3.0_wp*cosal2 + 4.0_wp)/rf + &
                            4.0_wp)*cosal2/rf/16.0_wp
                        if (abs(sinal - prev) < tol) exit longline
                        if (abs(cosal2) < eps) then
                            costm = -2.0_wp*(sinu1*sinu2/ &
                                            sign(eps, cosal2)) + cossig
                        else
                            costm = -2.0_wp*(sinu1*sinu2/cosal2) + &
                                    cossig
                        end if
                        costm2 = costm*costm
                        cycle antipodal
                    end if
                    if (((lam - test)*(test - prev)) < 0.0_wp .and. it > 5) &
                        lam = (2.0_wp*lam + 3.0_wp*test + prev)/6.0_wp
                    prev = test
                    test = lam
                    cycle longline
                end if
            else
                sinal = (lam - l)/d
                if (((sinal - test)*(test - prev)) < 0.0_wp .and. it > 5) &
                    sinal = (2.0_wp*sinal + 3.0_wp*test + prev)/6.0_wp
                prev = test
                test = sinal
                cosal2 = -sinal**2 + 1.0_wp
                sinlam = sinal*sinsig/(cosu1*cosu2)
                coslam = -sqrt(abs(-sinlam**2 + 1.0_wp))
                lam = atan2(sinlam, coslam)
                temp = cosu1*sinu2 - sinu1*cosu2*coslam
                sinsig = sqrt((cosu2*sinlam)**2 + temp**2)
                cossig = sinu1*sinu2 + cosu1*cosu2*coslam
                sig = atan2(sinsig, cossig)
                c = ((-3.0_wp*cosal2 + 4.0_wp)/rf + 4.0_wp)* &
                    cosal2/rf/16.0_wp
                if (abs(sinal - prev) >= tol) then
                if (abs(cosal2) < eps) then
                    costm = -2.0_wp*(sinu1*sinu2/ &
                                     sign(eps, cosal2)) + cossig
                else
                    costm = -2.0_wp*(sinu1*sinu2/cosal2) + &
                            cossig
                end if
                costm2 = costm*costm
                cycle antipodal
                end if
            end if

            exit longline

        end do antipodal

    end do longline

    if (kind == 2) then ! antipodal
        faz = sinal/cosu1
        baz = sqrt(-faz**2 + 1.0_wp)
        if (temp < 0.0_wp) baz = -baz
        faz = atan2(faz, baz)
        tem1 = -sinal
        tem2 = sinu1*sinsig - cosu1*cossig*baz
        baz = atan2(tem1, tem2)
    else ! long-line
        tem1 = cosu2*sinlam
        tem2 = cosu1*sinu2 - sinu1*cosu2*coslam
        faz = atan2(tem1, tem2)
        tem1 = -cosu1*sinlam
        tem2 = sinu1*cosu2 - cosu1*sinu2*coslam
        baz = atan2(tem1, tem2)
    end if

    if (faz < 0.0_wp) faz = faz + two_pi
    if (baz < 0.0_wp) baz = baz + two_pi
    ep2 = 1.0_wp/(boa*boa) - 1.0_wp
    bige = sqrt(1.0_wp + ep2*cosal2)
    bigf = (bige - 1.0_wp)/(bige + 1.0_wp)
    biga = (1.0_wp + bigf*bigf/4.0_wp)/(1.0_wp - bigf)
    bigb = bigf*(1.0_wp - 0.375_wp*bigf*bigf)
    z = bigb/6.0_wp*costm*(-3.0_wp + 4.0_wp* &
                           sinsig**2)*(-3.0_wp + 4.0_wp*costm2)
    dsig = bigb*sinsig*(costm + bigb/4.0_wp* &
                        (cossig*(-1.0_wp + 2.0_wp*costm2) - z))
    s = (boa*a)*biga*(sig - dsig)

end subroutine inverse

For the same ellipsoid, an example inverse calculation is:

$$ \begin{array}{rcll} \phi_1 &= & 29.97^{\circ} & \mathrm{(Latitude ~of ~ first ~ point)}\\ \lambda_1 &= & -95.35^{\circ} & \mathrm{(Longitude ~of ~ first ~point)} \\ \phi_2 &= & 40.77^{\circ} & \mathrm{(Latitude ~of ~ second ~ point)}\\ \lambda_2 &= & -73.98^{\circ} & \mathrm{(Longitude ~of ~ second ~point)} \end{array} $$

The azimuth and distance from the first to the second point is:

$$ \begin{array}{rcl} Az &= & 52.400056^{\circ} & \mathrm{(Azimuth)}\\ s &= & 2272.497~\mathrm{km} & \mathrm{(Distance)} \end{array} $$

map

References

Jul 03, 2014

Direct Geodetic Problem

The "direct" geodetic problem is: given the latitude and longitude of one point and the azimuth and distance to a second point, determine the latitude and longitude of that second point. The solution can be obtained using the algorithm by Polish American geodesist Thaddeus Vincenty [1]. A modern Fortran implementation is given below:

subroutine direct(a,f,glat1,glon1,glat2,glon2,faz,baz,s)

use, intrinsic :: iso_fortran_env, wp => real64

implicit none

real(wp),intent(in) :: a !semimajor axis of ellipsoid [m]
real(wp),intent(in) :: f !flattening of ellipsoid [-]
real(wp),intent(in) :: glat1 !latitude of 1 [rad]
real(wp),intent(in) :: glon1 !longitude of 1 [rad]
real(wp),intent(in) :: faz !forward azimuth 1->2 [rad]
real(wp),intent(in) :: s !distance from 1->2 [m]
real(wp),intent(out) :: glat2 !latitude of 2 [rad]
real(wp),intent(out) :: glon2 !longitude of 2 [rad]
real(wp),intent(out) :: baz !back azimuth 2->1 [rad]

real(wp) :: r,tu,sf,cf,cu,su,sa,csa,c2a,x,c,d,y,sy,cy,cz,e

real(wp),parameter :: pi = acos(-1.0_wp)
real(wp),parameter :: eps = 0.5e-13_wp

r = 1.0_wp-f
tu = r*sin(glat1)/cos(glat1)
sf = sin(faz)
cf = cos(faz)
baz = 0.0_wp
if (cf/=0.0_wp) baz = atan2(tu,cf)*2.0_wp
cu = 1.0_wp/sqrt(tu*tu+1.0_wp)
su = tu*cu
sa = cu*sf
c2a = -sa*sa+1.0_wp
x = sqrt((1.0_wp/r/r-1.0_wp)*c2a+1.0_wp)+1.0_wp
x = (x-2.0_wp)/x
c = 1.0_wp-x
c = (x*x/4.0_wp+1.0_wp)/c
d = (0.375_wp*x*x-1.0_wp)*x
tu = s/r/a/c
y = tu
do
    sy = sin(y)
    cy = cos(y)
    cz = cos(baz+y)
    e = cz*cz*2.0_wp-1.0_wp
    c = y
    x = e*cy
    y = e+e-1.0_wp
    y = (((sy*sy*4.0_wp-3.0_wp)*y*cz*d/6.0_wp+x)*d/4.0_wp-cz)*sy*d+tu
    if (abs(y-c)<=eps) exit
end do
baz = cu*cy*cf-su*sy
c = r*sqrt(sa*sa+baz*baz)
d = su*cy+cu*sy*cf
glat2 = atan2(d,c)
c = cu*cy-su*sy*cf
x = atan2(sy*sf,c)
c = ((-3.0_wp*c2a+4.0_wp)*f+4.0_wp)*c2a*f/16.0_wp
d = ((e*cy*c+cz)*sy*c+y)*sa
glon2 = glon1+x-(1.0_wp-c)*d*f
baz = atan2(sa,baz)+pi

end subroutine direct

References

  1. T. Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations", Survey Review XXII. 176, April 1975. [sourcecode from the U.S. National Geodetic Survey].