Algorithms for geodesics

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.


INTRODUCTION
The shortest path between two points on the earth, customarily treated as an ellipsoid of revolution, is called a geodesic.Two geodesic problems are usually considered: the direct problem of finding the end point of a geodesic given its starting point, initial azimuth, and length; and the inverse problem of finding the shortest path between two given points.Referring to Fig. 1, it can be seen that each problem is equivalent to solving the geodesic triangle NAB given two sides and their included angle (the azimuth at the first point, α 1 , in the case of the direct problem and the longitude difference, λ 12 , in the case of the inverse problem).The framework for solving these problems was laid down by Legendre (1806), Oriani (1806Oriani ( , 1808Oriani ( , 1810)), Bessel (1825), and Helmert (1880).Based on these works, Vincenty (1975a) devised algorithms for solving the geodesic problems suitable for early programmable desk 1 The ellipsoidal triangle NAB.N is the north pole, NAF and NBH are meridians, and AB is a geodesic of length s12.The longitude of B relative to A is λ12; the latitudes of A and B are φ1 and φ2.EFH is the equator with E also lying on the extension of the geodesic AB; and α0, α1, and α2 are the azimuths (in the forward direction) of the geodesic at E, A, and B.
calculators; these algorithms are in widespread use today.A good summary of Vincenty's algorithms and the earlier work in the field is given by Rapp (1993, Chap. 1).
The goal of this paper is to adapt the geodesic methods of Helmert (1880) and his predecessors to modern computers.The current work goes beyond Vincenty in three ways: (1) The accuracy is increased to match the standard precision of most computers.This is a relatively straightforward task of retaining sufficient terms in the series expansions and can be achieved at little computational cost.(2) A solution of the inverse problem is given which converges for all pairs of points.(Vincenty's method fails to converge for nearly antipodal points.)(3) Differential and integral properties of the geodesics are computed.The differential properties allow the behavior of nearby geodesics to be determined, which enables the scales of geodesic projections to be computed without resorting to numerical differentiation; crucially, one of the differential quantities is also used in the solution of the inverse problem.The integral properties provide a method for finding the area of a geodesic polygon, extending the work of Danielsen (1989).
Section 2 reviews the classical solution of geodesic problem by means of the auxiliary sphere and provides expansions of the resulting integrals accurate to O(f 6 ) (where f is the flattening of the ellipsoid).These expansions can be inserted into the solution for the direct geodesic problem presented by, for example, Rapp (1993) to provide accuracy to machine precision.Section 3 gives the differential properties of geodesics reviewing the results of Helmert (1880) for the reduced length and geodesic scale and give the key properties of these quantities and appropriate series expansions to allow them to be calculated accurately.Knowledge of the reduced length enables the solution of the inverse problem by Newton's method which is described in Sect. 4. Newton's method requires a good starting guess and, in the case of nearly antipodal points, this is provided by an approximate solution of the inverse problem by Helmert (1880), as given in Sect. 5.The computation of area between a geodesic and the equator is formulated in Sect.6, extending the work of Danielsen (1989).Some details of the implementation and present accuracy and timing data are discussed in Sect.7. As an illustration of the use of these algorithms, Sect.8 gives an ellipsoidal gnomonic projection in which geodesics are very nearly straight.This provides a convenient way of solving several geodesic problems.
For the purposes of this paper, it is useful to generalize the definition of a geodesic.The geodesic curvature, κ, of an arbitrary curve at a point P on a surface is defined as the curvature of the projection of the curve onto a plane tangent to the surface at P .All shortest paths on a surface are straight, defined as κ = 0 at every point on the path.In the rest of this paper, I use straightness as the defining property of geodesics; this allows geodesic lines to be extended indefinitely (beyond the point at which they cease to be shortest paths).
Several of the results reported here appeared earlier in a technical report, Karney (2011).

BASIC EQUATIONS AND DIRECT PROBLEM
I consider an ellipsoid of revolution with equatorial radius a, and polar semi-axis b, flattening f , third flattening n, eccentricity e, and second eccentricity e ′ given by As a consequence of the rotational symmetry of the ellipsoid, geodesics obey a relation found by Clairaut (1735), namely where β is the reduced latitude (sometimes called the parametric latitude), given by The geodesic problems are most easily solved by using an auxiliary sphere which allows an exact correspondence to be made between a geodesic and a great circle on a sphere.On the sphere, the latitude φ is replaced by the reduced latitude β, and azimuths α are preserved.From Fig. 2, it is clear that Clairaut's equation, sin α 0 = sin α cos β, is just the sine rule applied to the sides NE and NP of the triangle NEP and their opposite angles.The third side, the spherical arc length σ, and its opposite angle, the spherical longitude ω, are related to the equivalent quantities on the ellipsoid, the distance s and longitude λ, by (Rapp, 1993, Eqs. (1.28) and (1.170)) where FIG. 2 The elementary ellipsoidal triangle NEP mapped to the auxiliary sphere.NE and NP G are meridians; EG is the equator; and EP is the great circle (i.e., the geodesic).The corresponding ellipsoidal variables are shown in parentheses.Here P represents an arbitrary point on the geodesic EAB in Fig. 1.
See also Eqs. (5.4.9) and (5.8.8) of Helmert (1880).The origin for s, σ, λ, and ω is the point E, at which the geodesic crosses the equator in the northward direction, with azimuth α 0 .The point P can stand for either end of the geodesic AB in Fig. 1, with the quantities β, α, σ, ω, s, and λ acquiring a subscript 1 or 2. I also define s 12 = s 2 − s 1 as the length of AB, with λ 12 , σ 12 , and ω 12 defined similarly.(In this paper, α 2 is the forward azimuth at B. Several authors use the back azimuth instead; this is given by α 2 ± π.) Because Eqs. ( 7) and ( 8) depend on α 0 , the mapping between the ellipsoid and the auxiliary sphere is not a global mapping of one surface to another; rather the auxiliary sphere should merely be regarded as a useful mathematical technique for solving geodesic problems.Similarly, because the origin for λ depends on the geodesic, only longitude differences, e.g., λ 12 , should be used in converting between longitudes relative to the prime meridian and λ.
In solving the spherical trigonometrical problems, the following equations relating the sides and angles of NEP are useful, where i = √ −1 and ph(x + iy) is the phase of a complex number (Olver et al., 2010, §1.9(i)), typically given by the library function atan2(y, x).Equation (10) merely recasts Eq. ( 5) in a form that allows it to be evaluated accurately when α 0 is close to 1 2 π.The other relations are obtained by applying Napier's rules of circular parts to NEP .
The distance integral, Eq. ( 7), can be expanded in a Fourier series with the coefficients determined by expanding the integral in a Taylor series.It is advantageous to follow Bessel (1825, §5) and Helmert (1880, Eq. (5.5.1)) and use ǫ, defined by as the expansion parameter.This leads to expansions with half as many terms as the corresponding ones in k 2 .The expansion can be conveniently carried out to arbitrary order by a computer algebra system such as Maxima (2009) which yields This extends Eq. (5.5.7) of Helmert (1880) to higher order.These coefficients may be inserted into Eq.(1.40) of Rapp (1993) using where here, and subsequently in Eqs. ( 22) and ( 26), a script letter, e.g., B, is used to stand for Rapp's coefficients.

DIFFERENTIAL QUANTITIES
Before turning to the inverse problem, I present Gauss' solution for the differential behavior of geodesics.One differential quantity, the reduced length m 12 , is needed in the solution of the inverse problem by Newton's method (Sect.4) and an expression for this quantity is given at the end of this section.However, because this and other differential quantities aid in the solution of many geodesic problems, I also discuss their derivation and present some of their properties.
Consider a reference geodesic parametrized by distance s and a nearby geodesic separated from the reference by infinitesimal distance t(s).Gauss (1828) showed that t(s) satisfies the differential equation A geometric proof of Eq. ( 29) is shown in (c); here AB and A ′ B ′ are parallel at B and B ′ , BAB ′ = dα1, BB ′ = m12 dα1, AA ′ = M21m12 dα1, and finally AB ′ A ′ = M21 dα1, from which Eq. ( 29) follows.
where K(s) is the Gaussian curvature of the surface.As a second order, linear, homogeneous differential equation, its solution can be written as where A and B are (infinitesimal) constants and t A and t B are independent solutions.When considering the geodesic segment spanning and to write The quantity m 12 is the reduced length of the geodesic (Christoffel, 1868).Consider two geodesics which cross at s = s 1 at a small angle dα Several relations between m 12 and M 12 follow from the defining equation, Eq. ( 27).The reduced length obeys a reciprocity relation (Christoffel, 1868, §9), m 21 + m 12 = 0; the Wronskian is given by and the derivatives are The constancy of the Wronskian follows by noting that its derivative with respect to s 2 vanishes; its value is found by evaluating it at s 2 = s 1 .A geometric proof of Eq. ( 29) is given in Fig 3(c) and Eq. ( 30) then follows from Eq. ( 28).
With knowledge of the derivatives, addition rules for m 12 and M 12 are easily found, where points 1, 2, and 3 all lie on the same geodesic.Geodesics allow concepts from plane geometry to be generalized to apply to a curved surface.In particular, a geodesic circle may be defined as the curve which is a constant geodesic distance from a fixed point.Similarly, a geodesic parallel to a reference curve is the curve which is a constant geodesic distance from that curve.(Thus a circle is a special case of a parallel obtained in the limit when the reference curve degenerates to a point.)Parallels occur naturally when considering, for example, the "12-mile limit" for territorial waters which is the boundary of points lying within 12 nautical miles of a coastal state.
The geodesic curvature of a parallel can be expressed in terms of m 12 and M 12 .Let point 1 be an arbitrary point on the reference curve with geodesic curvature κ 1 .Point 2 is the corresponding point on the parallel, a fixed distance s 12 away.The geodesic curvature of the parallel at that point is found from Eqs. ( 29) and ( 30), The curvature of a circle is given by the limit κ 1 → ∞, If the reference curve is a geodesic (κ 1 → 0), then the curvature of its parallel is If the reference curve is indented, then the parallel intersects itself at a sufficiently large distance from the reference curve.This begins to happen when κ 2 → ∞ in Eq. ( 34).
The results above apply to general surfaces.For a geodesic on an ellipsoid of revolution, the Gaussian curvature of the surface is given by Helmert (1880, Eq. (6.5.1)) solves Eq. ( 27) in this case to give where Equation ( 39) may be obtained from Eq. (6.9.7) of Helmert ( 1880), which gives dm 12 /ds 2 ; M 12 may then be found from Eq. ( 29) with an interchange of indices.In the spherical limit, f → 0, Eqs. ( 38) and ( 39) reduce to m 12 = a sin σ 12 = a sin(s 12 /a), M 12 = cos σ 12 = cos(s 12 /a).
The integral I 2 (σ) in Eq. ( 40) may be expanded in a Fourier series in similar fashion to I 1 (σ), Eq. ( 15), where

INVERSE PROBLEM
The inverse problem is intrinsically more complicated than the direct problem because the given included angle, λ 12 in Fig. 1, is related to the corresponding angle on the auxiliary sphere ω 12 via an unknown equatorial azimuth α 0 .Thus, the inverse problem inevitably becomes a root-finding exercise.
I tackle this problem as follows.Assume that α 1 is known.Solve the hybrid geodesic problem: given φ 1 , φ 2 , and α 1 , find λ 12 corresponding to the first intersection of the geodesic with the circle of latitude φ 2 .The resulting λ 12 differs, in general, from the given λ 12 ; so adjust α 1 using Newton's method until the correct λ 12 is obtained.
I begin by putting the points in a canonical configuration, This may be accomplished swapping the end points and the signs of the coordinates if necessary, and the solution may similarly be transformed to apply to the original points.All geodesics with α 1 ∈ [0, π] intersect latitude φ 2 with λ 12 ∈ [0, π].Furthermore, the search for solutions can be restricted to α 2 ∈ [0, 1 2 π], because this corresponds to the first intersection with latitude φ 2 .
The behavior of λ 12 as a function of α 1 is shown in Fig. 4.
To apply Newton's method, an expression for dλ 12 /dα 1 is needed.Consider a geodesic with initial azimuth α 1 .If the azimuth is increased to α 1 + dα 1 with its length held fixed, then the other end of the geodesic moves by m 12 dα 1 in a direction 1 2 π + α 2 .If the geodesic is extended to intersect the parallel φ 2 once more, the point of intersection moves by m 12 dα 1 / cos α 2 ; see Fig. 5.The radius of this parallel is a cos β 2 ; thus the rate of change of the longitude difference is This equation can also be obtained from Eq. (6.9.8b) of Helmert (1880).Equation ( 46) becomes indeterminate when β 2 = ±β 1 and α 1 = 1 2 π, because m 12 and cos α 2 both vanish.In this case, it is necessary to let α 1 = 1 2 π + δ and to take the limit δ → ±0, which gives where sign(cos α 1 ) = − sign(δ).A numerical example of solving the inverse geodesic problem by this method is given at the end of the next section.Vincenty (1975a), who uses the iterative method of Helmert (1880, §5.13) to solve the inverse problem, was aware of its failure to converge for nearly antipodal points.In an unpublished report (Vincenty, 1975b), he gives a modification of his method which deals with this case.Unfortunately, this sometimes requires many thousands of iterations to converge, whereas Newton's method as described here only requires a few iterations.

STARTING POINT FOR NEWTON'S METHOD
To complete the solution of the inverse problem a good starting guess for α 1 is needed.In most cases, this is provided TABLE 3 First sample inverse calculation specified by φ1 = −30.12345 • , φ2 = −30.12344 • , and λ12 = 0.000 05 • .Because the points are not nearly antipodal, an initial guess for α1 is found assuming ω12 = λ12/ w.However, in this case, the line is short enough that the error in ω12 is negligible at the precision given and the solution of the inverse problem is completed by using s12 = a wσ12.More generally, the value of α1 would be refined using Newton's method.by assuming that ω 12 = λ 12 / w, where

Qty
and solving for the great circle on the auxiliary sphere, using (Vincenty, 1975a) An example of the solution of the inverse problem by this method is given in Table 3.This procedure is inadequate for nearly antipodal points because both the real and imaginary components of z 1 are small and α 1 depends very sensitively on ω 12 .In the corresponding situation on the sphere, it is possible to determine α 1 by noting that all great circles emanating from A meet at O, the point antipodal to A. Thus α 1 may be determined as the supplement of the azimuth of the great circle BO at O; in addition, because B and O are close, it is possible to approximate the sphere, locally, as a plane.
The situation for an ellipsoid is slightly different because the geodesics emanating from A, instead of meeting at a point, form an envelope, centered at O, in the shape of an astroid whose extent is O(f ) (Jacobi, 1891, Eqs. ( 16)-( 17)).The position at which a particular geodesic touches this envelope is given by the condition m 12 = 0. However elementary methods can be used to determine the envelope.Consider a geodesic leaving A (with . This first intersects the circle of opposite latitude, β 2 = −β 1 , with σ 12 = ω 12 = π and α 2 = π − α 1 .Equation ( 8) then gives Define a plane coordinate system (x, y) centered on the antipodal point where ∆ = f aπ cos 2 β 1 is the unit of length, i.e., In this coordinate system, Eq. ( 52) corresponds to the point x = − sin α 1 , y = 0 and the slope of the geodesic is − cot α 1 .Thus, in the neighborhood of the antipodal point, the geodesic may be approximated by where terms of order f 2 have been neglected.Allowing α 1 to vary, Eq. ( 54) defines a family of lines approximating the geodesics emanating from A. Differentiating this equation with respect to α 1 and solving the resulting pair of equations for x and y gives the parametric equations for the astroid, x = − sin 3 α 1 , and y = − cos 3 α 1 .Note that, for the ordering of points given by Eq. ( 44), x ≤ 0 and y ≤ 0. Given x and y (i.e., the position of point B), Eq. ( 54) may be solved to obtain a first approximation to α 1 .This prescription is given by Helmert (1880, Eq. (7.3.7)) who notes that this results in a quartic which may be found using the construction given in Fig. 6.Here COD and BED are similar triangles; if the (signed) length BC is µ, then an equation for µ can be TABLE 4 Second sample inverse calculation specified by φ1 = −30 • , φ2 = 29.9• , and λ12 = 179.8• .Because the points are nearly antipodal, an initial guess for α1 is found by solving the astroid problem.Here µ is the positive root of Eq. ( 55).If y = 0, then α1 is given by Eq. ( 57).The value of α1 is used in Table 5 found by applying Pythagoras' theorem to COD, which can be expanded to give a 4th-order polynomial in µ, Descartes' rule of signs shows that, for y = 0, there is one positive root (Olver et al., 2010, §1.11(ii)) and this is the solution corresponding to the shortest path.This root can be found by standard methods (Olver et al., 2010, §1.11(iii)).Equation (55) arises in converting from geocentric to geodetic coordinates, and I use the solution to that problem given by Vermeille (2002).The azimuth can then be determined from the triangle COD in Fig. 6, If y = 0, the solution is found by taking the limit y → 0, Tables 4-6 together illustrate the complete solution of the inverse problem for nearly antipodal points.

AREA
The last geodesic algorithm I present is for geodesic areas.Here, I extend the method of Danielsen (1989) to higher order so that the result is accurate to round-off, and I recast his series into a simple trigonometric sum.
Let S 12 be the area of the geodesic quadrilateral AFHB in Fig. 1.Following Danielsen (1989), this can be expressed as the sum of a spherical term and an integral giving the ellipsoidal correction, where is the authalic radius, Expanding the integrand in powers of e ′2 and k 2 and performing the integral gives TABLE 7 The calculation of the area between the equator and the geodesic specified by φ1 = 40 An example of the computation of S 12 is given in Table 7. Summing S 12 , Eq. ( 58), over the edges of a geodesic polygon gives the area of the polygon provided that it does not encircle a pole; if it does, 2πc 2 should be added to the result.The first term in Eq. ( 59) contributes c 2 (α 2 − α 1 ) to S 12 .This is the area of the quadrilateral AFHB on a sphere of radius c and it is proportional to its spherical excess, α 2 − α 1 , the sum of its interior angles less 2π.It is important that this term be computed accurately when the edge is short (and α 1 and α 2 are nearly equal).A suitable identity for α 2 − α 1 is given by Bessel (1825, §11),

IMPLEMENTATION
The algorithms described in the preceding sections can be readily converted into working code.The polynomial expansions, Eqs. ( 17), ( 18), ( 21), ( 24), ( 25), ( 42), (43), and ( 63), are such that the final results are accurate to O(f 6 ) which means that, even for f = 1 150 , the truncation error is smaller than the round-off error when using IEEE double precision arithmetic (with the fraction of the floating point number represented by 53 bits).For speed and to minimize round-off errors, the polynomials should be evaluated with the Horner method.The parenthetical expressions in Eqs. ( 24), (25), and (63) depend only on the flattening of the ellipsoid and can be computed once this is known.When determining many points along a single geodesic, the polynomials need be evaluated just once.Clenshaw (1955) summation should be used to sum the Fourier series, Eqs. ( 15), ( 23), (41), and (62).
There are several other details to be dealt with in implementing the algorithms: where to apply the two rules for choosing starting points for Newton's method, a slight improvement to the starting guess Eq. ( 56), the convergence cri-terion for Newton's method, how to minimize round-off errors in solving the trigonometry problems on the auxiliary sphere, rapidly computing intermediate points on a geodesic by using σ 12 as the metric, etc.I refer the reader to the implementations of the algorithms in GeographicLib (Karney, 2012) for possible ways to address these issues.The C++ implementation has been tested against a large set of geodesics for the WGS84 ellipsoid; this was generated by continuing the series expansions to O(f 30 ) and by solving the direct problem using with high-precision arithmetic.The round-off errors in the direct and inverse methods are less than 15 nanometers and the error in the computation of the area S 12 is about 0.1 m 2 .Typically, 2 to 4 iterations of Newton's method are required for convergence, although in a tiny fraction of cases up to 16 iterations are required.No convergence failures are observed.With the C++ implementation compiled with the g++ compiler, version 4.4.4,and running on a 2.66 GHz Intel processor, solving the direct geodesic problem takes 0.88 µs, while the inverse problem takes 2.34 µs (on average).Several points along a geodesic can be computed at the rate of 0.37 µs per point.These times are comparable to those for Vincenty's algorithms implemented in C++ and run on the same architecture: 1.11 µs for the direct problem and 1.34 µs for the inverse problem.(But note that Vincenty's algorithms are less accurate than those given here and that his method for the inverse problem sometimes fails to converge.)

ELLIPSOIDAL GNOMONIC PROJECTION
As an application of the differential properties of geodesics, I derive a generalization of the gnomonic projection to the ellipsoid.The gnomonic projection of the sphere has the property that all geodesics on the sphere map to straight lines (Snyder, 1987, §22).Such a projection is impossible for an ellipsoid because it does not have constant Gaussian curvature (Beltrami, 1865, §18); nevertheless, a projection can be constructed in which geodesics are very nearly straight.
The spherical gnomonic projection is the limit of the doubly azimuthal projection of the sphere, wherein the bearings from two fixed points A and A ′ to B are preserved, as A ′ approaches A (Bugayevskiy and Snyder, 1995).The construction of the generalized gnomonic projection proceeds in the same way; see Fig. 7. Draw a geodesic A ′ B ′ such that it is parallel to the geodesic AB at A. Its initial separation from AB is sin γ dt; at B ′ , the point closest to B, the separation becomes M 12 sin γ dt (in the limit dt → 0).Thus the difference in the azimuths of the geodesics A ′ B and A ′ B ′ at A ′ is (M 12 /m 12 ) sin γ dt, which gives γ + γ ′ = π − (M 12 /m 12 ) sin γ dt.Now, solving the planar triangle problem with γ and γ ′ as the two base angles gives the distance AB on the projection plane as m 12 /M 12 .
This leads to the following specification for the generalized gnomonic projection.Let the center point be A; for an arbitrary point B, solve the inverse geodesic problem between A and B; then B projects to the point the projection is undefined if M 12 ≤ 0. In the spherical limit, this becomes the standard gnomonic projection, ρ = a tan σ 12 (Snyder, 1987, p. 165).The azimuthal scale is 1/M 12 and the radial scale, found by taking the derivative dρ/ds 12 and using Eq. ( 28), is 1/M 2 12 .The reverse projection is found by computing α 1 = ph(y + ix), finding s 12 using Newton's method with dρ/ds 12 = 1/M 2 12 (i.e., the radial scale), and solving the resulting direct geodesic problem.
In order to gauge the usefulness of the ellipsoidal gnomonic projection, consider two points on the earth B and C, map these points to the projection, and connect them with a straight line in this projection.If this line is mapped back onto the surface of the earth, it will deviate slightly from the geodesic BC.To lowest order, the maximum deviation h occurs at the midpoint of the line segment BC; empirically, I find where l is the length of the geodesic, K is the Gaussian curvature, ∇K is evaluated at the center of the projection A, and t is the perpendicular vector from the center of projection to the geodesic.The deviation in the azimuths at the end points is about 4h/l and the length is greater than the geodesic distance by about 8 3 h 2 /l.In the case of an ellipsoid of revolution, the curvature is given by differentiating Eq. ( 37) with respect to φ and dividing by the meridional radius of curvature to give where φ is a unit vector pointing north.Bounding the magnitude of h, Eq. ( 66), over all the geodesics whose end points lie within a distance r of the center of projection, gives (in the limit that f and r are small) The maximum value is attained when the center of projection is at φ = ±45 • and the geodesic is running in an east-west direction with the end points at bearings ±45 • or ±135 • from the center.
Others have proposed different generalizations of the gnomonic projection.Bowring (1997) and Williams (1997) give a projection in which great ellipses project to straight lines; Letoval'tsev (1963) suggests a projection in which normal sections through the center point map to straight lines.Empirically, I find that h/r is proportional to r/a and r 2 /a 2 for these projections.Thus, neither does as well as the projection derived above (for which h/r is proportional to r 3 /a 3 ) at preserving the straightness of geodesics.
As an illustration of the properties of the ellipsoidal gnomonic projection, Eq. ( 65), consider Fig. 8 in which a projection of Europe is shown.The two circles are geodesic circles of radii 1000 km and 2000 km.If the geodesic between any two points within one of these circles is estimated by using a straight line on this figure, the deviation from the true geodesic is less than 1.7 m and 28 m, respectively.The maximum errors in the end azimuths are 1.1 ′′ and 8.6 ′′ and the maximum errors in the lengths are only 5.4 µm and 730 µm.
The gnomonic projection can be used to solve two geodesic problems accurately and rapidly.The first is the intersection problem: given two geodesics between A and B and between C and D, determine the point of intersection, O.This can be solved as follows.Guess an intersection point O (0)  Provided the given points lie within about a quarter meridian of the intersection or interception points (so that the gnomonic projection is defined), these algorithms converge quadratically to the exact result.

CONCLUSIONS
The classical geodesic problems entail solving the ellipsoidal triangle NAB in Fig. 1, whose sides and angles are represented by φ 1 , φ 2 , s 12 and α 1 , α 2 , λ 12 .In the direct problem φ 1 , α 1 , and s 12 are given, while in the inverse problem φ 1 , λ 12 , and φ 2 are specified; and the goal in each case is to solve for the remaining side and angles.The algorithms given here provide accurate, robust, and fast solutions to these problems; they also allow the differential and integral quantities m 12 , M 12 , M 21 , and S 12 to be computed.
Much of the work described here involves applying standard computational techniques to earlier work.However, at least two aspects are novel: (1) This paper presents the first complete solution to the inverse geodesic problem.(2) The ellipsoidal gnomonic projection is a new tool to solve various geometrical problems on the ellipsoid.
Furthermore, the packaging of these various geodesic capabilities into a single library is also new.This offers a straightforward solution of several interesting problems.Two geodesic projections, the azimuthal equidistant projection and the Cassini-Soldner projection, are simple to write and their domain of applicability is not artificially restricted, as would be the case, for example, if the series expansion for the Cassini-Soldner projection were used (Snyder, 1987, §13); the scales for these projections are simply given in terms of m 12 and M 12 .Several other problems can be readily tackled with this library, e.g., solving other ellipsoidal trigonometry problems and finding the median line and other maritime boundaries.These and other problems are explored in Karney (2011).The web page http://geographiclib.sf.net/geod.htmlprovides additional information, including the Maxima (2009) code used to carry out the Taylor expansions and a JavaScript implementation which allows geodesic problems to be solved on many portable devices.
1 , Fig 3(a); at s = s 2 , they will be separated by a distance m 12 dα 1 .Similarly I call M 12 the geodesic scale.Consider two geodesics which are parallel at s = s 1 and separated by a small distance dt 1 , Fig 3(b); at s = s 2 , they will be separated by a distance M 12 dt 1 .

FIG. 6
FIG.6The solution of the astroid equations by similar triangles.The scaled coordinates of B are (x, y); O is the point antipodal to A. The line BCD, which is given by Eq. (54), is the continuation of the geodesic from AB with C being its intersection with the circle β = −β1 and D its intersection with the meridian λ = λ1 + π.The envelope of lines satisfying CD = 1 gives the astroid, a portion of which is shown by the curves.
FIG.7The construction of the generalized gnomonic projection as the limit of a doubly azimuthal projection.
and use this as the center of the gnomonic projection; define a, b, c, d as the positions of A, B, C, D in the projection; find the intersection of AB and CD in the projection, i.e.,o = (c × d • ẑ)(b − a) − (a × b • ẑ)(d − c) (b − a) × (d − c) • ẑ ,(69)where ˆindicates a unit vector (â = a/a) and ẑ = x × ŷ is in the direction perpendicular to the projection plane.Project o back to geographic coordinates O(1) and use this as a new center of projection; iterate this process untilO (i) = O (i−1)which is then the desired intersection point.The second problem is the interception problem: given a geodesic between A and B, find the point O on the geodesic which is closest to a given point C. The solution is similar to that for the intersection problem; however the interception point in the projection is o = c • (b − a)(b − a) − (a × b • ẑ)ẑ × (b − a) |b − a| 2 .

TABLE 5
Second sample inverse calculation, continued.Here λ

TABLE 6
Second sample inverse calculation, concluded.Here the hybrid problem (φ1, φ2, and α1 given) is solved.The computed value of λ12 matches that given in the specification of the inverse problem in Table4.
(Wessel and Smith, 2010) and North Africa in the ellipsoidal gnomonic projection with center at (45 • N, 12 • E) near Venice.The graticule lines are shown at multiples of 10 • .The two circles are centered on the projection center with (geodesic) radii of 1000 km and 2000 km.The data for the coast lines is taken from GMT(Wessel and Smith, 2010)at "low" resolution.