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
The shortest path between two points on a ellipsoid at (lat1, lon1) and (lat2, lon2) is called the geodesic. Its length is s12 and the geodesic from point 1 to point 2 has azimuths azi1 and azi2 at the two end points. (The azimuth is the heading measured clockwise from north. azi2 is the "forward" azimuth, i.e., the heading that takes you beyond point 2 not back to point 1.) The reduced length of the geodesic, m12, is the rate at which the second point moves as the geodesic is rotated about the first point. Christoffel (1868) proved that this is equal to the rate at which the first point moves when the geodesic is rotated about the second point.
Given lat1, lon1, azi1, and s12, we can determine lat2, lon2, azi2, m12. This is the direct geodesic problem. (If s12 is sufficiently large that the geodesic wraps more than halfway around the earth, there will be another geodesic between the points with a smaller s12.)
Given lat1, lon1, lat2, and lon2, we can determine azi1, azi2, s12, m12. This is the inverse geodesic problem. Usually, the solution to the inverse problem is unique. In cases where there are multiple solutions (all with the same s12, of course), all the solutions can be easily generated once a particular solution is provided.
The framework for solving geodesic problems on the ellipsoid was laid down by Bessel (1826) who showed how the problem may be transferred to an "auxiliary sphere" where the latitude phi has been replaced by the reduced latitude beta where tan(beta) = (1 - f) tan(phi). On this sphere, the geodesic is a great circle and the azimuth is the same as on the ellipsoid. However, the ellipsoidal distance is related to the great circle distance by an integral; and the ellipsoidal longitude is similarly related to the longitude on the auxiliary sphere. In the case of small f, Bessel developed series expansions for the integrals and provided a systematic solution to the direct problem. Helmert (1880) extended Bessel's treatment by giving an expression for the reduced length of a geodesic for the ellipsoid. The "azimuthal scale" of the azimuthal equidistant projection about point 1 is given by s12/m12.
The algorithms included here are based on the work of Bessel (1826) and Helmert (1880). Key points are:
Vincenty (1975) has provided algorithms for the direct and inverse geodesic problems based in the series given by Rainsford (1955) (based, in turn, on Bessel). However Vincenty's methods were designed for use on the desk calculators of that era and are not optimal (in terms of accuracy and robustness) for a general purpose computer. The distinguishing features of the algorithms presented here compared to Vincenty's are:
The series expansions for the integrals are obtained using Maxima. In this release, 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 is evaluated by a combination of Horner's method and Clenshaw summation. For the direct problem, we use a reversion of the series for the distance to compute the spherical arc length and thereby complete the solution without the need to iterate.
The Maxima code, geod.mac, used to derive the series for the distance and longitude integrals. In addition this includes code to format the results for inclusion into C++ code. You will need to have Maxima installed to use this code. The comments at the top of geod.mac illustrate how to run it. It takes about 10s to generate the 8th order series.
For the inverse problem, we rearrange the points so that the point 1 is the one closest to a pole and is in the southern hemisphere and the point 2 is to the north and east of point 1. We consider geodesics leaving point 1 with azimuth alpha1 in [0, pi] and we determine the longitude where these geodesics first intersect the latitude = lat2. The longitude is an increasing function of alpha1 and we use Newton's method to determine the value of alpha1 for which the longitude matches lon2. This then gives the desired geodesic.
This class also returns the arc length a12 on the auxiliary sphere and allows distances to be specified in these terms. See the -a option to Geod.
A test set a geodesics is available at
This is about 34 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 f20 (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
gunzip -c GeodTest.dat.gz | cut -d' ' -f1,2,3,7 | ./Geod
gunzip -c GeodTest.dat.gz | cut -d' ' -f4,5,6,7 | sed "s/ \([^ ]*$\)/ -\1/" | ./Geod
gunzip -c GeodTest.dat.gz | cut -d' ' -f1,2,4,5 | ./Geod -i
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).
We regard all the test data as "exact". For each member of the test set, 6 error distances are computed. (The direct problem is solved from both ends and the maximum of the errors is used.)
The last check verifies that azi1 and azi2 returned by the inverse solution correspond to the same path and is important, e.g., when point 1 and point 2 are near opposite poles. To guard against a canceling error, the calculations in last error estimate use a long double version of the direct geodesic calculation (which is about 2000 times more accurate than the double version). The maximum of the resulting errors in less than 15 nm for the double version and 7 pm for the long double version.
We give here the series expansions for the various geodesic integrals valid to order f8.
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.
The series expanded to order f30 are given in geodseries30.html.
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) / (1 - eps);
C1[1] = - 1/2 * eps
+ 3/16 * eps^3
- 1/32 * eps^5
+ 19/2048 * eps^7;
C1[2] = - 1/16 * eps^2
+ 1/32 * eps^4
- 9/2048 * eps^6
+ 7/4096 * eps^8;
C1[3] = - 1/48 * eps^3
+ 3/256 * eps^5
- 3/2048 * eps^7;
C1[4] = - 5/512 * eps^4
+ 3/512 * eps^6
- 11/16384 * eps^8;
C1[5] = - 7/1280 * eps^5
+ 7/2048 * eps^7;
C1[6] = - 7/2048 * eps^6
+ 9/4096 * eps^8;
C1[7] = - 33/14336 * eps^7;
C1[8] = - 429/262144 * eps^8;
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;
C1'[2] = + 5/16 * eps^2
- 37/96 * eps^4
+ 1335/4096 * eps^6
- 86171/368640 * eps^8;
C1'[3] = + 29/96 * eps^3
- 75/128 * eps^5
+ 2901/4096 * eps^7;
C1'[4] = + 539/1536 * eps^4
- 2391/2560 * eps^6
+ 1082857/737280 * eps^8;
C1'[5] = + 3467/7680 * eps^5
- 28223/18432 * eps^7;
C1'[6] = + 38081/61440 * eps^6
- 733437/286720 * eps^8;
C1'[7] = + 459485/516096 * eps^7;
C1'[8] = + 109167851/82575360 * eps^8;
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) * (1 - eps);
C2[1] = + 1/2 * eps
+ 1/16 * eps^3
+ 1/32 * eps^5
+ 41/2048 * eps^7;
C2[2] = + 3/16 * eps^2
+ 1/32 * eps^4
+ 35/2048 * eps^6
+ 47/4096 * eps^8;
C2[3] = + 5/48 * eps^3
+ 5/256 * eps^5
+ 23/2048 * eps^7;
C2[4] = + 35/512 * eps^4
+ 7/512 * eps^6
+ 133/16384 * eps^8;
C2[5] = + 63/1280 * eps^5
+ 21/2048 * eps^7;
C2[6] = + 77/2048 * eps^6
+ 33/4096 * eps^8;
C2[7] = + 429/14336 * eps^7;
C2[8] = + 6435/262144 * eps^8;
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 = (2 - f) / (2 - del) * (1
- 1/16 * nu^2 * del^3
- 1/64 * nu^4 * del^4
+ 1/64 * nu^2 * del^5
+ (1/64 * nu^2 + 1/64 * nu^4 - 1/128 * nu^6) * del^6
+ (3/256 * nu^2 + 7/1024 * nu^4 + 15/2048 * nu^6) * del^7);
C3[1] = + 1/4 * eps
- 1/8 * eps*del
- (1/16 + 1/64 * nu^2) * eps*del^2
- (1/32 + 1/128 * nu^2) * eps*del^3
- (1/64 - 1/128 * nu^2 + 1/128 * nu^4) * eps*del^4
- (1/128 - 1/128 * nu^2 - 3/1024 * nu^4) * eps*del^5
- (1/256 - 5/1024 * nu^2 - 5/1024 * nu^4 + 61/16384 * nu^6) * eps*del^6;
C3[2] = + 1/16 * eps^2
- 3/64 * eps^2*del
- (1/64 + 1/128 * nu^2) * eps^2*del^2
- 1/256 * eps^2*del^3
+ (5/1024 * nu^2 - 7/2048 * nu^4) * eps^2*del^4
+ (1/1024 + 3/1024 * nu^2 + 37/16384 * nu^4) * eps^2*del^5;
C3[3] = + 5/192 * eps^3
- 3/128 * eps^3*del
- (1/192 + 7/1536 * nu^2) * eps^3*del^2
+ 1/768 * nu^2 * eps^3*del^3
+ (1/1024 + 35/12288 * nu^2 - 91/49152 * nu^4) * eps^3*del^4;
C3[4] = + 7/512 * eps^4
- 7/512 * eps^4*del
- (1/512 + 3/1024 * nu^2) * eps^4*del^2
+ (5/8192 + 23/16384 * nu^2) * eps^4*del^3;
C3[5] = + 21/2560 * eps^5
- 9/1024 * eps^5*del
- (3/4096 + 33/16384 * nu^2) * eps^5*del^2;
C3[6] = + 11/2048 * eps^6
- 99/16384 * eps^6*del;
C3[7] = + 429/114688 * eps^7;
When either point is at a pole, the azimuth is defined by keeping the longitude fixed and writing lat = 90 - eps or -90 + eps and taking the limit eps -> 0 from above.
The following prescription allows you to ascertain when there are multiple solutions to the inverse problem and to generate the complete set.
A rather large fraction of the code in GeographicLib::Geodesic deals merely with the spherical trigonometry required to do great circle calculations. When high accuracy is required, care has to be taken to use stable methods for computing angles. Typically, we represent angles internally with the pair [sin(theta), cos(theta)]. This expands the number of representable angles by about 4; in particular, it allows angles very close to all the cardinal points to be resolved. We use expressions for sin(theta) and cos(theta) which avoid large cancellations. Having a comprehensive test set as described in Test data for geodesics helps to identify problems. This representation allows input arguments which are multiples of 90o to be represented exactly.
The code for the inverse treats equatorial and meridional geodesics specially. All other case are handled by the same code (implementing Newton's method) with tests to select appropriate starting guess in difference regimes. Specifying a zero inverse flattening in Geodesic constructor is equivalent to zero flattening (i.e., a sphere). The method accommodates this case without any special coding.
In order to limits the enumeration of cases, the inverse problem is brought to a canonical form by interchanging the points and changing the signs of the coordinates. The effect of these changes are undone at the end by applying sign changes to the final calls to atan2 to give the azimuths.
The derivative needed for Newton's method for the inverse problem is given in terms of the reduced length. The reduced length can also be used to apply Newton's method to solve problems such as finding the distance between a point and a geodesic line.
Convergence failures with Newton's method for the inverse problem is signaled by negating the distance (and the arc distance and reduced length) and by reversing the azimuths. This convention was chosen as it allows both detection of the problem and subsequent use of the values (assuming that they are approximately correct). It is expected that the method converges for all models of the earth and for any two input points. Please report all cases where this does not happen.
Geographic coordinates which are close to zero are quantized (by adding and subtracting 1/16 deg). This allows certain cases of underflow to be avoided. This limits the accuracy to which coordinates can be specified to 0.7 pm. This should not be a problem since the overall accuracy of the algorithms is only 15 nm.
You can specify a negative flattening to get a prolate ellipsoid. The rule is that b is the radius of the axis of symmetry (the polar radius), that a is the equatorial radius, and that b = a (1 - f). The direct solution needed no modification. The inverse solution required a generalization in the way starting points for Newton's method are chosen; in addition a test are added to the case when the endpoints are on the same meridian to determine if a shorter non-meridional geodesic exists.
Most planets and satellites are approximated by oblate ellipsoids. However, satellites whose rotation frequency is tidally locked to a month may be prolate (as example is Saturn's moon, Tethys). Such satellites are typically considerably more eccentric that the earth and it may be necessary to increase the order of the series used in order to maintain accuracy.
GeodesicLine::Position provides a way of doing a series of direct calculations for a single geodesic. This is about 2.1 times faster than calling Geodesic::Direct. A speed-up by a further factor of 1.6 is obtained with the direct calculation if the distance is expressed in terms of the spherical arc length.
On a 2.6 GHz Intel machine (g++, version 4.3.0, -O3), timings are as follows:
The timing for Geodesic::Inverse is a mean. Individual calls may take a few times longer because of the need to do additional Newton iterations. Mean number of Newton iterations is about 3, but this increases (rarely) to 17 with some values of the arguments.
Future work:
1.6.1