GeographicLib: Geodesics on the ellipsoid

Back to Magnetic models. Forward to Transverse Mercator projection. 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 GeodesicProj utility provides an interface to these projections.

References

F. W. Bessel, The calculation of longitude and latitude from geodesic measurements (1825), Astron. Nachr. 331(8), 852-861 (2010); translated by C. F. F. Karney and R. E. Deakin. Preprint: arXiv:0908.1824.
F. R. Helmert, Mathematical and Physical Theories of Higher Geodesy, Part 1 (1880), Aeronautical Chart and Information Center (St. Louis, 1964), Chaps. 5–7.
J. Danielsen, The Area under the Geodesic, Survey Review 30 (232), 61–66 (1989).
C. F. F. Karney, Geodesics on an ellipsoid of revolution, Feb. 2011; preprint arxiv:1102.1215v1; resource page geod.html.
C. F. F. Karney, Algorithms for geodesics, Sept. 2011, preprint arXiv:1109.4448.
A collection of some papers on geodesics is available at http://geographiclib.sourceforge.net/geodesic-papers/biblio.html

Test data for geodesics

A test set a geodesics is available at

GeodTest.dat.gz

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

GeodTest-short.dat.gz

Each line of the test set gives 10 space delimited numbers

latitude for point 1, lat1 (degrees, exact)
longitude for point 1, lon1 (degrees, always 0)
azimuth for point 1, azi1 (clockwise from north in degrees, exact)
latitude for point 2, lat2 (degrees, accurate to 10^-18 deg)
longitude for point 2, lon2 (degrees, accurate to 10^-18 deg)
azimuth for point 2, azi2 (degrees, accurate to 10^-18 deg)
geodesic distance from point 1 to point 2, s12 (meters, exact)
arc distance on the auxiliary sphere, a12 (degrees, accurate to 10^-18 deg)
reduced length of the geodesic, m12 (meters, accurate to 0.1 pm)
the area under the geodesic, S12 (m², accurate to 1 mm²)

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 <= 180^o, 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 [0^o, 90^o], lon1 is taken to be zero, azi1 is chosen in [0^o, 180^o]. 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 [0^o, 180^o] and azi2 in [0^o, 180^o].

The direct calculation uses an expansion of the geodesic equations accurate to f³⁰ (approximately 1 part in 10⁵⁰) 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:

100000 entries randomly distributed
50000 entries which are nearly antipodal
50000 entries with short distances
50000 entries with one end near a pole
50000 entries with both ends near opposite poles
50000 entries which are nearly meridional
50000 entries which are nearly equatorial
50000 entries running between vertices (azi1 = azi2 = 90^o)
50000 entries ending close to vertices

(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

Direct from point 1:
```
  gunzip -c GeodTest.dat.gz | cut -d' ' -f1,2,3,7 | ./Geod
```
This should yield columns 4, 5, 6, and 9 of the test set.
Direct from point 2:
```
  gunzip -c GeodTest.dat.gz | cut -d' ' -f4,5,6,7 |
  sed "s/ $[^ ]*$$/ -\1/" | ./Geod
```
(The sed command negates the distance.) This should yield columns 1, 2, and 3, and the negative of column 9 of the test set.
Inverse between points 1 and 2:
```
  gunzip -c GeodTest.dat.gz | cut -d' ' -f1,2,4,5 | ./Geod -i
```
This should yield columns 3, 6, 7, and 9 of the test set.

Add, e.g., "-p 6", to the call to Geod to change the precision of the output. Adding "-f" causes Geod to print 12 fields specifying the geodesic; these include the 10 fields in the test set plus the geodesic scales M12 and M21 which are inserted between m12 and S12.

Expansions for geodesics

We give here the series expansions for the various geodesic integrals valid to order f¹⁰. 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 f³⁰ 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 Maxima code, geod.mac.

In the expansions below, we have

alpha is the azimuth
alpha₀ is the azimuth at the equator crossing
lambda is the longitude measured from the equator crossing
omega is the spherical longitude
sigma is the spherical arc length
a is the equatorial radius
b is the polar semi-axis
f is the flattening
e² = f (2 - f)
e'² = e²/(1 - e²)
k² = e'² cos² alpha₀ = 4 eps / (1 - eps)²
n = f / (2 - f)
c² = a²/2 + b²/2 (tanh^-1 e)/e
ep2 = e'²
k2 = k²

The formula for distance is

s/b = I1(sigma)

where

I1(sigma) = A1 (sigma + B1(sigma))
B1(sigma) = sum_{j = 1} C1_j 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 + sum_{j = 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 + k² sin²sigma₂) cos sigma₁ sin sigma₂
- sqrt(1 + k² sin²sigma₁) sin sigma₁ cos sigma₂
- cos sigma₁ cos sigma₂ (J(sigma₂) - J(sigma₁))

where

J(sigma) = I1(sigma) - I2(sigma)
I2(sigma) = A2 (sigma + B2(sigma))
B2(sigma) = sum_{j = 1} C2_j 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 alpha₀ I3(sigma)

where

I3(sigma) = A3 (sigma + B3(sigma))
B3(sigma) = sum_{j = 1} C3_j 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.1215v1 in terms of S,

S = c² alpha + e² a² cos alpha₀ sin alpha₀ I4(sigma)

where

I4(sigma) = sum_{j = 0} C4_j 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 Magnetic models. Forward to Transverse Mercator projection. Up to Contents.