Algorithms for geodesics

Charles F. F. Karney

Source-linked AI summary

Algorithms for geodesics

Charles F. F. Karney

arXiv:1109.4448v2Sep 21, 2011physics.geo-phphysics.comp-ph

TL;DR

The paper adapts classical geodesic methods for modern computers, addressing accuracy and convergence limitations in direct and inverse problems. It provides a complete inverse solution, differential and integral geodesic properties, and an ellipsoidal gnomonic projection, with accurate, robust, and fast computations.

Problem
Existing geodesic methods needed computer-precision accuracy, an inverse solution converging for all point pairs, and computation of differential and integral properties.
Method
The paper develops auxiliary-sphere series expansions, geodesic differential properties, Newton-based inverse solving, geodesic-area formulas, and an ellipsoidal gnomonic projection.
Results
The algorithms provide accurate, robust, and fast direct and inverse geodesic solutions while computing differential and integral quantities; projection errors remain small over tested radii.
Takeaways & Limitations
The complete inverse solution, new ellipsoidal gnomonic projection, and unified library support geodesic projections and other ellipsoidal geometry problems.
Takeaways & Limitations
The Newton starting procedure is inadequate for nearly antipodal points, where the azimuth depends very sensitively on the spherical separation.

Abstract

from arXiv · show

Algorithms for the computation of geodesics on an ellipsoid of revolution are given. These provide accurate, robust, and fast solutions to the direct and inverse geodesic problems and they allow differential and integral properties of geodesics to be computed.

1. INTRODUCTION

The paper modernizes classical algorithms for geodesics on an ellipsoid of revolution, targeting accurate direct and inverse solutions and additional differential, integral, and projection capabilities.

Geodesic algorithms solve direct problems from a starting point, initial azimuth, and length, and inverse problems between two given points.
The paper adapts Helmert’s methods to modern computers while retaining enough series terms to reach standard machine precision at little computational cost.
Its inverse solution converges for all point pairs, addressing Vincenty’s failure for nearly antipodal points.
Differential geodesic properties determine the behavior of nearby geodesics and enable projection scales without numerical differentiation.
Integral properties provide a method for computing the area of a geodesic polygon, while an ellipsoidal gnomonic projection makes geodesics nearly straight.

2. BASIC EQUATIONS AND DIRECT PROBLEM

The direct problem is formulated through an auxiliary sphere, where ellipsoidal geodesics correspond exactly to spherical great circles and distance and longitude are evaluated with series expansions.

An ellipsoid of revolution is parameterized by equatorial radius a, polar semi-axis b, flattening f, third flattening n, and eccentricities e and e′.
Clairaut’s relation sin α0 = sin α1 cos β1 = sin α2 cos β2 links the azimuths and reduced latitudes along a geodesic.
The auxiliary sphere replaces latitude φ with reduced latitude β while preserving azimuths α, giving an exact correspondence between geodesics and great circles.
The distance integral is expanded as a Fourier series using ε as the expansion parameter, producing half as many terms as corresponding expansions in k2.
The longitude integral is likewise expanded jointly in n and ε, both O(f), with coefficients inserted into established geodesic formulas.
The direct solution computes the initial spherical quantities, evaluates s1 and λ1, advances to s2 = s1 + s12, and recovers the endpoint coordinates and azimuth.

3. DIFFERENTIAL QUANTITIES

The paper develops differential quantities that describe how nearby geodesics separate and uses them to generalize geometric constructions on curved surfaces.

Differential quantities: The paper derives addition rules and series expansions for m12 and M12, including the spherical-limit relations m12 = a sin(s12/a) and M12 = cos(s12/a).
Differential quantities: A nearby-geodesic separation t(s) satisfies a second-order linear homogeneous differential equation governed by the surface’s Gaussian curvature K(s).
Differential quantities: The reduced length m12 measures separation at the endpoint between geodesics crossing initially at a small angle, while M12 measures separation between initially parallel geodesics.
Differential quantities: The reduced length obeys reciprocity, m21 + m12 = 0, and the Wronskian remains constant along the geodesic.
Differential quantities: Geodesic circles and parallels extend plane-geometry concepts to surfaces, but indented reference curves produce self-intersecting parallels beyond a sufficiently large distance.

4. INVERSE PROBLEM

The inverse problem is formulated as a root-finding problem because the ellipsoid’s longitude difference depends on an unknown equatorial azimuth. The paper solves it with a hybrid problem and Newton updates, including special handling for exceptional cases.

4. INVERSE PROBLEM: The inverse problem is intrinsically harder than the direct problem because λ12 depends on the unknown equatorial azimuth α0, making root-finding unavoidable.
4. INVERSE PROBLEM: The method assumes α1, solves for λ12 at the first intersection with latitude φ2, and adjusts α1 with Newton’s method until the specified longitude difference is matched.
4. INVERSE PROBLEM: For fixed endpoint latitudes, λ12 varies with α1 as shown in Fig. 4, and dλ12/dα1 is obtained from the endpoint displacement m12 dα1.
4. INVERSE PROBLEM: Meridional and equatorial geodesics are treated separately because their initial azimuths are known; the general case uses Newton’s method.
4. INVERSE PROBLEM: The hybrid solution computes reduced latitudes, azimuths, spherical quantities, and finally λ12 using the direct-problem solution.
4. INVERSE PROBLEM: When β2 = ±β1 and α1 = 1/2π, the derivative expression becomes indeterminate because both m12 and cos α2 vanish, so a limiting procedure is required.
4. INVERSE PROBLEM: Compared with Vincenty’s modified method, the described Newton method requires only a few iterations, whereas the modification can require many thousands for nearly antipodal points.

5. STARTING POINT FOR NEWTON’S METHOD

The paper constructs starting guesses for Newton’s method, using ordinary approximations for most cases and an astroid-based construction near antipodal configurations.

5. STARTING POINT FOR NEWTON’S METHOD: In most cases, a starting guess for α1 is obtained from the auxiliary-sphere formulation before Newton refinement.
5. STARTING POINT FOR NEWTON’S METHOD: For the short sample in Table 3, assuming ω12 = λ12/ w̄ gives negligible error, so the inverse solution is completed using s12 = a w̄σ12.
5. STARTING POINT FOR NEWTON’S METHOD: The ordinary auxiliary-sphere procedure is inadequate near antipodal points because α1 depends very sensitively on ω12 when both components of z1 are small.
5. STARTING POINT FOR NEWTON’S METHOD: On an ellipsoid, geodesics from A form an astroid-shaped envelope centered at the antipodal point, with extent O(f).
5. STARTING POINT FOR NEWTON’S METHOD: Given the scaled coordinates of B, the astroid construction reduces the starting-angle problem to a quartic whose positive root corresponds to the shortest path.
5. STARTING POINT FOR NEWTON’S METHOD: Near the antipodal point, geodesics are approximated by a family of lines whose envelope has parametric coordinates x = −sin^3 α1 and y = −cos^3 α1.
5. STARTING POINT FOR NEWTON’S METHOD: Tables 4–6 illustrate the complete inverse-problem solution for nearly antipodal points, while Table 3 provides a non-antipodal sample.

6. AREA

The paper extends geodesic-area computation to higher order, recasting the result as a trigonometric sum accurate to round-off and applicable to geodesic polygons.

6. AREA: The area algorithm extends Danielsen’s method to higher order and recasts its series as a simple trigonometric sum.The resulting computation is designed to achieve round-off accuracy.
6. AREA: Geodesic quadrilateral area is expressed as a spherical term plus an integral representing the ellipsoidal correction.
6. AREA: Tables 5–7 provide sample inverse-problem and geodesic-area calculations, including an area bounded by the equator and a specified geodesic.Table 6 reports agreement between the computed and specified longitude differences.
6. AREA: Expanding the integrand in e′2 and k2 and integrating produces the correction series used in the area calculation.
6. AREA: Summing edge contributions gives the area of a geodesic polygon, with 2πc2 added when the polygon encircles a pole.The spherical term contributes c2(α2 − α1) and requires accurate evaluation for short edges.

7. IMPLEMENTATION

The algorithms can be converted directly into working code using polynomial expansions, numerical-stability practices, and implementation guidance validated on WGS84 geodesics.

7. IMPLEMENTATION: The polynomial expansions produce final results accurate to O(f 6), with truncation error below IEEE double-precision round-off for f = 1/150.
7. IMPLEMENTATION: Horner evaluation improves speed and minimizes round-off errors, while flattening-dependent expressions can be precomputed.
7. IMPLEMENTATION: Implementation details include Newton starting rules, an improved initial guess, convergence criteria, auxiliary-sphere trigonometry, and intermediate-point computation.
7. IMPLEMENTATION: The C++ implementation in GeographicLib was tested against a large set of geodesics for the WGS84 ellipsoid.

8. ELLIPSOIDAL GNOMONIC PROJECTION

The paper derives an ellipsoidal generalization of the gnomonic projection in which geodesics are nearly straight, then uses it for accurate intersection and interception algorithms.

8. ELLIPSOIDAL GNOMONIC PROJECTION: Because an ellipsoid lacks constant Gaussian curvature, its geodesics cannot all map exactly to straight lines; the derived projection makes them very nearly straight.
8. ELLIPSOIDAL GNOMONIC PROJECTION: The generalized projection is constructed from a center A by solving the inverse geodesic problem to each point B and using the resulting geodesic quantities.
8. ELLIPSOIDAL GNOMONIC PROJECTION: The projection coordinates are x = ρ sin α1 and y = ρ cos α1, with radial coordinate ρ = m12/M12.
8. ELLIPSOIDAL GNOMONIC PROJECTION: The projection is undefined when M12 ≤ 0, and its reverse uses the projected azimuth, Newton’s method, and a direct geodesic solution.
8. ELLIPSOIDAL GNOMONIC PROJECTION: The derived projection has h/r proportional to r3/a3, outperforming alternative generalizations with h/r proportional to r/a or r2/a2.
8. ELLIPSOIDAL GNOMONIC PROJECTION: Within 1000-km and 2000-km projection circles, straight-line estimates deviate from true geodesics by less than 1.7 m and 28 m, respectively.The corresponding maximum endpoint azimuth errors are 1.1′′ and 8.6′′, while maximum length errors are 5.4 µm and 730 µm.
8. ELLIPSOIDAL GNOMONIC PROJECTION: The projection supports intersection and interception problems by iteratively projecting, solving a planar construction, and updating the geographic estimate.When points lie within about a quarter meridian, both algorithms converge quadratically to the exact result.

9. CONCLUSIONS

The paper provides accurate, robust, and fast solutions to the direct and inverse geodesic problems, including the first complete inverse solution. It also introduces an ellipsoidal gnomonic projection and packages geodesic capabilities into a single library for broader geometrical applications.

The algorithms solve the direct and inverse geodesic problems while computing differential and integral quantities m12, M12, M21, and S12.
The paper presents the first complete solution to the inverse geodesic problem.
The ellipsoidal gnomonic projection is introduced as a new tool for solving geometrical problems on the ellipsoid.
Packaging the geodesic capabilities into a single library supports projections, ellipsoidal trigonometry, and maritime-boundary problems.The library provides projection scales in terms of m12 and M12.