The quantities s, \(\sigma \), \(\lambda \), and \(\omega \) are measured from the node, the point where the geodesic crosses the equator in the northward direction (as in an ascending orbital node). A single subscript refers to a specific point, e.g., \(s_1\) measures the distance from the node to point 1; for s, \(\sigma \), and \(\lambda \), a double subscript denotes a difference, e.g., \(\lambda _{12} = \lambda _2 - \lambda _1\).

1 Introduction

In an earlier paper, Algorithms for geodesics (Karney 2013, henceforth referred to as AG), I presented algorithms for solving geodesic problems on an ellipsoid of revolution. This built on the classic work of Bessel (1825) and Helmert (1880) to give full double precision accuracy for ellipsoids with flattening satisfying \(|f| \le \frac{1}{100}\). (“Double precision” here refers to floating point numbers with 53 bits—about 16 decimal digits—of precision.) A small improvement, described in Appendix A.2, extended this range \(|f| \le \frac{1}{50}\). Perhaps the most consequential innovation of AG was the reliable solution of the inverse geodesic problem, computing the shortest path between two arbitrary points on the ellipsoid. (Previously, the state of the art was provided by Vincenty (1975), which sometimes fails to converge.) Another important advance was turning the formulation of Danielsen (1989) for the area between a geodesic segment and the equator into an algorithm for computing accurately the area of any geodesic polygon, a polygon whose edges are geodesics. These algorithms allow geodesic problems to be reliably solved and this, in turn, has meant that they have been widely adopted in a variety of disciplines.

The present paper seeks to extend the treatment to arbitrary ellipsoid (in this paper, the term “ellipsoid” should be understood to mean “ellipsoid of revolution”, i.e., a “spheroid”). The condition \(|f| \le \frac{1}{50}\) covers all terrestrial applications. However, for the giant planets in our solar system, f lies outside this range (for Saturn, we have \(f\approx \frac{1}{10}\)). For geodesic calculations, the solutions of the direct and inverse problems, this is a matter of realizing the formulation in terms of elliptic integrals, originally given by Legendre, as a usable algorithm. Generalizing the computation of the area of a geodesic polygon required a novel approach using the discrete sine transform.

Although the mathematical techniques used in this paper might be daunting for some readers, the concepts of geodesics with their interlocking definitions in terms of straightest lines and shortest paths are readily grasped. Furthermore, geodesics have well understood properties, for example:

  • The shortest path from a point to a line intersects the line at right angles;

  • Shortest geodesics satisfy the triangle inequality and the other conditions that define a metric space;

  • Consequently, nearest neighbor problems can be solved efficiently using a vantage-point tree (Uhlmann 1991; Yianilos 1993).

It is therefore possible, indeed it might be preferable, to tackle many problems involving geodesics by treating the implementation described here as a reliable “black box.”

The appendices cover various small improvements in the original algorithms given in AG and compare the algorithms presented here with recent work by Nowak and Nowak Da Costa (2022).

2 Geodesics in terms of elliptic integrals

2.1 Formulation of Legendre and Cayley

One of the earliest systematic studies of geodesics on an ellipsoid was given by Legendre (1811, §§126–129) who cast the solution in terms of his elliptic integrals. However, when Bessel later wished to obtain numerical results for the course of a geodesic, he stated “Because the tools to compute these special functions [elliptic integrals] are not yet sufficiently versatile, we instead develop a series solution which converges rapidly because \(e^2\) is so small.” Bessel’s approach, series expansions valid for small flattening, became the mainstay for solving the geodesic problem. For a long time, elliptic integrals continued to be difficult to compute conveniently; for example, Abramowitz and Stegun (1964, Chap. 19) only provides tables for the elliptic integrals for selected values of the modulus and parameter. This obstacle was only removed with the work of Carlson (1995) as we describe below.

Before we give Legendre’s solution, let us review the formulation in terms of the auxiliary sphere; see Fig. 2 of AG. This is a construct used implicitly by Legendre and explicitly by Bessel that maps the path of the geodesic on the ellipsoid onto a great circle on the sphere. The latitude \(\phi \) is mapped to the parametric latitude \(\beta = \tan ^{-1}\bigl ((1-f)\tan \phi \bigr )\), the azimuth \(\alpha \) is preserved, and \(\lambda \) and s are related to the corresponding spherical variables \(\omega \) and \(\sigma \) by (Karney 2011, Appendix A)

$$\begin{aligned} \frac{1}{a} \frac{\mathrm ds}{\mathrm d\sigma } = \frac{\mathrm d\lambda }{\mathrm d\omega } = \sqrt{1-e^2\cos ^2\beta }. \end{aligned}$$
(1)

Substituting the results from spherical trigonometry,

$$\begin{aligned} \frac{\mathrm d\omega }{\mathrm d\sigma }&= \frac{\sin \alpha }{\cos \beta } = \frac{\sin \alpha _0}{\cos ^2\beta }, \end{aligned}$$
(2)
$$\begin{aligned} \sin \beta&= \cos \alpha _0\sin \sigma , \end{aligned}$$
(3)

we obtain

$$\begin{aligned} \frac{s}{b}&= \int _0^\sigma \sqrt{1 + k^2\sin ^2\sigma '}\,\mathrm d\sigma ', \end{aligned}$$
(4)
$$\begin{aligned} \lambda&= (1-f)\sin \alpha _0 \int _0^\sigma \frac{\sqrt{1 + k^2\sin ^2\sigma '}}{1-\cos ^2\alpha _0\sin ^2\sigma '}\,\mathrm d\sigma ', \end{aligned}$$
(5)
$$\begin{aligned} J(\sigma )&= \int _0^\sigma \frac{k^2\sin ^2\sigma '}{\sqrt{1 + k^2\sin ^2\sigma '}} \,\mathrm d\sigma ', \end{aligned}$$
(6)

where \(k = e' \cos \alpha _0\). \(J(\sigma )\) is a term appearing in the expressions for the reduced length m and the geodesic scale M; both these quantities, which were introduced by Gauss (1828), describe the behavior of nearby geodesics; see Sect. 3 of AG. In brief, m is the ratio of the final linear separation to the initial angular separation for two nearby geodesics that are initially intersecting; M is the ratio of the final to initial linear separations of two nearby geodesics that are initially parallel. These expressions can be put in the form of Legendre’s elliptic integrals,

$$\begin{aligned} \frac{s}{b}&= E(\sigma , ik), \end{aligned}$$
(7)
$$\begin{aligned} \lambda&= (1 - f) \sin \alpha _0 G(\sigma , \cos ^2\alpha _0, ik), \end{aligned}$$
(8)
$$\begin{aligned} J(\sigma )&= E(\sigma , ik) - F(\sigma , ik), \end{aligned}$$
(9)

where

$$\begin{aligned} G(\phi ,\alpha ^2,k)&= \int _0^\phi \frac{\sqrt{1 - k^2\sin ^2\theta }}{1 - \alpha ^2\sin ^2\theta }\,\mathrm d\theta \nonumber \\&=\frac{k^2}{\alpha ^2}F(\phi , k) +\biggl (1-\frac{k^2}{\alpha ^2}\biggr )\Pi (\phi , \alpha ^2, k), \end{aligned}$$
(10)

and \(F(\phi , k)\), \(E(\phi , k)\), and \(\Pi (\phi , \alpha ^2, k)\), are incomplete elliptic integrals (Olver et al. 2010, §19.2(ii)). Equations (7) and (8) are similar to the expressions given by Legendre; the only significant difference is that our formulas involve an imaginary modulus, ik, because we choose our origin for geodesics, \(\sigma = 0\), to be at the node (rather than at the vertex, the point of maximum latitude, which was Legendre’s choice). This is inconsequential because the definitions of the elliptic integrals involve the square of the modulus; furthermore, for prolate ellipsoids, ik is real.

Equations (5) and (8) for \(\lambda \) are inconvenient to use in practice because the integrands are arbitrarily large as the geodesic passes close to a pole leading to a nearly discontinuous change in \(\lambda \). Of course, this is the behavior of great circles as they approach a pole, and it is easily described by spherical trigonometry; see Eq. (12) in AG. It is therefore advantageous to split the expression for \(\lambda \) into two terms, one involving \(\omega \) to describe the basic behavior plus another, proportional to f and involving special functions, to describe the ellipsoidal correction. This can be achieved by integrating

$$\begin{aligned} \frac{\mathrm d(\lambda -\omega )}{\mathrm d\omega } = \sqrt{1-e^2\cos ^2\beta } - 1, \end{aligned}$$
(11)

which gives

$$\begin{aligned} \lambda =\omega - e^2\sin \alpha _0 \int _0^\sigma \frac{1}{1+(1-f)\sqrt{1 + k^2\sin ^2\sigma '}} \,\mathrm d\sigma ', \nonumber \\ \end{aligned}$$
(12)

where the near-singular behavior of \(\lambda \) is captured by \(\omega \) with the difference given by a well behaved integral. This is the starting point for Bessel’s series solution.

Cayley (1870) found a similar way to transform Eq. (8) to give

$$\begin{aligned} \lambda = \chi - \frac{e'^2}{\sqrt{1+e'^2}}\sin \alpha _0 H(\sigma , -e'^2, ik), \end{aligned}$$
(13)

where

$$\begin{aligned} \tan \chi&= \sqrt{\frac{1+e'^2}{1+k^2\sin ^2\sigma }}\tan \omega , \end{aligned}$$
(14)
$$\begin{aligned} H(\phi , \alpha ^2, k)&= \int _0^\phi \frac{\cos ^2\theta }{(1-\alpha ^2\sin ^2\theta )\sqrt{1-k^2\sin ^2\theta }} \,\mathrm d\theta \nonumber \\&= \frac{1}{\alpha ^2} F(\phi , k) + \biggl (1 - \frac{1}{\alpha ^2}\biggr ) \Pi (\phi , \alpha ^2, k). \end{aligned}$$
(15)

Now, the main contribution to \(\lambda \) is \(\chi \), a variable that behaves similarly to the longitude near a pole, and the elliptic integral is relegated to a correction term; not only is it multiplied by \(e'^2\) but the parameter of the elliptic integral, \(-e'^2\) is also small. This allows the longitude to be computed more accurately. The equivalence of the two expressions for \(\lambda \), Eqs. (8) and (13), follows from Olver et al. (2010, Eq. 19.7.8).

2.2 Numerical evaluation of the elliptic integrals

Carlson (1995) provided high-quality algorithms that allowed symmetric elliptic integrals (Olver et al. 2010, §19.16) to be computed to arbitrary accuracy. Legendre’s elliptic integrals can then be simply evaluated using Olver et al. (2010, §19.25(i)). In our implementation, we also use the amendments to Carlson’s paper given in Olver et al. (2010, §19.36(i)).

With these methods in hand, we review the solutions of the direct and inverse geodetic problems. We first discuss the direct problem, determining the final position \(\phi _2\), \(\lambda _2\) and azimuth \(\alpha _2\), given a starting position \(\phi _1\), \(\lambda _1\), initial azimuth \(\alpha _1\), and distance \(s_{12}\). The starting information allows us to calculate the parameter characterizing the geodesic, \(\alpha _0\). Equations (7) and (13) provide a reliable way to follow a geodesic arbitrarily far given its starting position and initial azimuth. The independent parameter in these equations is \(\sigma \), the arc length on the auxiliary sphere. Because we are given the distance \(s_{12}\) to the second point, we need to invert Eq. (7) to find \(\sigma \) in terms of s. This presents no difficulty; s is a monotonically increasing function of \(\sigma \) so there’s a unique solution and, provided we are sufficiently close to the solution, Newton’s method can be used to obtain an accurate result; the derivative needed for Newton’s method, \(\mathrm ds/\mathrm d\sigma \) is given by Eq. (1). This allows the direct geodesic problem to be solved in the same way as in AG.

Fig. 1
figure 1

A simple closed geodesic for an ellipsoid with \(b/a=\frac{1}{4}\). The geodesic starts on the equator with azimuth \(\alpha _0 \approx 51.24052^\circ \) and completes two full oscillations about the equator before returning to the starting longitude: a a top view; b a view at an inclination of \(25^\circ \) to the equatorial plane; “N” marks the north pole. The solid and dotted lines show the visible and hidden portions of the geodesic

Full size image

Figure 1 illustrates a simple (i.e., non-self-intersecting) closed geodesic on an ellipsoid with \(f = \frac{3}{4}\) or \(b/a=\frac{1}{4}\). This is one of a class of “non-standard” closed geodesics that are neither equatorial nor meridional and which appear for sufficiently oblate ellipsoids, \(a > 2b\) (Klingenberg 1982, §3.5.19).

2.3 Solving the inverse problem

The inverse problem is more challenging. Here, we are given the endpoints \(\phi _1\), \(\lambda _1\) and \(\phi _2\), \(\lambda _2\), and wish to determine the shortest distance \(s_{12}\) and azimuths \(\alpha _1\) and \(\alpha _2\). Because we have no direct way to determine the azimuth at the node, \(\alpha _0\), we must resort to an iterative solution. The traditional method for solving this problem, given by Helmert (1880) and used by Vincenty (1975), is functional iteration which fails to converge in cases where the endpoints are nearly antipodal.

A key contribution of AG was to reduce the determination of \(\alpha _0\) to a simple one-dimensional root-finding problem. I briefly recapitulate the method here. Let us assume we have already dealt with equatorial and meridional geodesics. Equatorial geodesics are those with \(\phi _1 = \phi _2 = 0\) and, for oblate ellipsoids, the additional condition that \(|\lambda _{12}| \le (1-f)\pi \); for these we have \(\alpha _1 = \pm \frac{1}{2}\pi \). Meridional geodesics are those with \(\lambda _{12} = 0\) or \(\pm \pi \) and, for prolate ellipsoids, the additional condition that \(m_{12} \ge 0\); and for these geodesics, we have \(\alpha _1 = 0\) or \(\pm \pi \).

For the remaining, oblique, cases, we can, without loss of generality, specialize to the case

$$\begin{aligned} \phi _1 \le 0,\quad \phi _1 \le \phi _2 \le -\phi _1,\quad 0 < \lambda _{12} \le \pi . \end{aligned}$$
(16)

Note that \(\lambda _{12} = 0\) is excluded—this case corresponds to a meridional geodesic; the case \(\lambda _{12} = \pi \) only arises for \(f < 0\) (for \(f \ge 0\) this also corresponds to a meridional geodesic).

For a particular trial azimuth \(\alpha ^*_1\), we solve the hybrid problem: given \(\phi _1\), \(\alpha ^*_1\), and \(\phi _2\), compute the longitude difference \(\lambda _{12}(\alpha ^*_1)\) corresponding to the first intersection of the geodesic with the circle of latitude \(\phi _2\). The solution of the inverse problem is then given by finding \(\alpha ^*_1 = \alpha _1\) which solves

$$\begin{aligned} \lambda _{12}(\alpha _1) = \lambda _{12}, \end{aligned}$$
(17)

with \(0< \alpha _1 < \pi \). We use \(\alpha _1\) as the control variable for this root-finding problem; this is a more convenient choice than \(\alpha _0\) and, once this is found, we can readily obtain \(\alpha _0\). In AG, we solved this equation using Newton’s method. This converges very rapidly but requires a good starting guess for the solution and this, in turn, required a careful analysis of the case of nearly antipodal endpoints.

Fig. 2
figure 2

Solution to the hybrid problem for a \(\phi _1 = \phi _2 = 0\) and \(b/a = {\textstyle \frac{1}{2}}\) (marked with crosses); b \(\phi _1 = -30^\circ \), \(\phi _2 = 20^\circ \), and \(b/a = 2\) (marked with circles). This illustrates the functional dependence of \(\lambda _{12}(\alpha ^*_1)\)

Full size image

This technique worked well with ellipsoids with small flattening. However, in the general case, the starting guess was sometimes sufficiently far from the true solution that Newton’s method failed. Before describing how we might deal with this, let us first review the properties of \(\lambda _{12}(\alpha ^*_1)\). First of all, we have \(\lambda _{12}(0) = 0\) and \(\lambda _{12}(\pi ) = \pi \). Furthermore, with one exception, \(\lambda _{12}(\alpha ^*_1)\) varies continuously so Eq. 17 is guaranteed to have a solution. The exception is illustrated in Fig. 2(a); this is the case with endpoints on the equator for an oblate ellipsoid and \(\lambda _{12}(\alpha ^*_1)\) is discontinuous at \(\alpha ^*_1 = {\textstyle \frac{1}{2}}\pi \). Because equatorial and meridional geodesics have already been handled, this only occurs for \({(1-f)\pi }< \lambda _{12} < \pi \) and we can restrict the possible range of azimuth to \({\textstyle \frac{1}{2}}\pi< \alpha ^*_1 < \pi \). Within this range, \(\lambda _{12}(\alpha ^*_1)\) is continuous and a solution exists. Finally, there is exactly one root and hence the derivative \(\mathrm d\lambda _{12}(\alpha ^*_1)/\mathrm d\alpha ^*_1 = \lambda _{12}'(\alpha ^*_1)\) is non-negative at that root. In this context, it is important to consider the case of a prolate ellipsoid; see Fig. 2(b). Even though \(\lambda _{12}(\alpha ^*_1)\) has an interior maximum, \(\lambda _{12}(\alpha ^*_1) = \lambda _{12}\) still has only a single root provided that \(\lambda _{12} \le \pi \) (stipulated in the initial setup of the inverse problem) and that \(\alpha _1 < \pi \) (because we have already handled meridional geodesics).

With this background, even the crudest root-finding techniques, e.g., bisection or regula falsi, can be used to find the solution. We shall continue to use Newton’s method on account of its rapid convergence. However, we back it up with the bisection method to deal with two possible problems with Newton’s method: (a) \(\lambda _{12}'(\alpha ^*_1)\) is so small that the new value of \(\alpha ^*_1\) is outside the range \((0,\pi )\) or (b) \(\lambda _{12}'(\alpha ^*_1)\) is negative, in which case the new value of \(\alpha ^*_1\) is further from the true root. During the course of the Newton iteration, we maintain a bracketing range for \(\alpha _1 \in (\alpha _{1-}, \alpha _{1+})\). Initially, we have \(\alpha _{1-} = 0\), \(\alpha _{1+} = \pi \). On each iteration, we update one of \(\alpha _{1\pm }\) to \(\alpha ^*_1\) depending on the sign of \(\lambda _{12}(\alpha ^*_1) - \lambda _{12}\) provided that this results in a tighter bracket for \(\alpha _1\). Whenever either of the problems, (a) or (b) listed above, occurs, we restart Newton’s method with a new starting guess \(\alpha ^*_1 = {\textstyle \frac{1}{2}}(\alpha _{1-} + \alpha _{1+})\); and we will be able to bisect the bracket after the next iteration.

We can use the spherical solution to get a starting value for \(\alpha ^*_1\) with a value of 0 or \(\pi \) (i.e., outside the allowed range for \(\alpha _1\)) replaced by \(\frac{1}{2}\pi \). (See also Appendix A.3 for additional caveats for the case of equatorial endpoints.) Typically, the bisection safety mechanism is invoked at most once when solving the inverse problem. The net result is that this provides a reliable (always converging) and fast (rapidly converging) method of solving the inverse geodesic problem for arbitrary ellipsoids. In some cases, the inverse problem does not have a unique solution—there are multiple shortest geodesics; these cases are readily handled as discussed in Appendix A.1.

3 The area integral

3.1 Formulation of Danielsen

In AG, \(S_{12}\) was defined as the area between a geodesic segment \(s_{12}\) and the equator, i.e., the area of the geodesic quadrilateral \(AF\!H\!B\) in Fig. 1 of AG. Following Danielsen (1989), this was cast as an integral as follows:

$$\begin{aligned} S_{12} = S(\sigma _2) - S(\sigma _1),\quad S(\sigma ) = c^2\bigl (\alpha + p(\sigma )\bigr ), \end{aligned}$$
(18)

where

$$\begin{aligned} c^2 = \frac{a^2}{2} + \frac{b^2}{2} \frac{\tanh ^{-1}e}{e} \end{aligned}$$
(19)

is the authalic radius squared, and

$$\begin{aligned} p(\sigma )&= \int _{\pi /2}^\sigma q(\sigma ')\,\mathrm d\sigma ', \end{aligned}$$
(20)
$$\begin{aligned} q(\sigma )&= - A_4 \Delta t(e'^2, k^2\sin ^2\sigma ) \frac{\sin \sigma }{2}, \end{aligned}$$
(21)
$$\begin{aligned} A_4&= \frac{e^2a^2}{c^2}\cos \alpha _0 \sin \alpha _0, \end{aligned}$$
(22)
$$\begin{aligned} \Delta t(x,y)&= \frac{t(x)-t(y)}{x-y}, \end{aligned}$$
(23)
$$\begin{aligned} t(x)&= x + \sqrt{1+x}\, \frac{\sinh ^{-1}\sqrt{x}}{\sqrt{x}}. \end{aligned}$$
(24)

This is a change from the notation in Eq. (59) of AG; the relationship is

$$\begin{aligned} p(\sigma ) = A_4 I_4(\sigma ); \end{aligned}$$
(25)

this simplifies the assessment of the error in \(S_{12}\). The expression \(\Delta t(x,y)\) is the divided difference of t(x). In the limit \(y \rightarrow x\), we have \(\Delta t(x,x)=t'(x)\). If x and y are distinct and have the same sign, there is a systematic way to express \(\Delta t(x,y)\) in a form that avoids excessive roundoff errors (Kahan and Fateman 1999).

The area of an arbitrary geodesic polygon is easily computed by summing the signed area contributions \(S_{12}\) for each edge. An adjustment to this total needs to be made if the polygon encircles a pole as described in AG.

The challenge is that the integral, Eq. (20), cannot be carried out in closed form. However, both \(q(\sigma )\) and \(p(\sigma )\) are periodic functions and so may be written as Fourier series:

$$\begin{aligned} q(\sigma )&= \sum _{l=0}^\infty Q_l \sin \bigl (2(l+{\textstyle \frac{1}{2}})\sigma \bigr ), \end{aligned}$$
(26)
$$\begin{aligned} Q_l&= \frac{4}{\pi }\int _0^{\pi /2} q(\sigma ) \sin \bigl (2(l+{\textstyle \frac{1}{2}})\sigma \bigr )\,\mathrm d\sigma , \end{aligned}$$
(27)
$$\begin{aligned} p(\sigma )&= \sum _{l=0}^\infty P_l \cos \bigl (2(l+{\textstyle \frac{1}{2}})\sigma \bigr ), \end{aligned}$$
(28)
$$\begin{aligned} P_l&= -\frac{Q_l}{2(l+\frac{1}{2})}. \end{aligned}$$
(29)

Here, Eqs. (26) and (27) constitute a Fourier transform pair. In the next sections, we explore two ways to compute \(P_l\) approximately.

3.2 Using a Taylor series

In AG, we gave Taylor series expansions for \(P_l\) using \(e'^2\) and \(k^2\) as small parameters. For reasons explained in Appendix A.5, this is a poor choice of expansion parameters for eccentric ellipsoids; so we replace the expansion parameters by n and \(\epsilon \), Eqs. (43), resulting in the expansion Eqs. (44) which should be used instead of AG’s Eqs. (63). This expansion may be used to approximate \(p(\sigma )\) with

$$\begin{aligned} {\tilde{p}}^{(N)}(\sigma )&= \sum _{l=0}^\infty {\tilde{P}}_l^{(N)} \cos \bigl (2(l+{\textstyle \frac{1}{2}})\sigma \bigr ), \end{aligned}$$
(30)
$$\begin{aligned} {\tilde{P}}_l^{(N)}&= A_4 C_{4l}, \end{aligned}$$
(31)

where \(A_4\) is defined in Eq. (22) and the superscript (N) gives the order at which the Taylor series is truncated; thus Eqs. (44) with the ellipses dropped corresponds to \(N=6\). \(A_4\) is proportional to \(e^2\), so the highest-order terms in Eqs. (44) are \(O(f^6)\); furthermore, we have \(\tilde{P}_l^{(N)} = 0\) for \(l\ge N\).

One possibility for computing the area on eccentric ellipsoids is extending the Taylor series expansion to higher order. Thus, when I first implemented the solution of the geodesic problem in terms of elliptic integrals as described in Sect. 2, I computed the area using \({\tilde{p}}^{(30)}(\sigma )\); this gives full double precision accuracy for \({\textstyle \frac{1}{2}}\le b/a \le 2\). However, this approach quickly becomes cumbersome for extremely eccentric ellipsoids:

  • The value of N required to give full accuracy becomes very large;

  • The cost of determining the coefficients in the Taylor series becomes impractical—it takes Maxima (2022) 17 min for \(N=32\) (tolerable for generating code) but 19 h to extend this to \(N=64\) (which is starting to get prohibitive);

  • The storage required for the coefficients becomes excessive, scaling as \(N^3/6\);

  • Once the ellipsoid is selected, evaluating the coefficients of \(\epsilon ^j\) in Eqs. (44) requires \(O(N^3)\) operations;

  • Once \(\alpha _0\) is given, evaluating \({\tilde{P}}_l^{(N)}\) requires \(O(N^2)\) operations.

3.3 Quadrature using discrete sine transforms

To overcome the shortcomings of Taylor series for large eccentricities, I explored evaluating Eq. (20) using the quadrature method of Clenshaw and Curtis (1960). Given a definite integral,

$$\begin{aligned} \int _{-1}^1 g(x)\, \mathrm dx, \end{aligned}$$

this method entails expanding g(x) as Chebyshev polynomials, or, equivalently, \(g(\cos \theta )\) as a Fourier series. This allows the definite integral to be computed to high accuracy and, as a bonus, the indefinite integral can be computed. The accuracy is a result of the fast convergence of the trapezoidal rule for integrals of a periodic function over a full period; in this case, the integrals give the Fourier coefficients.

However, this is a roundabout procedure given that we are starting with a periodic function. Instead, we directly evaluate Eq. (27) using the trapezoidal rule and evaluate \(p(\sigma )\) by using the resulting Fourier coefficients in Eqs. (28) and (29). Thus, we approximate \(Q_l\) with

$$\begin{aligned} {\hat{Q}}_l^{(N)} = \frac{4h_N}{\pi }\mathop {\mathrm {\sum \nolimits ^\prime }}\limits _{j=0}^{N-1} q\left( (j+1)h_N\right) \sin \bigl (2(l+{\textstyle \frac{1}{2}})(j+1)h_N\bigr ),\nonumber \\ \end{aligned}$$
(32)

where \(h_N = \pi /(2N)\), we take \(0 \le l < N\), and the prime on the summation sign indicates that the last term in the sum should be halved. Because \(q(\sigma )\) is an analytic periodic function, \({\hat{Q}}_l^{(N)}\) converges to the true value \(Q_l\) faster than any power of N, as \(N\rightarrow \infty \). Equation (32) is, in fact, the discrete sine transform, DST, of type III (Wang and Hunt 1985). Its inverse is a type II DST,

$$\begin{aligned} q((j+1)h_N) = \sum _{l=0}^{N-1} {\hat{Q}}_l^{(N)} \sin \bigl (2(l+{\textstyle \frac{1}{2}})(j+1)h_N\bigr ); \end{aligned}$$
(33)

this holds for all integer j. Rewriting this as a continuous function,

$$\begin{aligned} {\hat{q}}^{(N)}(\sigma ) = \sum _{l=0}^{N-1} {\hat{Q}}_l^{(N)} \sin \bigl (2(l+{\textstyle \frac{1}{2}})\sigma \bigr ), \end{aligned}$$
(34)

we see that \({\hat{q}}^{(N)}(\sigma )\) is likely to be an excellent approximation to \(q(\sigma )\) given that the functions coincide wherever \(\sigma \) is a multiple of \(h_N\). In the same way, we approximate \(p(\sigma )\) by

$$\begin{aligned} {\hat{p}}^{(N)}(\sigma )&= \sum _{l=0}^\infty {\hat{P}}_l^{(N)} \cos \bigl (2(l+{\textstyle \frac{1}{2}})\sigma \bigr ), \end{aligned}$$
(35)
$$\begin{aligned} {\hat{P}}_l^{(N)}&= {\left\{ \begin{array}{ll} \displaystyle -\frac{{\hat{Q}}_l^{(N)}}{2(l+\frac{1}{2})},&{}\text {for }l < N,\\ 0, &{}\text {for }l \ge N. \end{array}\right. } \end{aligned}$$
(36)

A convenient way to compute DSTs is using the fast Fourier transform, FFT (Cooley and Tukey 1965). This reduces the cost of evaluating Eq. (32) for \(0\le l <N\) from \(O(N^2)\) to \(O(N\log N)\). Furthermore, Gentleman (1972) points out that the FFT results in smaller roundoff errors compared with other methods. Equation (35) can be evaluated using Clenshaw (1955) summation; this method is fast (the cosine terms are computed by a recurrence relation) and because \({\hat{Q}}_l^{(N)}\) decays rapidly with l and because the summation starts at high values of l, it is also accurate.

As we shall see, the error \(\max (| {\hat{q}}^{(N)}(\sigma ) -q(\sigma ) |)\) decays exponentially with N; thus doubling N roughly squares the error. Furthermore, it is possible to double N rather inexpensively. First, generate the type IV DST for \(q(\sigma )\),

$$\begin{aligned} {\check{Q}}_l^{(N)} = \frac{4h_N}{\pi }\sum _{j=0}^{N-1} q((j+{\textstyle \frac{1}{2}})h_N) \sin \bigl (2(l+{\textstyle \frac{1}{2}})(j+{\textstyle \frac{1}{2}})h_N\bigr ).\nonumber \\ \end{aligned}$$
(37)

The sums for computing \({\hat{Q}}_l^{(N)}\), Eq. (32), and \({\check{Q}}_l^{(N)}\) involve evaluating \(q(\sigma )\) at intervals of \({\textstyle \frac{1}{2}}h_N = h_{2N}\); it follows that \({\hat{Q}}_l^{(2N)}\) is just the mean of \({\hat{Q}}_l^{(N)}\) and \({\check{Q}}_l^{(N)}\). These coefficients for \(N\le l < 2N\) are

$$\begin{aligned} {\hat{Q}}_l^{(N)} = - {\hat{Q}}_{2N-1-l}^{(N)},\quad {\check{Q}}_l^{(N)} = + {\check{Q}}_{2N-1-l}^{(N)}. \end{aligned}$$
(38)

Thus, \({\hat{Q}}_l^{(2N)}\) is given by

$$\begin{aligned} {\hat{Q}}_{N-1/2\pm g}^{(2N)} = \frac{{\check{Q}}_{N-1/2-|g|}^{(N)} \mp {\hat{Q}}_{N-1/2-|g|}^{(N)}}{2}, \end{aligned}$$
(39)

for \(g = \frac{1}{2}, \frac{3}{2}, \ldots , N - \frac{1}{2}\). This method avoids wasting the results \({\hat{Q}}_l^{(N)}\) which have already been computed (Gentleman 1972).

3.4 Convergence of approximate Fourier series

Fig. 3
figure 3

Dependence of the errors in the individual Fourier coefficients, \(P_l\), on l for \(N=30\). The crosses, resp. squares, show the errors using the DST method, \(\delta {\hat{P}}_l^{(N)}\), resp. the Taylor series method, \(\delta \tilde{P}_l^{(N)}\). The case illustrated is for a geodesic on an ellipsoid with \(b/a=\frac{1}{4}\), \(\alpha _0=\frac{1}{4}\pi \). The line shows the exact value of \(|P_l|\)

Full size image

In assessing how well \({\hat{p}}^{(N)}(\sigma )\), Eq. (35), approximates \(p(\sigma )\), we start by looking at the error in the individual Fourier coefficients, namely \(\delta {\hat{P}}_l^{(N)} = | {\hat{P}}_l^{(N)} - P_l |\) for a given N. For illustrative purposes, we consider \(b/a = \frac{1}{4}\), \(\alpha _0 = \frac{1}{4} \pi \), and consider \(N = 30\); we compute the “exact” \(P_l\) using a much larger value N. The resulting errors \(\delta {\hat{P}}_l^{(30)}\) are shown in Fig. 3 as crosses. The figure also shows \(|P_l|\) (the solid line) decaying exponentially as a function of l. Obviously \(\delta {\hat{P}}_l^{(30)}\), which equals \(|P_l|\) for \(l\ge 30\), shows the same exponential decay. However, we also have exponential decay in \(\delta {\hat{P}}_l^{(30)}\) as l decreases from 29 to 0; indeed, \(\delta {\hat{P}}_l^{(30)}\) is roughly symmetric about \(l = 30\). This behavior is predicted by Eq. (39), since we can approximate \(\delta {\hat{P}}_l^{(N)}\) by \(|{\hat{P}}_l^{(N)} - {\hat{P}}_l^{(2N)}|\), with

$$\begin{aligned} {\hat{P}}_{N-1/2\pm g}^{(N)} - {\hat{P}}_{N-1/2\pm g}^{(2N)} = \frac{{\hat{Q}}_{N-1/2-|g|}^{(N)} - {\check{Q}}_{N-1/2-|g|}^{(N)}}{4(N\pm g)}; \end{aligned}$$
(40)

the variation of the numerator in this expression (symmetric in g) is much stronger than the variation of the denominator. We conclude that the error arising from using approximate Fourier coefficients, \({\hat{P}}_l^{(N)}\) instead of \(P_l\), is approximately the same as the error from truncating the Fourier sum in Eq. (35).

Figure 3 also shows the errors in the Taylor series coefficients \(\delta {\check{P}}_l^{(N)} = |{\check{P}}_l^{(N)} - P_l|\) for \(N=30\) as squares. The error is, of course, the same as for \(\delta {\hat{P}}_l^{(30)}\) for \(l \ge 30\); but for \(l < 30\), the error is several orders of magnitude bigger.

Fig. 4
figure 4

Dependence of the overall error in the Fourier series approximation to \(p(\sigma )\) on N. The crosses, resp. squares, show the errors for the DST method, \(\Delta {\hat{P}}^{(N)}\), resp. the Taylor series method, \(\Delta {\tilde{P}}^{(N)}\). This is for the same case as Fig. 3, namely \(b/a=\frac{1}{4}\), \(\alpha _0=\frac{1}{4}\pi \). The horizontal line indicates the roundoff limit for double precision \(2^{-53} \approx 10^{-16}\)

Full size image

Figure 3 gives the errors in the individual Fourier coefficients for a particular N. To gauge how the overall error in \({\hat{p}}^{(N)}(\sigma )\), Eq. (35), varies with N, we combine the errors in the coefficients, as follows:

$$\begin{aligned} \Delta {\hat{P}}^{(N)} = \sum _{l=0}^\infty \delta {\hat{P}}_l^{(N)}. \end{aligned}$$
(41)

We similarly define \(\Delta {\tilde{P}}^{(N)}\) for the coefficients obtained by Taylor expansion. Figure 4 shows how these error metrics vary as N is increased from 10 to 60 for \(b/a = \frac{1}{4}\) and \(\alpha _0 = \frac{1}{4}\pi \) (the same parameters used for Fig. 3). Both show exponential decay with \(\Delta {\hat{P}}^{(N)}\) being several orders of magnitude smaller and decaying at a slightly faster rate.

3.5 Determining the number of points in the quadrature

In practical applications, we will evaluate the area using Eqs. (18) and have to contend both with truncation errors, because we use a finite value of N in computing \(p(\sigma )\), and with roundoff errors, because of the limited precision of the floating point number system. We would like to pick N large enough that the truncation errors are less than the roundoff errors, and for this, we require that \(\Delta {\hat{P}}^{(N)} \le 2^{-m}\) where m is the number of bits in the fraction of a floating point number. For double precision, we typically have \(m = 53\) and, substituting the value of \(c\approx 6371\,\textrm{km}^2\) for the earth, this gives a truncation error in \(S(\sigma )\) of \(2^{-53}c^2 \approx 50\,\textrm{cm}^2\). The horizontal line in Fig. 4 labels this error threshold, \(2^{-53}\). This shows that for this case, we need \(N\ge 34\) for the DST method compared to \(N\ge 56\) when using a Taylor expansion.

Fig. 5
figure 5

Heavy line shows the minimum number of points N required in the DST to ensure double precision accuracy as a function of n for \(|n| \le 0.99\) and for arbitrary \(\alpha _0\). If we limit the possible values of N to \(2\times 2^j\) and \(3\times 2^j\) with integer \(j \ge 1\), then the result is the staircase shown as a light line

Full size image

In principle, we could reestimate N for each geodesic, i.e., for every n and \(\alpha _0\). (Since we are dealing with both oblate and prolate ellipsoids of arbitrary eccentricity, we use n as the measure of eccentricity instead of \(e = 2\sqrt{n}/(1+n)\) or \(f = 2n/(1+n)\).) However, that would require extra machinery to do the estimation at run time. Instead, we estimate N for \(-0.99 \le n \le 0.99\). For a particular n, we find the minimum N that satisfies the accuracy requirement for all \(\alpha _0\). The result is shown in Fig. 5 as the heavy line. The light line shows the result of restricting the allowed values of N to \(2\times 2^j\) and \(3\times 2^j\) with integer \(j \ge 1\); now, N is the product of small factors (2 and 3), a requirement for the efficient implementation of the FFT. The light line is converted to a simple table allowing N to be looked up at run time based on the value of n.

The conclusion is that the discrete sine transform method provides a convenient and accurate way to evaluate the area integral. In comparison with the Taylor series method: the same accuracy is obtained with substantially fewer terms (\(N = 34\) for the DST versus \(N = 56\) for the Taylor series in the example illustrated in Fig. 4). Just as important, N can be easily adjusted for a particular n, Fig. 5 (no need to compute a new Taylor series), the memory required is O(N) instead of \(O(N^3)\), and the setup time for a particular geodesic is \(O(N\log N)\) instead of \(O(N^3)\).

Table 1 Results for a geodesic segment with \(\phi _1 = 0\), \(\lambda _1 = 0\), \(\alpha _1 = \frac{1}{4}\pi \), \(\sigma _{12} = {\textstyle \frac{1}{2}}\pi \) on an oblate ellipsoid with \(a = 6\,400\,\textrm{km}\) and with various values of the third flattening, n
Full size table
Table 2 Analogous results to Table 1, but for prolate ellipsoids
Full size table

4 Implementation

The algorithms detailed in Sects. 2 and 3 have been implemented in version 2.1.2 of the open-source C++ library, GeographicLib (Karney 2022). The solution of the geodesic problems in terms of elliptic integrals was added in version 1.25 (2012). At that time, the area was computed with a 30th-order Taylor series as described in Sect. 3.2. This was replaced by the implementation using DSTs, Sect. 3.3, in version 2.1 (2022). These algorithms supplemented the earlier ones described in AG since, as we see below, these remain the preferred methods for terrestrial ellipsoids.

We use the KISS FFT package (Boergerding 2021) for performing the DST. Even though this package has no direct support for the DSTs, as a “header-only” package, it is easier to integrate into GeographicLib than the more capable FFTW package (Frigo and Johnson 2021).

The main geodesic calculations are also exposed with two utility programs provided with GeographicLib, GeodSolve for geodesic calculations, and Planimeter for measuring the perimeter and area of geodesic polygons. With both utilities, the -E flag uses the algorithms described here for arbitrary ellipsoids.

In developing numerical algorithms, especially those that offer high accuracy, it is very useful to be able to distinguish truncation errors from roundoff errors. For this reason, high precision calculations were added to GeographicLib in version 1.37 (2014). This allowed the library to be compiled with various alternatives to standard double precision (53 bits of precision or 16 decimal digits) for floating point numbers:

  • Extended precision (64 bits or 19 decimal digits) available as the “long double” type on some platforms;

  • Quadruple precision (113 bits or 34 decimal digits) using the boost multi-precision library (Maddock and Kormanyos 2022);

  • Arbitrary precision (selected at the start of execution) using the MPREAL C++ interface (Holoborodko 2022) to the MPFR library (Fousse et al. 2007).

The results in Figs. 3, 4, and 5 and in Tables 1 and 2 were computed with the last option setting the precision to 256 bits (77 decimal digits).

5 Results

A starting point for testing the new algorithms is the test data set given in Karney (2010). This is specific to the WGS84 ellipsoid and as such allows us to compare the new algorithms with the series solutions described in AG. The maximum errors in geodesic calculations when converted to a distance are about \(0.03\,\mu \mathrm m\) which is 2–3 times larger than the corresponding figures for the series solution. The error in the area is about \(0.1\,\mathrm m^2\) in both cases. On this data set, the new routines run about 2.5 times slower. The lower accuracy and speed are to be expected given the additional complexity of the new algorithms.

The lesson is clear: For terrestrial ellipsoids, the algorithms described in AG are to be preferred—they are faster and more accurate than the new routines. In fact, the new routines are still very accurate and fast. Once the eccentricity of the ellipsoid is such that \(|f| > \frac{1}{50}\), the accuracy of the series solution degrades, and the new routines should be used. Let us, therefore, turn to assessing their performance on eccentric ellipsoids.

We start by considering the geodesic that starts on the equator, i.e., at its node \(\phi _1 = 0\), with azimuth \(\alpha _0 = \alpha _1 = \frac{1}{4}\pi \) and compute the longitude at the vertex (the point of maximum latitude), the distance to the vertex, and the area under this geodesic segment. We shall carry out this calculation with various values of the third flattening n. The arc length on the auxiliary sphere, in this case, is \(\sigma _{12} = {\textstyle \frac{1}{2}}\pi \), and, because \(\tan \alpha _0 = 1\), we have \(\tan \phi _2 = (1+n)/(1-n)\).

The results are shown in Tables 1 and 2 for oblate and prolate ellipsoids, respectively. For each n, the first row shows the results with high precision arithmetic rounded to 17 decimal digits. These allow independent verification of the algorithms. The second row shows the results of carrying out the calculation in double precision. In all cases, the roundoff errors are very small—the changes take place in no more than the last 5 bits of the fraction of the binary representation; with \(|n| \le 0.9\), the errors amount to at most 7 units in the last place in the binary representation.

A more challenging test is presented with the data for the administrative boundaries of Poland. (This was kindly provided to me in 2013 by Paweł Pȩdzich of the Warsaw University of Technology.) This set includes a polygonal representation of the boundary of Poland itself. This consists of \(67\,801\) edges with mean length \(53\,\mathrm m\), minimum length \(62\,\textrm{mm}\), and maximum length \(9.7\,\textrm{km}\); the perimeter of the polygon is about \(3\,600\,\textrm{km}\) and its area is about \(313\,000\,\textrm{km}^2\). This is a good test for geodesic algorithms as it represents a realistic “large” problem allowing the overall accuracy and speed of the algorithms to be assessed.

Fig. 6
figure 6

Absolute errors in the calculation of the area \(\Delta A\) (plus signs, scale on the left) and perimeter \(\Delta P\) (diamonds, scale on the right) of a polygonal model of Poland. Here the model is mapped to ellipsoids with various eccentricities (denoted by the third flattening n) by preserving the value of the authalic latitude

Full size image

In performing tests on ellipsoids with varying eccentricities, we need to stipulate how the points are mapped onto the ellipsoid. If we merely keep the latitude fixed, then for extreme oblate, resp. prolate, ellipsoids, all the points end up concentrated near the equator, resp. the poles, making it difficult to interpret the results. We instead map the points keeping the authalic latitude fixed. If, in addition, we choose a so that the area of the ellipsoid matches that of the Earth, the area of the polygon is nearly constant as n is varied, and we assess the errors by computing the difference between the computed area and the true area. For this exercise, we find the true area using quadruple precision arithmetic.

The errors in the area and the perimeter are shown in Fig. 6. Generally speaking, the errors in the area are between \(1\,\mathrm m^2\) and \(10\,\mathrm m^2\) and the errors in the perimeter are between \(1\,\mu \mathrm m\) and \(3\,\mu \mathrm m\).

The elapsed time for these calculations is approximated by

$$\begin{aligned} \bigl (4.5 + N(40+\log _2N)/450\bigr )\,\mathrm {\mu s} \text { per edge}, \end{aligned}$$

where the first term represents the time to compute the perimeter of the polygon, which involves solving the inverse geodesic problem for an edge, and the second term is the additional time required to compute its area. Here, N is the number of terms used in the DST which is given in Fig. 5. The time to compute the area begins to dominate for \(|n| > 0.2\). The area calculation could be made faster if the DSTs were computed using the FFTW package instead of KISS FFT.

6 Discussion

In this paper, we have shown how to solve geodesic problems on an ellipsoid of arbitrary eccentricity and an implementation of these methods is available in GeographicLib. The work naturally divides into two parts:

The method of solving the standard direct and inverse problems was an exercise in applying Carlson’s algorithms for the evaluation of elliptic integrals to the work of Legendre and Cayley. These then replaced the Taylor series expansions for the elliptic integrals used in AG. In addition, the method for solving the inverse problem needed to be adjusted to ensure convergence in all cases.

The evaluation of the area integral, Eq.(20), is simplicity itself: write the periodic integrand as a Fourier series, evaluate the Fourier coefficients by trapezoidal quadrature, and then the integral of the resulting series is trivial. This is an accurate and fast method because (1) the Fourier series of an analytic function converges rapidly, (2) quadrature is very accurate when applied to the full period of a periodic function, (3) the FFT can be used to evaluate all the coefficients efficiently. This method will be useful in other applications, e.g., to evaluate the Abelian integrals appearing in the solution given by Jacobi (1839) for geodesics on a triaxial ellipsoid.

The algorithms have been tested for a wide range of eccentricities, \(|n| \le 0.99\) or \(1/199 \le b/a \le 199\); this encompasses ellipsoids flatter than a pancake and as thin as spaghetti. In principle, the methods would work at more extreme eccentricities; however, it is likely a more specialized approach would be better in these cases. Despite this, the Taylor series method described in AG is preferred for small flattening, \(|f| \le \frac{1}{50}\). This is somewhat more accurate and somewhat faster than the general method given here.

The obvious domain of application for these algorithms is studying planets and other celestial bodies with high eccentricity. But, they also play a role in terrestrial applications. First of all, they allow the errors in the Taylor series method to be measured accurately. Secondly, because the flattening of terrestrial ellipsoids is so small, a starting approximation to the solution is given by spherical trigonometry and this solution can be subsequently refined using ellipsoidal geodesics. Unfortunately, any specifically ellipsoidal behavior might be overlooked, precisely because the flattening is small. Therefore, it behooves the developer to test the algorithms on more eccentric ellipsoids, say, with \(f = \frac{1}{5}\), because then any ellipsoid effects will be much more apparent.

I should point out an alternative method of treating arbitrary eccentricity, namely to integrate the ordinary differential equations for the geodesic numerically. This is the approach taken by Panou and Korakitis (2017), and it has the advantage that it can be generalized to more complex bodies, e.g., a triaxial ellipsoid (Panou and Korakitis 2019). However, for the problem at hand, our approach is preferable: it is fast, it provides high accuracy, the computation time is independent of the length of the geodesic, and it gives the solution of the inverse problem. For bodies that are well-approximated by a triaxial ellipsoid, an attractive approach would be to recast the solution of Jacobi (1839) into a workable numerical algorithm; however, this (non-trivial!) task is beyond the scope of the present paper.