Geodesics on the Ellipsoid

Back to Transverse Mercator Projection. Forward to Geocentric coordinates. Up to Contents.

GeographicLib::Geodesic and GeographicLib::GeodesicLine provide accurate solutions to the direct and inverse geodesic problems. The Geod utility provides an interface to these classes. GeographicLib::AzimuthalEquidistant implements the azimuthal equidistant projection in terms of geodesics. GeographicLib::CassiniSoldner implements a transverse cylindrical equidistant projection in terms of geodesics. The EquidistantTest utility provides an interface to these projections.

References

Test data for geodesics

A test set a geodesics is available at

This is about 39 MB (compressed). This consists of a set of geodesics for the WGS84 ellipsoid. A subset of this (consisting of 1/50 of the members — about 690 kB, compressed) is available at

Each line of the test set gives 9 space delimited numbers

These are computed using as direct geodesic calculations with the given lat1, lon1, azi1, and s12. The distance s12 always corresponds to an arc length a12 <= 180o, so the given geodesics give the shortest paths from point 1 to point 2. For simplicity and without loss of generality, lat1 is chosen in [0o, 90o], lon1 is taken to be zero, azi1 is chosen in [0o, 180o]. Furthermore, lat1 and azi1 are taken to be multiples of 10-12 deg and s12 is a multiple of 0.1 um in [0 m, 20003931.4586254 m]. This results lon2 in [0o, 180o] and azi2 in [0o, 180o].

The direct calculation uses an expansion of the geodesic equations accurate to f30 (approximately 1 part in 1050) and is computed with with Maxima's bfloats and fpprec set to 100 (so the errors in the data are probably 1/2 of the values quoted above).

The contents of the file are as follows:

(a total of 500000 entries). The values for s12 for the geodesics running between vertices are truncated to a multiple of 0.1 pm and this is used to determine point 2.

This data can be fed to the Geod utility as follows

Add, e.g., "-p 6", to the call to Geod to change the precision of the output. Adding "-f" causes Geod to print all 9 fields specifying the geodesic (in the same order as the test set).

Expansions for geodesics

We give here the series expansions for the various geodesic integrals valid to order f10. In this release of the code, we use a 6th-order expansions. This is sufficient to maintain accuracy for doubles for the SRMmax ellipsoid (a = 6400 km, f = 1/150). However, the preprocessor macro GEOD_ORD can be used to select any order up to 8. (If using long doubles, with a 64-bit fraction, the default order is 7.) The series expanded to order f30 are given in geodseries30.html.

In the formulas below ^ indicates exponentiation (f^3 = f*f*f) and / indicates real division (3/5 = 0.6). The equations need to be converted to Horner form, but are here left in expanded form so that they can be easily truncated to lower order. These expansions were obtained using the the Maxima code, geod.mac.

In the expansions below, we have

The formula for distance is

    s/b = I1(sigma)

where

    I1(sigma) = A1 (sigma + B1(sigma))
    B1(sigma) = sumj = 1 C1j sin(2 j sigma)

and

A1 = (1 + 1/4 * eps^2
        + 1/64 * eps^4
        + 1/256 * eps^6
        + 25/16384 * eps^8
        + 49/65536 * eps^10) / (1 - eps);
C1[1] = - 1/2 * eps
        + 3/16 * eps^3
        - 1/32 * eps^5
        + 19/2048 * eps^7
        - 3/4096 * eps^9;
C1[2] = - 1/16 * eps^2
        + 1/32 * eps^4
        - 9/2048 * eps^6
        + 7/4096 * eps^8
        + 1/65536 * eps^10;
C1[3] = - 1/48 * eps^3
        + 3/256 * eps^5
        - 3/2048 * eps^7
        + 17/24576 * eps^9;
C1[4] = - 5/512 * eps^4
        + 3/512 * eps^6
        - 11/16384 * eps^8
        + 3/8192 * eps^10;
C1[5] = - 7/1280 * eps^5
        + 7/2048 * eps^7
        - 3/8192 * eps^9;
C1[6] = - 7/2048 * eps^6
        + 9/4096 * eps^8
        - 117/524288 * eps^10;
C1[7] = - 33/14336 * eps^7
        + 99/65536 * eps^9;
C1[8] = - 429/262144 * eps^8
        + 143/131072 * eps^10;
C1[9] = - 715/589824 * eps^9;
C1[10] = - 2431/2621440 * eps^10;

The function tau(sigma) = s/(b A1) = sigma + B1(sigma) may be inverted by series reversion giving

    sigma(tau) = tau + sumj = 1 C1'j sin(2 j tau)

where

C1'[1] = + 1/2 * eps
         - 9/32 * eps^3
         + 205/1536 * eps^5
         - 4879/73728 * eps^7
         + 9039/327680 * eps^9;
C1'[2] = + 5/16 * eps^2
         - 37/96 * eps^4
         + 1335/4096 * eps^6
         - 86171/368640 * eps^8
         + 4119073/28311552 * eps^10;
C1'[3] = + 29/96 * eps^3
         - 75/128 * eps^5
         + 2901/4096 * eps^7
         - 443327/655360 * eps^9;
C1'[4] = + 539/1536 * eps^4
         - 2391/2560 * eps^6
         + 1082857/737280 * eps^8
         - 2722891/1548288 * eps^10;
C1'[5] = + 3467/7680 * eps^5
         - 28223/18432 * eps^7
         + 1361343/458752 * eps^9;
C1'[6] = + 38081/61440 * eps^6
         - 733437/286720 * eps^8
         + 10820079/1835008 * eps^10;
C1'[7] = + 459485/516096 * eps^7
         - 709743/163840 * eps^9;
C1'[8] = + 109167851/82575360 * eps^8
         - 550835669/74317824 * eps^10;
C1'[9] = + 83141299/41287680 * eps^9;
C1'[10] = + 9303339907/2972712960 * eps^10;

The reduced length is given by

    m/b = sqrt(1 + k2 sin2sigma2) cos sigma1 sin sigma2
            - sqrt(1 + k2 sin2sigma1) sin sigma1 cos sigma2
            - cos sigma1 cos sigma2 (J(sigma2) - J(sigma1))

where

    J(sigma) = I1(sigma) - I2(sigma)
    I2(sigma) = A2 (sigma + B2(sigma))
    B2(sigma) = sumj = 1 C2j sin(2 j sigma)

A2 = (1 + 1/4 * eps^2
        + 9/64 * eps^4
        + 25/256 * eps^6
        + 1225/16384 * eps^8
        + 3969/65536 * eps^10) * (1 - eps);
C2[1] = + 1/2 * eps
        + 1/16 * eps^3
        + 1/32 * eps^5
        + 41/2048 * eps^7
        + 59/4096 * eps^9;
C2[2] = + 3/16 * eps^2
        + 1/32 * eps^4
        + 35/2048 * eps^6
        + 47/4096 * eps^8
        + 557/65536 * eps^10;
C2[3] = + 5/48 * eps^3
        + 5/256 * eps^5
        + 23/2048 * eps^7
        + 191/24576 * eps^9;
C2[4] = + 35/512 * eps^4
        + 7/512 * eps^6
        + 133/16384 * eps^8
        + 47/8192 * eps^10;
C2[5] = + 63/1280 * eps^5
        + 21/2048 * eps^7
        + 51/8192 * eps^9;
C2[6] = + 77/2048 * eps^6
        + 33/4096 * eps^8
        + 2607/524288 * eps^10;
C2[7] = + 429/14336 * eps^7
        + 429/65536 * eps^9;
C2[8] = + 6435/262144 * eps^8
        + 715/131072 * eps^10;
C2[9] = + 12155/589824 * eps^9;
C2[10] = + 46189/2621440 * eps^10;

The longitude is given in terms of the spherical longitude by

    lambda = omega - f sin alpha0 I3(sigma)

where

    I3(sigma) = A3 (sigma + B3(sigma))
    B3(sigma) = sumj = 1 C3j sin(2 j sigma)

and

A3 = 1 - (1/2 - 1/2 * n) * eps
       - (1/4 + 1/8 * n - 3/8 * n^2) * eps^2
       - (1/16 + 3/16 * n + 1/16 * n^2 - 5/16 * n^3) * eps^3
       - (3/64 + 1/32 * n + 5/32 * n^2 + 5/128 * n^3 - 35/128 * n^4) * eps^4
       - (3/128 + 5/128 * n + 5/256 * n^2 + 35/256 * n^3 + 7/256 * n^4) * eps^5
       - (5/256 + 15/1024 * n + 35/1024 * n^2 + 7/512 * n^3) * eps^6
       - (25/2048 + 35/2048 * n + 21/2048 * n^2) * eps^7
       - (175/16384 + 35/4096 * n) * eps^8
       - 245/32768 * eps^9;
C3[1] = + (1/4 - 1/4 * n) * eps
        + (1/8 - 1/8 * n^2) * eps^2
        + (3/64 + 3/64 * n - 1/64 * n^2 - 5/64 * n^3) * eps^3
        + (5/128 + 1/64 * n + 1/64 * n^2 - 1/64 * n^3 - 7/128 * n^4) * eps^4
        + (3/128 + 11/512 * n + 3/512 * n^2 + 1/256 * n^3 - 7/512 * n^4) * eps^5
        + (21/1024 + 5/512 * n + 13/1024 * n^2 + 1/512 * n^3) * eps^6
        + (243/16384 + 189/16384 * n + 83/16384 * n^2) * eps^7
        + (435/32768 + 109/16384 * n) * eps^8
        + 345/32768 * eps^9;
C3[2] = + (1/16 - 3/32 * n + 1/32 * n^2) * eps^2
        + (3/64 - 1/32 * n - 3/64 * n^2 + 1/32 * n^3) * eps^3
        + (3/128 + 1/128 * n - 9/256 * n^2 - 3/128 * n^3 + 7/256 * n^4) * eps^4
        + (5/256 + 1/256 * n - 1/128 * n^2 - 7/256 * n^3 - 3/256 * n^4) * eps^5
        + (27/2048 + 69/8192 * n - 39/8192 * n^2 - 47/4096 * n^3) * eps^6
        + (187/16384 + 39/8192 * n + 31/16384 * n^2) * eps^7
        + (287/32768 + 47/8192 * n) * eps^8
        + 255/32768 * eps^9;
C3[3] = + (5/192 - 3/64 * n + 5/192 * n^2 - 1/192 * n^3) * eps^3
        + (3/128 - 5/192 * n - 1/64 * n^2 + 5/192 * n^3 - 1/128 * n^4) * eps^4
        + (7/512 - 1/384 * n - 77/3072 * n^2 + 5/3072 * n^3 + 65/3072 * n^4) * eps^5
        + (3/256 - 1/1024 * n - 71/6144 * n^2 - 47/3072 * n^3) * eps^6
        + (139/16384 + 143/49152 * n - 383/49152 * n^2) * eps^7
        + (243/32768 + 95/49152 * n) * eps^8
        + 581/98304 * eps^9;
C3[4] = + (7/512 - 7/256 * n + 5/256 * n^2 - 7/1024 * n^3 + 1/1024 * n^4) * eps^4
        + (7/512 - 5/256 * n - 7/2048 * n^2 + 9/512 * n^3 - 21/2048 * n^4) * eps^5
        + (9/1024 - 43/8192 * n - 129/8192 * n^2 + 39/4096 * n^3) * eps^6
        + (127/16384 - 23/8192 * n - 165/16384 * n^2) * eps^7
        + (193/32768 + 3/8192 * n) * eps^8
        + 171/32768 * eps^9;
C3[5] = + (21/2560 - 9/512 * n + 15/1024 * n^2 - 7/1024 * n^3 + 9/5120 * n^4) * eps^5
        + (9/1024 - 15/1024 * n + 3/2048 * n^2 + 57/5120 * n^3) * eps^6
        + (99/16384 - 91/16384 * n - 781/81920 * n^2) * eps^7
        + (179/32768 - 55/16384 * n) * eps^8
        + 141/32768 * eps^9;
C3[6] = + (11/2048 - 99/8192 * n + 275/24576 * n^2 - 77/12288 * n^3) * eps^6
        + (99/16384 - 275/24576 * n + 55/16384 * n^2) * eps^7
        + (143/32768 - 253/49152 * n) * eps^8
        + 33/8192 * eps^9;
C3[7] = + (429/114688 - 143/16384 * n + 143/16384 * n^2) * eps^7
        + (143/32768 - 143/16384 * n) * eps^8
        + 429/131072 * eps^9;
C3[8] = + (715/262144 - 429/65536 * n) * eps^8
        + 429/131072 * eps^9;
C3[9] = + 2431/1179648 * eps^9;

The formula for area between the geodesic and the equator is given in Sec. 15 of arxiv:1102.1215 in terms of S,

    S = c2 alpha + e2 a2 cos alpha0 sin alpha0 I4(sigma)

where

    I4(sigma) = sumj = 0 C4j cos((2j + 1) sigma)

with

C4[0] = + (2/3 - 1/15 * ep2 + 4/105 * ep2^2 - 8/315 * ep2^3 + 64/3465 * ep2^4 - 128/9009 * ep2^5 + 512/45045 * ep2^6 - 1024/109395 * ep2^7 + 16384/2078505 * ep2^8 - 32768/4849845 * ep2^9)
        - (1/20 - 1/35 * ep2 + 2/105 * ep2^2 - 16/1155 * ep2^3 + 32/3003 * ep2^4 - 128/15015 * ep2^5 + 256/36465 * ep2^6 - 4096/692835 * ep2^7 + 8192/1616615 * ep2^8) * k2
        + (1/42 - 1/63 * ep2 + 8/693 * ep2^2 - 80/9009 * ep2^3 + 64/9009 * ep2^4 - 128/21879 * ep2^5 + 2048/415701 * ep2^6 - 4096/969969 * ep2^7) * k2^2
        - (1/72 - 1/99 * ep2 + 10/1287 * ep2^2 - 8/1287 * ep2^3 + 112/21879 * ep2^4 - 1792/415701 * ep2^5 + 512/138567 * ep2^6) * k2^3
        + (1/110 - 1/143 * ep2 + 4/715 * ep2^2 - 56/12155 * ep2^3 + 896/230945 * ep2^4 - 768/230945 * ep2^5) * k2^4
        - (1/156 - 1/195 * ep2 + 14/3315 * ep2^2 - 224/62985 * ep2^3 + 64/20995 * ep2^4) * k2^5
        + (1/210 - 1/255 * ep2 + 16/4845 * ep2^2 - 32/11305 * ep2^3) * k2^6
        - (1/272 - 1/323 * ep2 + 6/2261 * ep2^2) * k2^7
        + (1/342 - 1/399 * ep2) * k2^8
        - 1/420 * k2^9;
C4[1] = + (1/180 - 1/315 * ep2 + 2/945 * ep2^2 - 16/10395 * ep2^3 + 32/27027 * ep2^4 - 128/135135 * ep2^5 + 256/328185 * ep2^6 - 4096/6235515 * ep2^7 + 8192/14549535 * ep2^8) * k2
        - (1/252 - 1/378 * ep2 + 4/2079 * ep2^2 - 40/27027 * ep2^3 + 32/27027 * ep2^4 - 64/65637 * ep2^5 + 1024/1247103 * ep2^6 - 2048/2909907 * ep2^7) * k2^2
        + (1/360 - 1/495 * ep2 + 2/1287 * ep2^2 - 8/6435 * ep2^3 + 112/109395 * ep2^4 - 1792/2078505 * ep2^5 + 512/692835 * ep2^6) * k2^3
        - (1/495 - 2/1287 * ep2 + 8/6435 * ep2^2 - 112/109395 * ep2^3 + 1792/2078505 * ep2^4 - 512/692835 * ep2^5) * k2^4
        + (5/3276 - 1/819 * ep2 + 2/1989 * ep2^2 - 32/37791 * ep2^3 + 64/88179 * ep2^4) * k2^5
        - (1/840 - 1/1020 * ep2 + 4/4845 * ep2^2 - 8/11305 * ep2^3) * k2^6
        + (7/7344 - 7/8721 * ep2 + 2/2907 * ep2^2) * k2^7
        - (2/2565 - 4/5985 * ep2) * k2^8
        + 1/1540 * k2^9;
C4[2] = + (1/2100 - 1/3150 * ep2 + 4/17325 * ep2^2 - 8/45045 * ep2^3 + 32/225225 * ep2^4 - 64/546975 * ep2^5 + 1024/10392525 * ep2^6 - 2048/24249225 * ep2^7) * k2^2
        - (1/1800 - 1/2475 * ep2 + 2/6435 * ep2^2 - 8/32175 * ep2^3 + 112/546975 * ep2^4 - 1792/10392525 * ep2^5 + 512/3464175 * ep2^6) * k2^3
        + (1/1925 - 2/5005 * ep2 + 8/25025 * ep2^2 - 16/60775 * ep2^3 + 256/1154725 * ep2^4 - 1536/8083075 * ep2^5) * k2^4
        - (1/2184 - 1/2730 * ep2 + 1/3315 * ep2^2 - 16/62985 * ep2^3 + 32/146965 * ep2^4) * k2^5
        + (1/2520 - 1/3060 * ep2 + 4/14535 * ep2^2 - 8/33915 * ep2^3) * k2^6
        - (7/20400 - 7/24225 * ep2 + 2/8075 * ep2^2) * k2^7
        + (14/47025 - 4/15675 * ep2) * k2^8
        - 1/3850 * k2^9;
C4[3] = + (1/17640 - 1/24255 * ep2 + 2/63063 * ep2^2 - 8/315315 * ep2^3 + 16/765765 * ep2^4 - 256/14549535 * ep2^5 + 512/33948915 * ep2^6) * k2^3
        - (1/10780 - 1/14014 * ep2 + 2/35035 * ep2^2 - 4/85085 * ep2^3 + 64/1616615 * ep2^4 - 384/11316305 * ep2^5) * k2^4
        + (5/45864 - 1/11466 * ep2 + 1/13923 * ep2^2 - 16/264537 * ep2^3 + 32/617253 * ep2^4) * k2^5
        - (1/8820 - 1/10710 * ep2 + 8/101745 * ep2^2 - 16/237405 * ep2^3) * k2^6
        + (1/8976 - 1/10659 * ep2 + 2/24871 * ep2^2) * k2^7
        - (1/9405 - 2/21945 * ep2) * k2^8
        + 1/10010 * k2^9;
C4[4] = + (1/124740 - 1/162162 * ep2 + 2/405405 * ep2^2 - 4/984555 * ep2^3 + 64/18706545 * ep2^4 - 128/43648605 * ep2^5) * k2^4
        - (1/58968 - 1/73710 * ep2 + 1/89505 * ep2^2 - 16/1700595 * ep2^3 + 32/3968055 * ep2^4) * k2^5
        + (1/41580 - 1/50490 * ep2 + 8/479655 * ep2^2 - 16/1119195 * ep2^3) * k2^6
        - (7/242352 - 7/287793 * ep2 + 2/95931 * ep2^2) * k2^7
        + (7/220077 - 2/73359 * ep2) * k2^8
        - 1/30030 * k2^9;
C4[5] = + (1/792792 - 1/990990 * ep2 + 1/1203345 * ep2^2 - 16/22863555 * ep2^3 + 32/53348295 * ep2^4) * k2^5
        - (1/304920 - 1/370260 * ep2 + 4/1758735 * ep2^2 - 8/4103715 * ep2^3) * k2^6
        + (7/1283568 - 7/1524237 * ep2 + 2/508079 * ep2^2) * k2^7
        - (2/268983 - 4/627627 * ep2) * k2^8
        + 1/110110 * k2^9;
C4[6] = + (1/4684680 - 1/5688540 * ep2 + 4/27020565 * ep2^2 - 8/63047985 * ep2^3) * k2^6
        - (1/1516944 - 1/1801371 * ep2 + 2/4203199 * ep2^2) * k2^7
        + (2/1589445 - 4/3708705 * ep2) * k2^8
        - 1/520520 * k2^9;
C4[7] = + (1/26254800 - 1/31177575 * ep2 + 2/72747675 * ep2^2) * k2^7
        - (1/7335900 - 1/8558550 * ep2) * k2^8
        + 1/3403400 * k2^9;
C4[8] = + (1/141338340 - 1/164894730 * ep2) * k2^8
        - 1/34714680 * k2^9;
C4[9] = + 1/737176440 * k2^9;
Back to Transverse Mercator Projection. Forward to Geocentric coordinates. Up to Contents.