Theory, strict formula derivation and algorithm development for the computation of a geodesic polygon area

Abstract

The paper presents a coherent, original theory and an effective algorithm to calculate the surface area of a geodesic polygon, i.e. the area on an ellipsoid of revolution limited by geodesics, connecting pairs of subsequent vertices of known geographical coordinates. Using the cylindrical coordinate system and introducing an equatorial geodesic triangle, the novel strict recursive formulas were derived for the longitude on the ellipsoid and the distance along the geodesic and the area of a geodesic polygon. The area is determined by constant coefficients, not the values of finite series. The orthodromic properties of geodesics, which led to the detection of bifurcation lines composed of the double solutions of the inverse geodetic problem were analysed. The novel algorithm presented, integrating the problem with geodesic polygon area calculation, is based on the newly proposed parameterisation with two alternative working unknowns and a new formula to obtain the initial value of the geodesic C parameter. It allows solutions for all cases, including near-antipodal or near-vertices, and offers two solutions if bifurcation occurs. The algorithm presented was tested and verified using numerical examples from known publications. The numerical tests proved the accuracy of 1m² of the computed geodesic area using standard IEEE 754 Double Floating Point PC arithmetic, for geodesic lengths up to ten thousand km.

1 Introduction

The computation of the surface area of polygons, whose sides are arcs of geodesics, is related to the knowledge of the exact parameters of geodesics equations. A geodesic is a unique, locally minimising length, curve joining two points on any surface, e.g. a straight line is a geodesic on a plane. Geodesics are used in classical geodetic problems of computing the geodetic coordinates, azimuths and distances on the surface of an ellipsoid of revolution. The formula of the accuracy of ${e}^{\prime 8}$ for calculating the area of a figure between the geodesic and the equator bounded by two meridians (referred to as the area under the geodesic) was formulated by Danielsen in 1989, but the practical application of this formula is difficult in terms of calculation. As a consequence, although storing and calculating data by their original geographical coordinates – currently being achieved easily and more and more accurately, e.g. from Global Navigation Satellite Systems – seems the ideal solution (as noticed also by many others researchers i.e. Canters and De Genst 1997, Leung et al. 2004, Chrisman 2017), users of data that define spatial location are mostly using coordinates expressed in traditional planar coordinate systems, based on cartographic map projections that allow for quick calculations of the area of polygons bounded by straight line segments on projected ellipsoids. The pragmatic choice is an equal-area projection, such as the Albers equal-area conic projection – used in the USA (Snyder 1987), the Lambert Azimuthal Equal-Area projection – recommended in the EU for statistical analysis and display (Annoni et al. 2001) or even Mollweide projection – still in use for world or sky maps (Snyder 1987). Even if an equal-area projection is used, however, the polygon border segments, i.e. the projected geodesics, do not coincide with straight lines, which is an important source of error of the surface area calculated in this manner. Moreover, many of these projections cover only a part of the globe, which generates additional calculation problems when calculating areas outside these mapping zones or outside the boundaries of the projected ellipsoid surface.

As an introduction to the topic of calculating the surface area on an ellipsoid, let us analyse the task of calculation of the surface area of a triangle on a sphere of radius $R=$ 7⋅10⁶ = 7,000,000 m with the accuracy of 1m², using the formula $\mathrm{S}={R}^{2}\upvarepsilon$, i.e. $\mathrm{S}=4.91\cdot {10}^{13}\cdot \upvarepsilon$, where ε is the spherical excess. The errors are propagated according to the formula σ_s = 4.9⋅10¹³ · σ_ε [m²], σ_s = 1 [m²], hence σ_ε = 2⋅10^–14. Consequently, the angles must be determined with at least fourteen digits of precision. In contrast, the determination of the position $x=R \varphi =7\cdot {10}^{6} \varphi$ with the accuracy of σ_x = 1⋅10^–3 = 1 mm, means that the geographical coordinates need to be calculated with a precision of ten digits. In this context, the issue of the determination of the surface area of a polygon delimited by geodesics on an ellipsoid casts new light on the classical geodetic inverse problem and requires a decidedly higher computational precision.

The research presented here describes a new theory and approach to obtaining the explicit equation in cylindrical coordinates and the arc length of a geodesic on an ellipsoid of revolution and the surface area of a geodesic polygon, i.e. the area bounded by arcs of geodesics connecting vertices with known geographical coordinates such as Global Navigation Satellite System (GNSS) derived coordinates. The goal was to derive strict formulas and develop an algorithm that can quickly find the solution with high accuracy. Furthermore, the solution should be applicable to all cases, including near-antipodal and near-vertices cases.

The focus of this research includes the determination of the border path and surface area of a contour delineated by any feature on the surface of the Earth or any surface that is modelled by an ellipsoid of revolution. The feature’s boundaries and surface area are determined unambiguously, regardless of the mapping projection and the accuracy of the computed result is only limited by the floating-point representation of the computing system used. For example, high accuracy of the boundary and surface area is needed not only for the maintenance of a cadastre at the national level, but also at the international level or to support the calculation of charges associated with the use of a given contour, such as territorial seas and the leases for the exploitation of the shelf's natural resources.

The three most important novel contributions of this theory and approach to accurately calculate the surface area of a geodesic polygon are as follows.

• 1. A new geometry on an ellipsoid is introduced, namely an equatorial geodesic triangle (see Sect. 5.1). This surface shape, bounded by three geodesics, simplifies determination of the area under the geodesic and, consequently, the area of a geodesic polygon on an ellipsoid that equals the total areas under the geodesics stretched on each pair of successive polygon vertices.

• 2. The difference between the reverse and forward azimuths (${\alpha }_{12}$) of the arc of a geodesic is directly determined, based on Napier's rules, thus avoiding subtraction which reduces the accuracy of the surface area calculation (see Sect. 7.2).

• 3. An innovative parameterisation of the inverse problem with two alternative working unknowns (see Sect. 7.1) to solve all cases is presented.

We extend the work by Karney (2013), Danielsen (1989) and Vincenty (1975, 1976) in several ways:

• 1 A coherent original theory of ellipsoid geometry is presented; the novel, strict recursive formulas for the equation in cylindrical coordinates and the arc length of a geodesic, the area under the geodesic and the surface area of the geodesic polygon are derived; the area is determined by constant coefficients, not the values of finite series.

• 2 A cylindrical coordinate system is used and the integration is performed with respect to the radius of the circle of latitude, not to the frequently used reduced geodetic latitude.

• 3 The solution is applicable to any ellipsoid of revolution adopted as the reference surface, whether oblate or prolate.

• 4 The phenomenon of bifurcation, i.e. the existence of two different geodesic arcs of the same length connecting a given pair of points, leading to two solutions of the inverse geodetic problem, is detected.

• 5 A novel algorithm for solving the inverse problem, extended by the area under the geodesic calculation, is developed; it uses Newton's iterative process, and is based on the newly proposed parameterisation with two alternative working unknowns, providing a solution for all cases, including near-antipodal or near-vertices ones, and offers two solutions if bifurcation occurs.

• 6 A new formula is given to estimate the initial value of the geodesic C parameter, covering all cases, including near-antipodal or near-vertices ones, which replaces the fourth degree equation for Jacobi’s (1891) astroid frequently used.

• 7 The accuracy of the resulting values, i.e. the geodesic arc length, longitude and the area under the geodesic, are only limited by the computer hardware.

2 Related works

The two fundamental geodetic operations are the problems of the geodesic on an ellipsoid: the direct problem is to find the position along the geodesic given a starting point, an initial heading (azimuth) and a distance; while the inverse problem is to find the distance and the forward and reverse azimuths given the latitude and longitude of the two points. Although these problems are fundamentally important, only few researchers have attempted to solve them from a theoretical standpoint, i.e. to derive the formulas and fully examine the existence of a solution. The aim of the published formulas or methods for solving these problems was usually to offer specific practical applications (e.g. triangulation), which typically referred to geodesics of a few tens of kilometres. Long geodesic arcs, antipodal and near-antipodal cases were rarely analysed. This section presents a non-exhaustive selection of significant contributions in the research of the length of geodesics and area on an ellipsoid.

The issue of geodesics first appeared in the works of Clairaut (1733) and it was applied in geodesy in the form of the direct and inverse problem by Bessel (1826). The latter was the first to deal with long geodesics. Bessel developed a series expansion of elliptic integrals (distance and longitude) using an auxiliary sphere, while Helmert (1880) perfected the technology of the practical computation of the classical geodetic problems. Modifying the original Bessel method, Vincenty (1975, 1976) provided the formulas and algorithm for solving the problems that arose when technological progress brought the beginning of the computer era and improved accuracy of measurements. Although Vincenty did not present a derivation of the formulas, they have been published on the official website of the American National Geodetic Survey agency (www.ngs.noaa.gov) along with the original Fortran code, and they are widely known. Vincenty used the expansion of elliptical integrals into a series to modify Helmert’s iterative method and achieved high computational efficiency, which was extremely important for military applications, to which the formulas were then dedicated. Using Vincenty's formulas, however, a failure to converge for near-antipodal points was noticed, which is proven by a very high number of iterations. This fact was mentioned by numerous authors, e.g. Rapp (1993), Karney (2013) or Tseng (2014).

The studies by Kivioja (1971) mark the beginning of a new trend, namely numerical integration of elliptical integrals. Its supporters, e.g. Thomas and Featherstone (2005), developed the numerical integration of the differential equation of the geodesic. The parallel numerical approach was initiated by Saito (1970) and then developed by Sjöberg (2006, 2012). This involves the numerical integration of elliptical integrals, which are usually functions of the elements of the ellipsoid and the auxiliary sphere.

As far as the solutions for the longitude and distance along the geodesic are concerned, it is also worth noting the approaches using a recursion relation: by Levallois and Dupuy (1952) and by Pittman (1986). Pittman expressed the distance integral by a recursion formula consisting of base integrals. The delta lambda formula (longitude) involves double adding up (i.e. a series within a series). These are the only solutions that resemble some of the presented in this study, except that the recursion was used only for deriving formulas, not for calculations. Moreover, Pittman’s formula, in particular, is very complex in computational terms, and its assumption that the initial approximation of the parameter is of spherical origin, makes it impossible to solve antipodal tasks. Pittman’s significant contribution was to accelerate the convergence of the inverse geodetic problem’s solution by using Newton’s iterative process. Unfortunately, this method does not provide a solution if the end point of the geodesic arc is localised near its turning point, a problem of which Pittman is also aware. Moreover, Pittman's use of an approximate derivative increases the number of iterations; and this derivative equals zero if the end points of the geodesic arc lie on the same latitude.

The only formula to determine the surface area under a geodesic in geodetic coordinates on an ellipsoid was presented by Danielsen (1989), but this was not widely accepted at first. Danielsen derived a series to develop the formula to calculate the surface area under the geodesic and provided its six elements, each of which is a short finite series. The application of the formula requires the subtraction of two large, close numbers, calculated by a series, with a consequent sharp drop in the numerical accuracy of the result. Despite its imperfections, Danielsen's work represents a milestone in this field alongside that of Karney (2013), who popularised and extended it. Karney also used Newton’s iterative process to accelerate the inverse geodetic problem convergence.

3 Equations for geodesics on a surface of revolution

When $r,z,\lambda$, cylindrical coordinates are used on the surface of revolution, a curve is determined by the arc length differential and azimuth (Fig. 1) defined by formulas (1) and (2), respectively:

$$ {{ds}}^{2} = {{dm}}^{2} + {{r}}^{2} {{d}}\lambda ^{2}$$

(1)

$$ \sin \alpha = {{r}}\frac{{d\lambda }}{{ds}} $$

(2)

Fig. 1

Differential elements of a curve in cylindrical coordinates

where: r is the radius of a parallel (the distance of a point of the surface of revolution from the axis of rotation), $\lambda$ is the longitude (also referred to as geodetic longitude if datum is specified), m is the meridian arc length, α is the azimuth of a curve in a given point.

If a curve on a surface of revolution is a geodesic, it is characterised by Clairaut's relation (Clairaut 1733; Pressley 2010), also known as Clairaut's constant:

$$ r \sin \alpha = c = {\text{const}} $$

(3)

It should be noted that this relation discloses the sinusoidal character of geodesics on a surface of revolution and the ${\text{r}} \ge \left| {\text{c}} \right|$ constraint for the radius of a parallel. It can be demonstrated that the radius takes the minimum value, i.e. ${\text{r}}_{P} = \left| {\text{c}} \right|$, if and only if ${\text{P}} = {\text{P}}_{0}$ is a turning point of a geodesic.

Substituting (2) into (3) we get a new equation which determines the relationship between the longitude and the distance along a geodesic:

$$ {{r}}^{2} {{d}}\lambda = c\,{{ds}} $$

(4)

From (4) it is important to note that the sign of the constant c determines the sign of dλ.

Substituting ${\text{ds }}$ from Eq. (4) into Eq. (1) we get consecutively

$$ {{dm}}^{2} + {{r}}^{2} {{d}}\lambda ^{2} = \frac{{{{r}}^{4} }}{{{{c}}^{2} }}{{d}}\lambda ^{2} ,\;{{dm}}^{2} = \frac{{{{r}}^{2} \left( {{{r}}^{2} - {{c}}^{2} } \right)}}{{{{c}}^{2} }}{{d}}\lambda ^{2} $$

$$ {{d}}\lambda = \frac{{c{ }}}{{r\sqrt {{{r}}^{2} - {{c}}^{2} } }}{{dm }} $$

(5)

Substituting ${\text{d}}\lambda$ from Eq. (4) into Eq. (5) we get successively

$$ {{ds}} = \frac{{r{ }}}{{\sqrt {{{r}}^{2} - {{c}}^{2} } }}{{dm }} $$

(6)

After integrating Eq. (5) and (6), we obtain the explicit equation of a geodesic in cylindrical coordinates (7), i.e. longitude as a function of radius (often simply referred to as longitude), and the geodesic arc length (8) (often referred to as distance):

$$ {\lambda }\left( {c,{{r}}} \right) = {{c}}\int \frac{{{{r}}^{ - 1} {{dm}}}}{{\sqrt {{{r}}^{2} - c^{2} } }} $$

(7)

$$ {{s}}\left( {c,{{r}}} \right) = \int \frac{{{{r}}^{ + 1} {{dm}}}}{{\sqrt {{{r}}^{2} - c^{2} } }} $$

(8)

The similarity of both integrals was exposed in the notation, which consequently leads to a very efficient algorithm of their simultaneous determination.

4 Equations of geodesics on the ellipsoid of revolution

An ellipsoid of revolution in cylindrical coordinates is described with equation

$$ \frac{{r^{2} }}{{a^{2} }} + \frac{{z^{2} }}{{b^{2} }} = 1,~\quad - \pi < \lambda \le \pi ~ $$

(9)

Hence

$$ \frac{{{{dz}}}}{{{{dr}}}} = - { }\frac{{\sqrt {1 - e^{2} }\,r }}{{\sqrt {a^{2} - {{r}}^{2} } }}$$

$$ {{dm}} = \sqrt {1 + \left( \frac{dz}{dr} \right)^{2} } {{dr}} = - \sqrt {\frac{{a^{2} - e^{2} {{r}}^{2} }}{{a^{2} - {{r}}^{2} }}} {{ dr}} $$

(10)

where ${{e}}^{2} = 1 - \frac{{{{b}}^{2} }}{{{{a}}^{2} }} = {{f}}\left( {2 - {{f}}} \right), f = \frac{a - b}{a}$, a – equatorial semi-axis, b – polar semi-axis.

It should be noted that the parameter which is commonly marked as $e^{2}$ can have a negative value, i.e. any ellipsoid of revolution can be adopted as the reference surface, whether oblate or prolate.

Substituting Eq. (10) into Eq. (7) and Eq. (8) and introducing dimensionless values ${\uprho } = \frac{r}{a}, {\text{S}} = \frac{{\text{s}}}{a}, \mu = \frac{m}{a}, C = \frac{c}{a}$, we will get a description of the geodesic arc from the equator to the turning point (also referred to as the furthermost point of a geodesic), i.e. the longitude (11) and the distance (12) of the point of the radius $\rho = \rho_{P}$, respectively:

$$ {\uplambda }\left( {C;\rho_{P} } \right) = - C\mathop \int \limits_{1}^{{\rho_{P} }} \frac{{\sqrt {1 - e^{2} \rho^{2} } }}{{\sqrt {1 - \rho^{2} } }}\frac{{\rho^{ - 1} }}{{\sqrt {\rho^{2} - C^{2} } }}d\rho $$

(11)

$$ {\text{S}}\left( {C;\rho_{P} } \right) = - \mathop \int \limits_{1}^{{\rho_{P} }} \frac{{\sqrt {1 - e^{2} \rho^{2} } }}{{\sqrt {1 - \rho^{2} } }}\frac{{\rho^{ + 1} }}{{\sqrt {\rho^{2} - C^{2} } }}d\rho $$

(12)

Note that the following property occurs: $\left| C \right| \le \rho_{P} \le 1$.

4.1 Integrals describing the arc of a geodesic from the equator to the turning point

Introducing an auxiliary binomial series

$$ \begin{aligned} \left( {1 - e^{2} ~} \right)^{{1/2}} = & \mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i~~}}\!; \\ g_{{1,0}} = & 1;~~~~g_{{1,i > 0}} = \frac{{2i - 3}}{{2i}}e^{2} ~g_{{1,i - 1}} \end{aligned} $$

(13)

of which the constant coefficients are expressible by a Newton’s generalised binomial coefficient as $g_{1,i} = \left( {\begin{array}{*{20}c} {1/2} \\ i \\ \end{array} } \right)\left( { - e^{2} } \right)^{i}$, e.g. in Tseng (2014), however this requires more number of calculations than the proposed recursive formula (13).

We will convert (11) and (12) integrals into a linear combination of basic integrals

$$ \begin{aligned} \lambda \left( {C;\rho _{P} } \right) = & - C\mathop \int \limits_{1}^{{\rho _{P} }} \frac{{\rho ^{{ - 1}} \sqrt {1 - e^{2} \rho ^{2} } }}{{E\left( {C^{2} ;\rho ^{2} } \right)}}d\rho \\ = & - C\mathop \int \limits_{1}^{{\rho _{P} }} \frac{{\rho ^{{ - 1}} }}{{E\left( {C^{2} ;\rho ^{2} } \right)}}\mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i~~}} \rho ^{{2i}} d\rho = \\ = & - C\mathop \sum \limits_{{i = 0}}^{\infty } \mathop \int \limits_{1}^{{\rho _{P} }} \frac{{g_{{1,i~~}} \rho ^{{2i - 1}} d\rho }}{{E\left( {C^{2} ;\rho ^{2} } \right)}} = C\mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i~~}} G_{{1,i}} \left( {C;\rho _{P} } \right) \\ \end{aligned} $$

(14)

$$ \begin{aligned} S\left( {C;\rho _{P} } \right) = & - \mathop \int \limits_{1}^{{\rho _{P} }} \frac{{\rho ^{{ + 1}} \sqrt {1 - e^{2} \rho ^{2} } }}{{E\left( {C^{2} ;\rho ^{2} } \right)}}d\rho \\ = & - \mathop \int \limits_{1}^{{\rho _{P} }} \frac{{\rho ^{1} }}{{E\left( {C^{2} ;\rho ^{2} } \right)}}\mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i~~}} \rho ^{{2i + 1}} d\rho \\ = & - \mathop \sum \limits_{{i = 0}}^{\infty } \mathop \int \limits_{1}^{{\rho _{P} }} \frac{{g_{{1,i~~}} \rho ^{{2i + 1}} d\rho }}{{E\left( {C^{2} ;\rho ^{2} } \right)}} = \mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i~~}} G_{{1,i + 1}} \left( {C;\rho _{P} } \right) \\ \end{aligned} $$

(15)

where basic integrals are as follows:

$$ G_{{1,i}} \left( {C;\rho _{P} } \right) = - \mathop \int \limits_{1}^{{\rho _{P} }} \frac{{\rho ^{{2i - 1}} ~d\rho }}{{E\left( {C^{2} ;\rho ^{2} } \right)}} $$

(16)

$$ \begin{aligned} {\text{E}}\left( {C^{2} ;\rho ^{2} } \right) = & \sqrt {\left( {1 - \rho ^{2} } \right)\left( {\rho ^{2} - C^{2} } \right)} \\ = & \sqrt { - \rho ^{4} + \left( {1 + C^{2} } \right)\rho ^{2} - C^{2} } \\ \end{aligned} $$

(17)

Separating the derivative of the radicand, we will get

$$ \begin{aligned} G_{{1,i}} \left( {C;\rho _{P} } \right) = & \mathop \int \limits_{1}^{{\rho _{P} }} \frac{{ - 4\rho ^{3} + 2\left( {1 + C^{2} } \right)\rho }}{{2E\left( {C^{2} ;\rho ^{2} } \right)}}~~\frac{{\rho ^{{2i - 4}} }}{2}d\rho \\ - & \frac{{1 + C^{2} }}{2}\mathop \int \limits_{1}^{{\rho _{P} }} \frac{{\rho ^{{2i - 3}} d\rho }}{{E\left( {C^{2} ;\rho ^{2} } \right)}} \\ \end{aligned} $$

Integrating the first integrand by parts, we will get the following expression

$$ \begin{aligned} & G_{{1,i}} \left( {C;\rho _{P} } \right)\\ & \quad= \frac{{\rho _{P} ^{{2i - 4}} }}{2}{\text{E}}\left( {C^{2} ;~\rho _{P}^{2} ~} \right) \\ &\qquad - \frac{{2i - 4}}{2}\mathop \int \limits_{1}^{{\rho _{P} }} \rho ^{{2i - 5}} \frac{{ - \rho ^{4} + \left( {1 + C^{2} } \right)\rho ^{2} - C^{2} }}{{{\text{E}}\left( {C^{2} ;\rho ^{2} } \right)}}d\rho \\ &\quad= \frac{{\rho _{P}^{{2i - 4}} }}{2}{\text{E}}\left( {C^{2} ;\rho _{P}^{2} } \right) \\ &\qquad - \left( {i - 2} \right)\mathop \int \limits_{1}^{{\rho _{P} }} \frac{{ - \rho ^{{2i - 1}} + \left( {1 + C^{2} } \right)\rho ^{{2i - 3}} - C^{2} \rho ^{{2i - 5}} }}{{{\text{E}}\left( {C^{2} ;\rho ^{2} } \right)}}d\rho \\ &\quad= \frac{{\rho _{P}^{{2i - 4}} }}{2}{\text{E}}\left( {C^{2} ;\rho _{P}^{2} } \right) \\ & \qquad - {\text{~}}\left( {i - 2} \right)\left( {{\text{~G}}_{{1,i}} - \left( {1 + C^{2} } \right){\text{~G}}_{{1,i - 1}} + {\text{~}}C^{2} {\text{~G}}_{{1,i - 2}} } \right) \\ \end{aligned} $$

Now, the following functional equation can be written as

$$ \begin{aligned} {\text{G}}_{{1,i}} = & \frac{{\rho _{P}^{{2i - 4}} }}{2}{\text{E}}\left( {C^{2} ;\rho _{P}^{2} } \right) \\ - & {\text{~}}\left( {i - 2} \right)\left( {{\text{~G}}_{{1,i}} - \left( {1 + C^{2} } \right){\text{~G}}_{{1,i - 1}} + {\text{~}}C^{2} {\text{~G}}_{{1,i - 2}} } \right) \\ & \quad + \frac{{1 + C^{2} }}{2}{\text{G}}_{{1,i - 1}} \\ \end{aligned} $$

which, once rearranged, leads to a recursive relation enabling effective determination of basic integrals.

$$ {\text{G}}_{1,0} \left( {{\text{C}};\rho_{P} } \right) = \frac{1}{C}{\text{sin}}^{ - 1} \left( {\frac{{{\text{C}}\sqrt {1 - \rho_{P}^{2} } }}{{\rho_{P} \sqrt {1{ } - C^{2} } }}} \right) $$

(18)

$$ {\text{G}}_{1,1} \left( {{\text{C}};\rho_{P} } \right) = {\text{sin}}^{ - 1} \left( {\frac{{\sqrt {1 - \rho_{P}^{2} } }}{{\sqrt {1 - {\text{C}}^{2} } }}} \right){ } $$

(19)

$$\begin{aligned}& {\text{G}}_{{1,i > 1}}\\ & = \frac{{{\text{E}}\left( {C^{2} ;\rho _{P}^{2} } \right){\text{~}}\rho _{P} ^{{2i - 4}} + {\text{~}}\left( {2i - 3} \right)\left( {1 + C^{2} } \right){\text{~G}}_{{1,i - 1}} - {\text{~}}\left( {2i - 4} \right){\text{C}}^{2} {\text{~G}}_{{1,i - 2}} }}{{2i - 2}}\end{aligned} $$

(20)

The two initial basic integrals can be obtained by substituting $e^{2} = 0$ into formulas (11) and (12). They show the longitude (18) $C \cdot {\text{G}}_{{1,0}} \left( {{\text{C}};\rho _{P} } \right)\;\underline{\underline{{{\text{def}}}}} \;\omega$ and arc length (19) ${\text{G}}_{{1,1}} \left( {{\text{C}};\rho _{P} } \right)\;\underline{\underline{{{\text{def}}}}} \;\sigma$ of the geodesic on the sphere, which can be referred to as a special (boundary) case of an ellipsoid.

Since the basic integrals are positive and $0 \le \rho \le 1$, the integrand (16) cannot increase with an increase in the i index, and hence ${\text{G}}_{1,i + 1} \le {\text{G}}_{1,i}$. This, in turn, allows us to show the convergence of series (14) and (15) based on Eq. (13).

$$ {\uplambda }\left( {C;\rho_{P} } \right) \le C \cdot {\text{G}}_{1,0} \mathop \sum \limits_{i = 0}^{\infty } g_{1,i } = C \cdot {\text{G}}_{1,0} \left( {1 - e^{2} } \right)^{1/2} $$

$$ {\text{S}}\left( {C;\rho_{P} } \right) \le {\text{G}}_{1,1} \mathop \sum \limits_{i = 0}^{\infty } g_{1,i } = {\text{G}}_{1,1} \left( {1 - e^{2} } \right)^{1/2} $$

4.2 Decomposition of formulas which determine the longitude and distance of a geodesic

The recursive formula (18, 19, 20) suggests that the basic integrals (20) depend on the two initial integrals, i.e. (18, 19). Let us calculate the next integral:

$$ G_{1,2} = \frac{1}{2}E\left( {C^{2} ;\rho_{P}^{2} } \right) + \frac{{1 + C^{2} }}{2}G_{1,1} $$

It turns out that the third term is eliminated and this and successive integrals are functions with the following form

$$ G_{{1,i > 1}} = \gamma _{i} \left( {C^{2} } \right)~G_{{1,1}} + \varepsilon _{i} \left( {C^{2} ,\rho _{P}^{2} } \right)~E\left( {C^{2} ;\rho _{P}^{2} } \right) $$

(21)

The new coefficients are determined from the converted recursive formula (20)

$$ \begin{aligned} & \gamma _{0} = \;0,~\gamma _{1} = 1,\\ & \gamma _{{i > 1}} = ~\frac{{2i - 3}}{{2i - 2}}\left( {1 + C^{2} } \right)~\gamma _{{i - 1}} - ~\frac{{2i - 4}}{{2i - 2}}C^{2} ~\gamma _{{i - 2}} \end{aligned} $$

(22)

$$ \begin{aligned} &\varepsilon _{1} = 0,~~\varepsilon _{2} = \frac{1}{2}, \\ &\varepsilon _{{i > 2}} = \frac{{{\text{{\rm P}}}^{{2i - 4}} }}{{2i - 2}} + ~\frac{{2i - 3}}{{2i - 2}}\left( {1 + C^{2} } \right)~\varepsilon _{{i - 1}} - \frac{{2i - 4}}{{2i - 2}}C^{2} ~\varepsilon _{{i - 2}} \end{aligned} $$

(23)

Substituting (21) into formulas (14) and (15) we get the formulas which decompose integrals into spherical distance and position impacts (Eq. 24–27).

$$ \lambda \left( {C;\rho _{P} } \right) = \omega + C~\left( {\sigma ~Ls\left( {C^{2} } \right) + ~E\left( {C^{2} ;\rho _{P}^{2} } \right)~Le\left( {C^{2} ,\rho _{P}^{2} } \right)~} \right)\\ $$

(24)

$$ ~S\left( {C;\rho _{P} } \right) = \sigma ~Ss\left( {C^{2} } \right) + E\left( {C^{2} ;\rho _{P}^{2} } \right)~Se\left( {C^{2} ,\rho _{P}^{2} } \right)~ $$

(25)

$$ \begin{aligned} ~Ls\left( {C^{2} } \right) = \mathop \sum \limits_{{i = 1}}^{\infty } g_{{1,i}} ~\gamma _{i} ;~\quad Ss\left( {C^{2} } \right) = \mathop \sum \limits_{{i = 1}}^{\infty } g_{{1,i - 1}} ~\gamma _{i} \\ \end{aligned} $$

(26)

$$ \begin{aligned} Le\left( {C^{2} ,\rho _{P}^{2} } \right) = \mathop \sum \limits_{{i = 2}}^{\infty } g_{{1,i}} ~\varepsilon _{i} ;\quad Se\left( {C^{2} ,\rho _{P}^{2} } \right) = \mathop \sum \limits_{{i = 2}}^{\infty } g_{{1,i - 1}} ~\varepsilon _{i} \\ \end{aligned} $$

(27)

In this way, the elliptical integrals are expressed by combinations of elementary functions and polynomial series. Tseng (2014) obtained a similar character of the formulas for $\lambda$ and S.

The practical application of formulas (24) and (25) is facilitated by the properties of the multipliers. $Ls, Le, Se$ change in very narrow ranges (Table 1), while the $\gamma_{i} = 1$ (for $i > 1$) relation occurs on the equator (i.e. C = 1). Therefore, based on Eq. (13):

$$ \begin{aligned} Ss\left( {C^{2} = 1} \right) = & \mathop \sum \limits_{{i = 1}}^{\infty } g_{{1,i - 1}} \\ = & \mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i}} = \sqrt {1 - e^{2} } = 1 - f; \\ Ls\left( {C^{2} = 1} \right) = & \mathop \sum \limits_{{i = 1}}^{\infty } g_{{1,i}} = - f \\ \end{aligned} $$

Formulas (24), (25) express the effect of the change of the polar axis (in respect to the equatorial axis) on the deformation of the sphere. The main influence of the flattening $f$ propagates proportionally to the length $\sigma$. The remaining nonlinear part is proportional to the function ${\text{E}}\left( {C^{2} ;\rho_{P}^{2} } \right)$. It follows from Eq. (17) that the function E equals zero at the equator ${\text{E}}\left( {C^{2} ;\rho_{P}^{2} = 1} \right) = 0$ and at the turning points of a geodesic ${\text{E}}\left( {C^{2} ;\rho_{P}^{2} = C^{2} } \right) = 0$, while its maximum is

$$ \max {\text{E}}\left( {C^{2} ;\rho_{P}^{2} } \right) = {\text{E}}\left( {C^{2} ;\rho_{P}^{2} = \frac{{1 + C^{2} }}{2}} \right) = \frac{{1 - C^{2} }}{2}. $$

4.3 Solution in spherical coordinates

Section 4.1 presents integrals providing the solution for the arc of a geodesic from the equator to the turning point of a geodesic. Let us introduce new variables leading to the extension of the solution to the entire geodesic.

A sphere with geographical coordinates, marked by Karney as $\omega = \lambda^{{{\text{sfer}}}}$ and $\beta = \varphi^{{{\text{sfer}}}}$, is a boundary case of $e^{2} \to 0$ ellipsoid. In order to achieve full correspondence between the geodesics on the ellipsoid and the great circle on the sphere, we assume identical radii and azimuths, following Bessel (1826); hence

$$ \rho = {\text{cos}}\left( \beta \right) $$

(28)

where $\beta$ is the parametric (also referred to as reduced) latitude related to the geodetic latitude $\varphi$ by equation

$$ {\text{tan}}\left( \beta \right) = \left( {1 - f} \right) {\text{tan}}\left( \varphi \right) $$

(29)

The essence of extending the solution beyond the part of the arc of the geodesic line from the equator to its turning point, is the following substitution:

$$ \sqrt {1 - \rho^{2} } = \zeta = {\text{sin}}\left( \beta \right) $$

(30)

which allows to obtain the solution also on the southern hemisphere.

Using Clairaut's relationship, we get the extension

$$ \sqrt {\rho^{2} - C^{2} } = \rho \,{\text{cos}}\left( \alpha \right) $$

(31)

$$ {\text{E}}\left( {\alpha ,\beta } \right) = \zeta \,\rho\, {\text{cos}}\left( \alpha \right) = {\text{sin}}\left( \beta \right){\text{ cos}}\left( \beta \right)\, {\text{cos}}\left( \alpha \right) $$

(32)

that (a) has a geometric interpretation and (b) determines the signs of the decomposition components. The extension defined by formulas (31–32) is fundamentally important since it makes it possible to apply formulas (24) and (25) to any arcs of geodesic. Moreover, it replaces the artificial determination of combinations of signs of individual components that occurs in Pittman’s (1986) and Tseng’s (2014) methods.

Focusing more closely on the function E, which constitutes the nonlinear component of the longitude and the distance of the geodesic line defined by formulas (24) and (25), it is noted that it changes its sign when a geodesic crosses the equator (β = 0) or passes through its turning point (Fig. 2). This is a direct consequence of the properties of the turning point of a geodesic, which can be made apparent by differentiating relation (3)

$$ dr\;\sin (\alpha) + r\;\cos (\alpha) \;d\alpha = 0 $$

Fig. 2

Points of a geodesic representing some characteristic properties of the function E

hence

$$ \frac{{d{{r}}}}{{d{\upalpha }}} = - { }\frac{{{{r}}^{2} \cos {(\alpha) }}}{c}{ } $$

The sign of the cosine changes when passing through the point(s) ${\text{P}} = {\text{P}}_{0}$ where ${\upalpha } = \pm \pi /2$, which is precisely the turning point of a geodesic, thus ${\text{r}}_{P} = \left| {\text{c}} \right|$ ⇔ ${\uprho }_{P} = \left| {\text{C}} \right|$. Considering Eq. (28), β takes extreme values there, marking the northernmost or southernmost points of a geodesic.

Let examine the extreme values of the nonlinear terms in formulas (39) and (40). If $\rho_{1}^{2} = \rho_{2}^{2}$, then, according to formula (42), we have ${\text{Le}}_{1} = {\text{Le}}_{2} { };{\text{ Se}}_{1} = {\text{Se}}_{2}$, and from formula (17), it follows that $\left| {{\text{E}}_{1} } \right|$= $\left| {{\text{E}}_{2} } \right|$. As a result, the absolute values of the nonlinear terms of both formulas become equal, i.e. $|{\text{E}}_{1} {\text{ Le}}_{1} |= |{\text{E}}_{2} {\text{ Le}}_{2}|$, $|{\text{E}}_{1} {\text{ Se}}_{1}| = |{\text{E}}_{2} {\text{ Se}}_{2}|$. Let us also note that the signs of the function E are determined by the rule derived from Eq. (32). As a result, the maximum of nonlinear influences occurs when points 1, 2 are situated on the same parallel, while for points located on antipodes the influences disappear (Fig. 2).

The remaining relevant properties of the function E for the points 1 to 6 lying on the geodesic line illustrated on Fig. 2 are as follows: ${\text{E}}_{2} = {\text{E}}_{4} = {\text{E}}_{6} = 0$, ${\text{E}}_{3} = - {\text{E}}_{1}$, ${\text{E}}_{5} = {\text{E}}_{1}$.

4.4 Extension of the solutions to the entire geodesic

Using Eqs. (28–32) we reformulate Eqs. (18) and (19), i.e. the longitude and arc length of the geodesic on the sphere, in the following way

$$ {\text{ sin}}\left( {\omega } \right) = \frac{{{{C}} \,{\text{sin}}\left( \beta \right)}}{{\sqrt {1 - {{C}}^{2} } \,{\text{cos}}\left( \beta \right)}} $$

(33)

$$ {\text{sin}}\left( \sigma \right) = \frac{{{\text{sin}}\left( \beta \right)}}{{\sqrt {1 - {{C}}^{2} } }} $$

(34)

All solutions of the equations are determined by formulas for any integer k ($k \in {\mathbb{Z}}$):

$$ { }\sigma\left( {C,\beta ,k} \right) = k\pi + \left( { - 1} \right)^{k} {\text{ sin}}^{ - 1} \left( {\frac{{ {\text{sin}}\left( \beta \right)}}{{\sqrt {1 - {{C}}^{2} } }}} \right) $$

(35)

$$ \omega \left( {C,\beta ,k} \right) = k\pi + \left( { - 1} \right)^{k} {\text{sin}}^{ - 1} \left( {\frac{{{{C}}\, \,{\text{sin}}\left( \beta \right)}}{{\sqrt {1 - {{C}}^{2} } \,{\text{cos}}\left( \beta \right)}}} \right) $$

(36)

Substituting them into (24) and (25), we will get

$$ \begin{aligned} ~S\left( {C;\rho _{P} ,k} \right) &= \sigma \left( {C,\beta ,k} \right)~Ss\left( {C^{2} } \right) \\ &\quad + E\left( {C^{2} ;\rho _{P}^{2} } \right)~Se\left( {C^{2} ,\rho _{P}^{2} } \right)~~ \\ \end{aligned} $$

(37)

$$ \begin{aligned} & \lambda \left( {C;\rho _{P} ,k} \right) = \omega \left( {C,\beta ,k} \right) \\ &\quad + C~\left( {\sigma \left( {C,\beta ,k} \right)~Ls\left( {C^{2} } \right) + E\left( {C^{2} ;\rho _{P}^{2} } \right)~Le\left( {C^{2} ,\rho _{P}^{2} } \right)~} \right)~ \\ \end{aligned} $$

(38)

Hence, it can be concluded that a geodesic is infinite and, analogous to a sinusoid, it oscillates around the equator with a $2\pi \left( {1 + \left| C \right| Ls} \right)$ period, intersecting the equator in points $\beta = 0; \lambda_{k} = k\pi \left( {1 + \left| C \right| Ls} \right)$. In practical applications, the geodesics do not exceed half of the period and finding the solution for k = 0 or k = 1 is not usually problematic.

4.5 Geodesic arc of length not exceeding half of the period

An arc of the geodesic with C parameter connects two points with reduced latitudes $\beta_{1} ,\beta_{2}$ if $\left| C \right| \le \rho_{i} = {\text{cos}}\left( {\beta_{i} } \right)$ for i = 1,2. Then

$$ \begin{aligned} \lambda _{{12}} = & \;\lambda \left( {C;\beta _{2} } \right) - \lambda \left( {C;\beta _{1} } \right) \\ = & ~\;\omega _{{12}} + C~\left( {\sigma _{{12}} ~Ls\left( {C^{2} } \right) + ~E_{2} ~Le_{2} - E_{1} ~Le_{1} ~} \right)~~ \\ \end{aligned} $$

(39)

$$ \begin{aligned} S_{{12}} = & \;S\left( {C;\beta _{2} } \right) - S\left( {C;\beta _{1} } \right) \\ = & \;\sigma _{{12}} ~Ss\left( {C^{2} } \right) + E_{2} ~Se_{2} - E_{1} ~Se_{1} ~ \\ \end{aligned} $$

(40)

$$ \omega _{{12}} = \omega _{2} - \omega _{1} ;~~\sigma _{{12}} = \sigma _{2} - \sigma _{1} ;~ $$

(41)

$$ \begin{aligned} E_{i} &= E\left( {\alpha _{i} ,\beta _{i} } \right);~~ Le_{i} = Le\left( {C^{2} ,\rho _{i}^{2} ~_{~} } \right);~~ Se_{i} = Se\left( {C^{2} ,\rho _{i}^{2} ~} \right)~ \\ \end{aligned} $$

(42)

There are two possible cases: if $\sigma_{12} < \pi$ then the spherical triangle meets the requirements for the Euler triangle, which enables the use of spherical trigonometry, or if $\sigma_{12} = \pi$ then the spherical triangle degenerates to a biangle.

5 Area on an ellipsoid of revolution

Using the area differential on the ellipsoid ${d{\text{Area}}}={r \,d\lambda\, dm}={{a}}^{2}\,{\rho \,d\mu\, d\lambda }$ (Fig. 1), we determine the area distribution T along the meridian, which is the area (dA) measured from the equator to the parallel with radius ${\rho }_{P}$, between two differentially close ($d\lambda$) meridians

$$ d{\text{A}} = a^{2} ~{\text{T}}\left( \rho \right)~d\lambda $$

(43)

$$ ~{\text{T}}\left( {\rho _{P} } \right) = \mathop \int \limits_{1}^{{\rho _{P} }} \rho ~d\mu ~ = - \mathop \int \limits_{1}^{{\rho _{P} }} \frac{{\sqrt {1 - e^{2} \rho ^{2} } }}{{\sqrt {1 - \rho ^{2} } }}\rho ~~d\rho ~ $$

(44)

Using the binomial series (13) we transform the integral (44) into a following form:

$$ \begin{aligned} &{\text{T}}\left( {\rho _{P} } \right) = - \mathop \int \limits_{1}^{{\rho _{P} }} \frac{{\rho \sqrt {1 - e^{2} \rho ^{2} } }}{{\sqrt {1 - \rho ^{2} } }}~d\rho \\ & = - \mathop \int \limits_{1}^{{\rho _{P} }} \frac{\rho }{{\sqrt {1 - \rho ^{2} } }}~~\mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i}} \rho ^{{2i}} \,d\rho \\ & = - \mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i}} ~\mathop \int \limits_{1}^{{\rho _{P} }} \rho ^{{2i}} \frac{\rho }{{\sqrt {1 - \rho ^{2} } }}\,d\rho \\ \end{aligned} $$

In consequence, a linear combination of basic integrals is obtained as follows:

$$ {\text{ T}}\left( {\rho_{P} } \right) = \mathop \sum \limits_{{{\text{i}} = 0}}^{\infty } {\text{g}}_{1,i} {\text{ J}}_{i} \left( {\rho_{P} } \right) $$

(45)

$$ ~{\text{J}}_{i} \left( {\rho _{P} } \right) = \mathop \int \nolimits_{1}^{{\rho _{P} }} \rho ^{{2i}} \frac{{\left( { - \rho } \right)\,d\rho }}{{\sqrt {1 - \rho ^{2} } }} $$

(46)

Integrating by parts, we have

$$ \begin{aligned} {\text{J}}_{i} \left( {\rho _{P} } \right) &= \rho _{P}^{{2i}} ~\sqrt {1 - \rho _{P}^{2} } - 2i\mathop \int \limits_{1}^{{\rho _{P} }} \rho ^{{2i - 1}} \frac{{\left( {1 - \rho ^{2} } \right)~d\rho }}{{\sqrt {1 - \rho ^{2} } }} \\ & = \rho _{P}^{{2i}} \sqrt {1 - \rho _{P}^{2} } + 2i~{\text{J}}_{{i - 1}} \left( {\rho _{P} } \right) - 2i~{\text{J}}_{i} \left( {\rho _{P} } \right) \\ \end{aligned} $$

hence the recursive formula

$$ {\text{ J}}_{i > 0} \left( {\rho_{P} } \right) = \frac{1}{2i + 1}\left( {\rho_{P}^{2i} \sqrt {1 - \rho_{P}^{2} } + { }2i{\text{ J}}_{i - 1} \left( {\rho_{P} } \right)} \right) $$

(47)

$$ {\text{ J}}_{0} \left( {\rho_{P} } \right) = \sqrt {1 - \rho_{P}^{2} } $$

(48)

Since all these integrals depend on the initial one, any of these integrals can be written as

$$ {\text{ J}}_{i} \left( {\rho_{P} } \right) = {\text{J}}_{0} \left( {\rho_{P} } \right) \mathop \sum \limits_{k = 0}^{i} \tau_{i,k}\, \rho_{P}^{2k} $$

(49)

where

$$ \tau_{i,i} = \frac{1}{2i + 1}; \tau_{i,k < i} = \frac{2i}{{2i + 1}}\tau_{i - 1,k} $$

(50)

Substituting (49) into (45), we get

$$ {\text{ T}}\left( {\rho_{P} } \right) = \sqrt {1 - \rho_{P}^{2} } \mathop \sum \limits_{k = 0}^{\infty } t_{k} \,\rho_{P}^{2k} $$

(51)

$$ t_{k} = \mathop \sum \limits_{i = k}^{\infty } g_{1,i} \tau_{i,k} $$

(52)

Another method of calculating $t_{k}$ constants is worth analysing. Coefficients $\tau_{i,k}$ compose a triangular table where $0 \le k \le i; i \ge 0$. The k^th column of the table is necessary in order to compute $t_{k}$. Let us compare the coefficients of the k^th and k-1^st column in the same line. Using formula (50), we calculate successively

$$ \tau_{k,k} = \frac{1}{2k + 1}; \tau_{k,k - 1} = \frac{2k}{{\left( {2k - 1} \right)\left( {2k + 1} \right)}} $$

(53)

$$ \begin{aligned} ~&\tau _{{i,k}} = \frac{{2i}}{{2i + 1}}\tau _{{i - 1,k}} ; \\&\tau _{{i,k - 1}} = \frac{{2i}}{{2i + 1}}\tau _{{i - 1,k - 1}} \quad{\text{for}}\;i = k + 1,k + 2,...~~~ \\ \end{aligned} $$

Hence, after division

$$ \tau_{i,k - 1} = \tau_{i,k} \frac{{\tau_{i - 1,k - 1} }}{{\tau_{i - 1,k} }} \quad{\text{for}}\;i = k + 1,k + 2,... $$

By using the previous recursive relation we get

$$ \tau_{i,k - 1} = \tau_{i,k} \frac{{\tau_{k,k - 1} }}{{\tau_{k,k} }} \quad{\text{for}}\;i = k,k + 1,... $$

After substituting (53), this becomes

$$ \tau_{i,k - 1} = \tau_{i,k} \frac{2k}{{2k - 1}}\quad {\text{for}}\,\, i = k,k + 1,... $$

(54)

Substituting into (52) we get

$$ \begin{aligned} t_{{k - 1}} &= \mathop \sum \limits_{{i = k - 1}}^{\infty } g_{{1,i}} ~\tau _{{i,k - 1}} \\ &\quad= g_{{1,k - 1}} ~\tau _{{k - 1,k - 1}} + \mathop \sum \limits_{{i = k}}^{\infty } g_{{1,i}} ~\tau _{{i,k - 1}} \\ &\quad = \frac{{g_{{1,k - 1}} }}{{2k - 1}} + \frac{{2k}}{{2k - 1}}\mathop \sum \limits_{{i = k}}^{\infty } g_{{1,i}} ~\tau _{{i,k}} \\ &\quad= \frac{{g_{{1,k - 1}} }}{{2k - 1}} + \frac{{2k}}{{2k - 1}}t_{k} \\ \end{aligned} $$

Hence, the recursive formula which replaces the computation of series (52) is

$$ t_{k} = \frac{{\left( {2k - 1} \right)\, t_{k - 1} - g_{1,k - 1} }}{2k} $$

(55)

We calculate the initial constant as

$$ t_{0} = \mathop \sum \limits_{i = 0}^{\infty } g_{1,i} \,\tau_{i,0} = \mathop \sum \limits_{i = 0}^{\infty } \psi_{i} $$

After taking into consideration formulas (13) and (50), we get

$$ \psi_{0} = g_{1,0}\,\tau_{0,0} = 1 $$

$$ \begin{aligned} \psi _{{i > 0}} = & g_{{1,i}} ~\tau _{{i,0}} = \frac{{2i - 3}}{{2i}}e^{2} g_{{1,i - 1}} \frac{{2i}}{{2i + 1}}\tau _{{i - 1,0}} \\ = & \frac{{2i - 3}}{{2i + 1}}e^{2} \psi _{{i - 1}} \\ \end{aligned} $$

Extending the sequence for i > 3 results in the reduction in odd numbers in the numerator and denominator.

$$ \begin{aligned} \psi _{1} = & \frac{{ - 1}}{3}\psi _{0} = - \frac{1}{3},\; \\ \psi _{2}= & \frac{1}{5}\psi _{1} = - \frac{1}{{3~ \cdot ~5}},\; \\ \psi _{3} = & \frac{3}{7}\psi _{2} = - \frac{3}{{3~ \cdot ~5 \cdot ~7}} = - \frac{1}{{5~ \cdot 7}} \\ \end{aligned} $$

Hence, the final value of the constant

$$\begin{aligned} ~t_{0} &= 1 - \frac{{e^{2} }}{{1 \cdot 3}} - \frac{{e^{4} }}{{3 \cdot 5}} - \frac{{e^{6} }}{{5 \cdot 7}} - \ldots\\ & \quad - \frac{{e^{{2n}} }}{{\left( {2n - 1} \right)\left( {2n + 1} \right)}} - \ldots ~\end{aligned} $$

(56)

For the oblate ellipsoid ($e^{2} > 0$) the constant is determined with the well-known Mollweide’s formula (Mollweide 1805; Biernacki 1949, p.254;p.266)

$$ t_{0} = \frac{1}{2} + \frac{{\left( {1 - e^{2} } \right)}}{4e} {\text{ln}}\left( {\frac{1 + e}{{1 - e}}} \right) $$

and for the prolate ellipsoid ($e^{2} < 0$) the constant is determined with Wrinch’s (1932) formula

$$ t_{0} = \frac{1}{2} + \frac{{{\text{sin}}^{ - 1}\epsilon }}{{2\epsilon\sqrt {1 -\epsilon^{2} } }};{\text{ where }}\epsilon^{2} = \frac{{b^{2} - a^{2} }}{{b^{2} }} = \frac{{ - e^{2} }}{{1 - e^{2} }} $$

Obtaining the exact values of the constant coefficients $t_{k}$ is the key aspect of the theory of the area on the ellipsoid of revolution. The function ${\text{T}}\left( \rho \right)$ (defined at the beginning of Chapter 5) makes direct determination of the area of a geographic rectangle (A_rec) possible as:

$$\begin{aligned}& {\text{A}}_{{{\text{rec}}}} \left( {\rho_{1} < \rho < \rho_{2} , \lambda_{1} < \lambda < \lambda_{2} } \right)\\ & = a^{2} \left( {\lambda_{2} - \lambda_{1} } \right)\left( {{\text{T}}\left( {\rho_{2} } \right) - {\text{T}}\left( {\rho_{1} } \right)} \right)\end{aligned} $$

(57)

and, consequently, the area of an ellipsoid (denoted as A_ell) composed of two geographic rectangles

$$ \begin{aligned} {\text{A}}_{{{\text{ell}}}} = \, 2~{\text{A}}_{{{\text{rec}}}} \left( {0 < \rho < 1, - \pi < \lambda < \pi } \right) = \,4\,\pi \,a^{2}\,\left( {{\text{T}}\left( 0 \right) - {\text{T}}\left( 1 \right)} \right) = 4\,\pi\,a^{2} \,t_{0} \\ \end{aligned} $$

(58)

5.1 Area of the equatorial geodesic triangle

Let us reformulate Eq. (14) as

$$ ~d\lambda = \frac{{ - C}}{{\sqrt {1 - \rho ^{2} } ~\sqrt {\rho ^{2} - C^{2} } }}\mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i}} ~\rho ^{{2i - 1}} ~d\rho ~ $$

(59)

Substituting Eqs. (51) and (59) into (43), we will get the area of the equatorial geodesic triangle (denoted as A_egt), a geometry on an ellipsoid delimited by the geodesic, meridian and equator (e.g. shape 0,1,1’, Fig. 3)

$$ \begin{aligned} {\text{A}}_{{{\text{egt}}}}& \left( {C,\rho _{P} } \right) \\ = & - a^{2} C\mathop \int \limits_{1}^{{\rho _{P} }} \left( {\mathop \sum \limits_{{j = 0}}^{\infty } t_{j} \rho ^{{2j}} } \right) \left( {\mathop \sum \limits_{{i = 0}}^{\infty } g_{{1,i}} \rho ^{{2i - 1}} } \right)\frac{{d\rho }}{{\sqrt {\rho ^{2} - C^{2} } }} \\ = & \,a^{2} ~{\text{A}}\left( {C,\rho _{P} } \right) \\ \end{aligned} $$

(60)

Fig. 3

The area under the geodesic (hatched area) and the equatorial geodesic triangle (1,0,1’ or 2,0,2’ figure)

where ${\text{A}}\left( {C,\rho_{ P}} \right)$ is the normalised area of the equatorial geodesic triangle that we use subsequently.

After multiplying the two series in (60) and ordering using the variable $k = i + j$, we get

$$ {\text{A}}\left( {C,\rho } \right) = \mathop \sum \limits_{k = 0}^{\infty } \left( {\mathop \sum \limits_{j = 0}^{k} t_{j} g_{1,k - j} } \right)\left( { - C\mathop \int \limits_{1}^{{\rho_{P} }} \rho^{2k - 1} \frac{ d\rho }{{\sqrt {\rho^{2} - C^{2} } }}} \right) $$

(61)

Then, denoting

$$ g_{2,k } = \mathop \sum \limits_{j = 0}^{k} t_{j} g_{1,k - j} $$

(62)

$$ {\text{G}}_{2,k} \left( {C,\rho_{P} } \right) = - C\mathop \int \limits_{1}^{{\rho_{P} }} \rho^{2k - 1} \frac{ d\rho }{{\sqrt {\rho^{2} - C^{2} } }} = I_{k} \left( {\rho_{P} } \right) - I_{k} \left( 1 \right) $$

(63)

we get the area of interest as a linear combination of basic integrals

$$ {\text{A}}\left( {C,\rho_{P} } \right) = \mathop \sum \limits_{k = 0}^{\infty } g_{2,k } {\text{ G}}_{2,k} \left( {C,\rho_{P} } \right) $$

(64)

Let us calculate the integral

$$ I_{k} \left( \rho \right) = - C\int \rho^{2k - 1} \frac{ d\rho }{{\sqrt {\rho^{2} - C^{2} } }} $$

We convert the integrand to a form convenient for integration by parts (aiming at obtaining the derivative of the radicand in the numerator). We start with

$$ I_{k} \left( \rho \right) = I_{k} = - C\int \rho^{2k - 2} \frac{\rho \,d\rho }{{\sqrt {\rho^{2} - C^{2} } }} $$

Then, integrating by parts, we get

$$ \begin{aligned} I_{k} = & - C \,\rho^{2k - 2} \sqrt {\rho^{2} - C^{2} } + C \left( {2k - 2} \right)\int \rho^{2k - 3} \frac{{ \rho^{2} - C^{2} }}{{\sqrt {\rho^{2} - C^{2} } }}d\rho \\ = & - C \,\rho^{2k - 2} \sqrt {\rho^{2} - C^{2} } + C \left( {2k - 2} \right)\int \frac{{ \rho^{2k - 1} - C^{2} \rho^{2k - 3} }}{{\sqrt {\rho^{2} - C^{2} } }}d\rho \\ = & - C \,\rho^{2k - 2} \sqrt {\rho^{2} - C^{2} } + \left( {2k - 2} \right)\left( { - I_{k} + C^{2} I_{k - 1} } \right) \\ \end{aligned} $$

Hence, the recursive formula for the integral I_k is

$$ I_{k} = \frac{{ - C \,\rho^{2k - 2} \sqrt {\rho^{2} - C^{2} } + \left( {2k - 2} \right)C^{2} I_{k - 1} }}{2k - 1} $$

(65)

Since the integral $I_{1}$ does not depend on $I_{0}$ (i.e. $I_{1} = - C\sqrt {\rho^{2} - C^{2} }$), the subsequent integrals have the form

$$ I_{k > 0} = C\sqrt {\rho^{2} - C^{2} }\, \theta_{k} \left( {C^{2} , \rho^{2} } \right) $$

(66)

where the function ${\theta }_{k}$ is described by recursive formula

$$ \theta_{k} = \frac{{ - \rho^{2k - 2} + \left( {2k - 2} \right)C^{2} \,~\theta_{k - 1} }}{2k - 1} $$

(67)

Substituting (66) into (63), we get

$$\begin{aligned} {\text{G}}_{2,0} \left( {C,\rho_{P} } \right) & = - C\mathop \int \limits_{1}^{{\rho_{P} }} \frac{ d\rho }{{\rho \sqrt {\rho^{2} - C^{2} } }}\\ & = {\text{sin}}^{ - 1} \left( {\frac{C}{{\rho_{P} }}} \right) - {\text{sin}}^{ - 1} \left( C \right) = \alpha - \alpha_{0}\end{aligned} $$

$$ \begin{aligned} & {\text{G}}_{{2,k > 0}} \left( {C,\rho _{P} } \right) = \;I_{k} \left( {\rho _{P} } \right) - I_{k} \left( 1 \right) \\ &\quad = C\left( {\sqrt {\rho _{P}^{2} - C^{2} } ~~\theta _{k} \left( {C,\rho _{P} } \right) - \sqrt {1 - C^{2} } ~~\theta _{k} \left( {C,1} \right)} \right) \\ &\quad = {\text{sin}}\left( \alpha \right)~{\text{cos}}\left( \alpha \right)~\rho _{P}^{2} ~\theta _{k} \left( {C,\rho _{P} } \right) \\ &\qquad - {\text{sin}}\left( {\alpha _{0} } \right)~{\text{cos}}\left( {\alpha _{0} } \right)~\theta _{k} \left( {C,1} \right) \\ \end{aligned} $$

(68)

where $\alpha_{0} = {\text{sin}}^{ - 1} \left( C \right)$ is the azimuth under which a geodesic crosses the equator.

Substituting them subsequently into (64), we get

$$ \begin{aligned} {\text{A}}\left( {C,\rho _{P} } \right) = & \mathop \sum \limits_{{k = 0}}^{\infty } g_{{2,k~}} {\text{~G}}_{{2,k}} \left( {C,\rho _{P} } \right) \\ = & \;g_{{2,0~}} {\text{~G}}_{{2,0}} \left( {C,\rho _{P} } \right) + \mathop \sum \limits_{{k = 1}}^{\infty } g_{{2,k~}} {\text{~G}}_{{2,k}} \left( {C,\rho _{P} } \right) \\ = & \;g_{{2,0~}} \left( {\alpha - \alpha _{0} } \right) \\ + & \;{\text{sin}}\left( {2\alpha } \right)\frac{{\rho _{P}^{2} }}{2}\mathop \sum \limits_{{k = 1}}^{\infty } g_{{2,k~}} ~\theta _{k} \left( {C^{2} ,\rho _{P}^{2} } \right) \\ - & \;{\text{sin}}\left( {2\alpha _{0} } \right)\frac{1}{2}\mathop \sum \limits_{{k = 1}}^{\infty } g_{{2,k~}} ~\theta _{k} \left( {C^{2} ,1} \right) \\ \end{aligned} $$

And introducing the auxiliary function (Ar), which is the primitive function of the determined integral (60),

$$ {\text{Ar}}\left( {C^{2} ,\rho_{P}^{2} } \right) = \frac{{\rho_{P}^{2} }}{2} \mathop \sum \limits_{k = 1}^{\infty } g_{2,k } \,\theta_{k} \left( {C^{2} ,\rho_{P}^{2} } \right) $$

(69)

we will get the normalised area of the equatorial geodesic triangle

$$ \begin{aligned} {\text{A}}\left( {C,\rho _{P} } \right) = & g_{{2,0~}} \left( {\alpha - \alpha _{0} } \right) + \;{\text{sin}}\left( {2\alpha } \right)~{\text{Ar}}\left( {C^{2} ,\rho _{P}^{2} } \right) \\ &\quad - \;{\text{sin}}\left( {2\alpha _{0} } \right)~{\text{Ar}}\left( {C^{2} ,1} \right)~ \\ \end{aligned} $$

(70)

The use of the Ar function defined by the series (69), with fixed coefficients (62) and terms defined recursively (67), distinguishes formula (70) from the approximate formulas of Danielsen (1989) and Karney (2013) that use the coefficients defined by finite series.

The estimated values of constant coefficients (Table 2) indicate that the constants $g_{1,k} \approx e^{ - 2k}$. These eight coefficients are sufficient to determine the surface area of the Earth's ellipsoid (a² $\approx$ 4,068 ⋅10⁺¹³ [m²]) with an accuracy of 1 cm², as long as the resolution of the floating point implementation allows it.

Table 1 Variation range of the introduced $Ls,Ss, Le, Se$ multipliers

Since new formulas have been obtained, let us test them on the example of the surface area of the ellipsoid.

A biangle is a boundary case of an equatorial triangle. It is formed when the geodesic reaches the equator having crossed the turning point, then $\alpha = \pi - \alpha_{0}$. For the meridian $C = 0;{ }\alpha_{0} = 0;$ the biangle holds ¼ of the ellipsoid. This observation allows to check the calculations by computing the ellipsoid area using the following alternative formula

$$ {\text{A}}_{{{\text{ell}}}} = 4\,\pi\, g_{{2,0{ }}}\, a^{2} $$

(71)

Based on formula (62), we obtain $g_{{2,0{ }}} = t_{0} ,{ }g_{{1,0{ }}} = t_{0}$, which in turn confirms that the result is consistent with formula (58).

5.2 Area under the geodesic

Danielsen (1989) called the area between a geodesic and the equator and the two meridians connecting to the end points of the arc of geodesic (Fig. 3) the area under the geodesic.

It should be noted that the area under the geodesic connecting points 1 and 2 on an ellipsoid (${\text{A}}_{12}$) can be defined as the difference between the areas of two equatorial triangles, i.e. 2,0,2’ and 1,0,1’ (Fig. 3). Hence

$$ \begin{aligned} {\text{A}}_{{12}} = & {\text{A}}\left( {C,\beta _{2} } \right) - {\text{A}}\left( {C,\beta _{1} } \right) \\ = & g_{{2,0~}} \alpha _{{12}} + {\text{sin}}\left( {2\alpha _{2} } \right)\,{\text{Ar}}\left( {C^{2} ,\rho _{2}^{2} } \right) \\ & - {\text{sin}}\left( {2\alpha _{1} } \right)\,{\text{Ar}}\left( {C^{2} ,\rho _{1}^{2} } \right)~ \\ \end{aligned} $$

(72)

The sign of the area depends on the arc orientation, i.e. ${\text{A}}_{21} = - {\text{A}}_{12}$. The area is positive for arcs oriented according to the direction of the azimuth increase. The accuracy of the area calculation mainly depends on the accuracy of the difference between the reverse and forward azimuths $\alpha_{12} = { }\alpha_{2} - \alpha_{1}$ and therefore it should be calculated directly in the polar triangle.

Similarly to formulas (24) and (25), formula (72) shows the effect of the change of the polar axis on the deformation of the sphere, except that the main effect of the flattening $f$ is transferred proportionally by the coefficient $g_{2,0 }$. The azimuth difference $\alpha_{12}$ changes stepwise on the meridian, i.e. it takes the value π if points 1,2 are on either side of the pole, otherwise it equals 0. It follows that the pole causes problems in calculating the area.

5.3 Area of a geodesic polygon

The area of the geodesic polygon (A_gp) equals the total areas under the geodesics, which mark its border. Hence

$$ \begin{aligned} ~{\text{A}}_{{{\text{gp}}}} \left( {1,2,3, \ldots ,n,1} \right) = & \;\left( {{\text{A}}_{{12}} + {\text{A}}_{{23}} + \ldots + {\text{A}}_{{n1}} } \right)~a^{2} \\ & = \;a^{2} \mathop \sum \limits_{{i = 1}}^{n} {\text{A}}_{{i,i + 1}} ~ \\ \end{aligned} $$

(73)

where n stands for the number of polygon vertices, and in the last component of the sum, the polygon vertex $n + 1$ is equal to the first vertex ($n + 1\;\underline{\underline{{{\text{def}}}}} \;1)$. The surface area is positive if the polygon's border is travelled clockwise, i.e. according to the direction of the azimuth increase.

Let us analyse in detail the area of the geodesic triangle shown in Fig. 4. It can be seen from the figure that

Fig. 4

The area of a geodesic triangle (grey area)

$${\text{A}}_{{{\text{gp}}}} \left( {1,2,3,1} \right) = {\text{A}}_{{{\text{gp}}}} \left( {1,2,2{^{\prime}} ,1^{\prime},1} \right) + {\text{A}}_{{{\text{gp}}}} \left( {2,3,3^{\prime} ,2^{\prime},2} \right) - {\text{A}}_{{{\text{gp}}}} \left( {1,3,3^{\prime} ,1^{\prime},1} \right)$$

While, by changing the last area description to anticlockwise, we get

$$ \begin{aligned} & {\text{A}}_{{{\text{gp}}}} \left( {1,2,3,1} \right) = \;{\text{A}}_{{{\text{gp}}}} \left( {1,2,2^{\prime} ,1^{\prime},1} \right) \\ &\quad + \;{\text{A}}_{{{\text{gp}}}} \left( {2,3,3^{\prime} ,2^{\prime},2} \right) + {\text{A}}_{{{\text{gp}}}} \left( {3,1,1^{\prime}, 3^{\prime},3} \right) \\ & = \;\left( {A_{{12}} + A_{{23}} + A_{{31}} } \right)~a^{2} \\ \end{aligned} $$

which is in accordance with formula (73).

Worth analysing also is the case of the geodesic polygon delimiting the area, which includes the pole, namely a circumpolar polygon (Fig. 5). Applying formula (73) to the circumpolar polygon, we obtain the following negative result:

$$ \begin{aligned} {\text{A}}_{{{\text{gp}}}} \left( {1,2,3,4,1} \right) &= - {\text{A}}_{{{\text{gp}}}} \left( {1,2,2^{\prime} ,1^{\prime},1} \right) - {\text{A}}_{{{\text{gp}}}} \left( {2,3,3^{'} ,2^{\prime},2} \right)\\ &\quad - {\text{A}}_{{{\text{gp}}}} \left( {3,4,4^{\prime}, 3^{\prime},3} \right) - {\text{A}}_{{{\text{gp}}}} \left( {4,1,1^{'} ,4^{\prime},4} \right) \\ \end{aligned} $$

Fig. 5

A circumpolar polygon (grey area)

It is the only case when we get a negative result travelling along the polygon circumference according to the direction of the azimuth increase. As a matter of fact, this value applies to an area, which completes half of the ellipsoid area (the area marked in white in Fig. 5). While the area of the circumpolar polygon (${\text{A}}_{{{\text{cmp}}}}$) is defined as:

$$ {\text{A}}_{{{\text{cmp}}}} \left( {1,2,3, \ldots ,n,1} \right) = \frac{1}{2} {\text{A}}_{{{\text{ell}}}} + a^{2} \mathop \sum \limits_{i = 1}^{n} {\text{A}}_{i, i + 1} $$

(74)

6 Orthodromic path of a geodesic on an ellipsoid: The bifurcation phenomenon

Analysing the path of a geodesic connecting two equator points lying on the antipodes, one can see that two geodesics connect these two points, i.e. the equator and a meridian. On an oblate ellipsoid of revolution, the curve of the meridian is shorter and constitutes an orthodrome. The problem of the unequivocalness of a geodesic line, therefore, comes down to the question of when a geodesic constitutes an orthodrome.

6.1 Bifurcation of Clairaut’s equation. Bifurcate antipodal nodes

Considering that a geodesic is defined by Clairaut’s equation ($\rho \sin \alpha = {\text{C}}$), it can be noted that there are two corresponding geodesics with the same parameter C, passing through a given point with $\beta_{0} ,\lambda_{0}$ coordinates. Their azimuths are $\alpha_{1} = {\text{sin}}^{ - 1} \left( {{\text{C}}/\rho_{0} } \right), \alpha_{1^{\prime}} = \pi - \alpha_{1}$, respectively. According to the theory of differential equations (the azimuth is a function of the first derivative, see formula (2)), this phenomenon is termed bifurcation (Kuszta and Bailey 1982).

Let us analyse the path of geodesic lines in Fig. 6:

• i. The geodesics, with the initial azimuth $- \pi /2{ } < \alpha_{1} < \pi /2$, rise towards the north to reach the northern turning point for a minimal radius of $\rho = \left| {\text{C}} \right|$, after which they descend.

• ii. The other geodesics descend towards the south to reach the southern turning point for a minimal radius of $\rho = \left| {\text{C}} \right|$, after which they rise.

• iii. The corresponding geodesic lines intersect in the nodes lying on the on opposite parallel $( - \beta_{0} )$ on the other (antipodal) side of the equator plane; the curve of each of the lines corresponds to half a period.

Fig. 6

A bundle of orthodromic geodesics starting in the origin point creates a line of bifurcation nodes on the circle of latitude passing through the antipodal point. The Plate Carrée projection (normal cylindrical, equidistant along the meridians) was applied; the cylinder is developed by unrolling from a cut along the image of the meridian passing through the starting point of the bundle

Starting from the origin point, the geodetic arcs:

• i. smaller than half of the period, do not intersect and are orthodromes, because they uniquely determine the shortest path from the origin point to a given point of the ellipsoid;

• ii. equal to half a period, intersect at the nodes with their bifurcation counterparts, i.e. they constitute two arcs of the same length from the origin point to the node, and both are orthodromes;

• iii. greater than half a period, are not orthodromes because they always intersect with the above-described orthodromes covering the entire ellipsoid.

6.2 Bifurcation line of the nodes for an oblate ellipsoid

The nodes of all geodesics coming from a specific point form a bifurcation line assigned to the point, shown in Fig. 6 as the bold part of the parallel. The bifurcation line is the geodesic uniqueness border: the starting point and a node link two orthodromes with the same length and azimuths, which complete $\pi$; for instance $\alpha _{1} = 30^{{\text{o}}} ,~~\alpha _{2} = 150^{{\text{o}}}$ and $\alpha _{{1'}} = 150^{{\text{o}}} ,~~~\alpha _{{2'}} = 30^{{\text{o}}}$.

The beginning of the bifurcation line lies at the turning point of the geodesic with azimuths $\alpha _{1} = 90^{{\text{o}}} ,~~\alpha _{2} = 90^{{\text{o}}}$; and its parameter $C = \rho_{0}$. Substituting it into formula (39), we will get

$$ \lambda_{{{\text{min}}}} = \pi \left( {1 + \left| C \right| Ls\left( {C^{2} } \right)} \right) $$

Hence, the Boolean variable, which determines whether bifurcation occurs, is determined as:

$$ B\!I\!F = \left( {\beta _{2} = - \beta _{1} } \right)~{\text{and}}~\Big( {\left| {\lambda _{{12}} } \right| \ge \pi \left( {1 + \left| C \right|~Ls( {C^{2} } )} \right)} \Big) $$

(75)

Let us consider the case of the bifurcation $\sigma_{12} = \omega_{12} = \pi$ at the equator. Then $C = \pm 1, \alpha_{1} = \alpha_{2} = \pm \frac{\pi }{2}, E_{1} = E_{2} = 0, Ls = - f, Ss = 1 - f$. And from Eq. (40), we obtain $s_{12} = a S_{12} = a\left( {1 - f} \right)\pi = b\pi$. It follows that the equator is an orthodrome, provided its curve segment does not exceed the length of a semicircle with a radius equal to the semi-minor axis of an ellipsoid. This can be illustrated by imagining an ellipsoid in which the semi-minor axis decreases to zero. In such a case, the ellipsoid becomes a flat disc, and the shortest line between the points on the equator becomes a chord.

6.3 Bifurcation line of the nodes for prolate ellipsoid

In case the equator semi-axis is smaller than the polar semi-axis, the half-period is greater than $\pi$, since the parameter f for such an ellipsoid is negative, as can be seen from Eq. (10). Before reaching the parallel $- \beta_{0}$, the geodesic crosses the symmetrical geodesic with the azimuth $\alpha = - \alpha_{1}$. The intersection point lies on the antipodal meridian $\lambda_{0} + \pi$, which is a bifurcation line for the prolate ellipsoid ($e^{2} < 0$).

If the semi-axes are equal, the bifurcation line is reduced to the antipodal point, where all geodesics of a bundle meet.

7 The inverse problem extended with the area under the geodesic

The analysis of the problematic cases of the inverse problem, i.e. those involving the antipodes and the turning points of a geodesic, revealed issues related to the choice of parameters constituting the main equation, which manifest themselves as: (i) the lack of monotonicity of the function forming the equation; (ii) the consequent limitation of the applicability of Newton’s method, which requires the existence of a nonzero derivative; and (iii) the difficulty of determining the initial estimate for the iterative process.

The algorithms for solving the inverse problem are numerical solutions of an equation with one unknown. The classical methods of Bessel (1826) and their modifications, including the popular Vincenty's (1975) formulas, use $\omega_{12}$ as a working variable, while Karney (2013) assumes, after Helmert (1880), that the azimuth $\alpha_{1}$ constitutes a variable of the main equation. None of the methods used to date solves all cases of the inverse problem, and the reason for this is the selection of the controlled variable used in the iterations. The methods based on $\omega_{12}$ variable encounter problems in near-antipodal areas, while the methods based on $\alpha_{1}$, do not provide a solution if the case involves the area near the turning point of a geodesic. Therefore, we will revise the parameterisation of the indirect problem with the aim to obtain solutions in all cases. We will use two functionally connected variables for parameterisation. The inverse problem is symmetrical with respect to the end points, consequently, the algorithm should not prefer any of them.

In our algorithm, the initial data of our extended inverse geodetic problem are the same as for the classical inverse geodetic problem (i.e. latitude and longitude of two points: $\varphi_{1} ,\lambda_{1} ,\varphi_{2} ,\lambda_{2}$), while the outputs (i.e. $\alpha_{1} ,\alpha_{2} ,{ }s_{12} ,{\text{ Ar}}_{12} ,{ }T\!ype,{ }B\!I\!F$) additionally include the numerical value of the surface area under geodesic ($Ar_{12}$), the Type variable specifying the result of dynamic parameter choice, taking the value C or S (discussed further in ch.7.1) and a Boolean flag indicating the occurrence of bifurcation (BIF).

7.1 Analysis of the inverse problem. Novel parameterisation

Considering the above reasoning, we decided to base the parameterisation of the inverse geodetic problem on two variables, namely ${\upomega }_{12}$ and C. The former is Bessel’s (1826) well-known variable, and it represents the longitude on the auxiliary sphere (Fig. 7). It has been used by most authors of geodetic algorithms, including Vincenty (1975).

Fig. 7

Polar spherical triangle 1,2,P and path of the geodesic. The coordinate system of the sphere is defined by longitude $\omega$ and latitude $\beta$

In contrast, the latter variable (i.e. C) was introduced by Helmert (1880) as: $C = \rho_{{{\text{min}}}}$. The geometric interpretation of this variable is shown in Fig. 7 and is described by the following formula.

$$ C = \rho_{{{\text{min}}}} = \sin \alpha_{0} = \cos \beta_{{{\text{max}}}} $$

Karney (2013) and Pittman (1986) used variants of this variable. Karney replaced it by ${\upalpha }_{0}$ that $\alpha_{0} = {\text{sin}}^{ - 1} C$. The variable $\cos \beta_{{{\text{max}}}}$ was used by Pittman (1986), denoted as ${\text{cos}}\psi _{0}$.

The fundamental equation is composed of the given values difference (i.e. $\lambda_{12} = \lambda_{2} - \lambda_{1}$) and its value is expressed by formula (39)

$$ \delta = \lambda_{12} - \omega_{12} - C \left( {\sigma_{12} Ls\left( {C^{2} } \right) + E_{2} Le_{2} - E_{1} Le_{1} } \right) = 0 $$

(76)

Equation (76) is interpreted differently depending on the type of relations between the C and $\omega_{12}$ variables. If bifurcation occurs, then the spherical triangle (Fig. 7) degenerates into a digon; $\sigma_{12} = \left| {\omega_{12} } \right| = \pi$ and the main equation becomes a function of the $C$ parameter. The inverse problem has two solutions that meet Clairaut’s relation. Such a case is assigned the discriminant BIF = 1 in the output of the authors' implementations. Otherwise, the polar spherical triangle links C and $\omega_{12}$. Then, the inverse problem has one solution, and it is assigned the discriminant BIF = 0.

The fundamental Eq. (76) includes three variables $C,{ \omega }_{12}, {\sigma }_{12}$, which describe a polar spherical triangle (Fig. 7). In order to solve all cases, we need formulas which describe a biangle as a boundary case for triangle ${\sigma }_{12}=\left|{\omega }_{12}\right|=\pi$. We use a sine formula for polar spherical triangle converted into the form of Clairaut's equation’s

$$ \rho_{1} {\text{sin}}\left( {\alpha_{1} } \right) = \rho_{2} {\text{sin}}\left( {\alpha_{2} } \right) = \frac{{\rho_{1} \rho_{2} {\text{sin}}\left( {\omega_{12} } \right)}}{{{\text{sin}}\left( {\sigma_{12} } \right)}} = C $$

(77)

hence

$$ C {\text{sin}}\left( {\sigma_{12} } \right) - \rho_{1} \rho_{2} {\text{sin}}\left( {\omega_{12} } \right) = 0 $$

(78)

and cosine formula for the same triangle

$$ \cos \left( {\sigma_{12} } \right) = \zeta_{1} \zeta_{2} + \rho_{1} \rho_{2} {\text{cos}}\left( {\omega_{12} } \right) $$

(79)

An inverse problem will be solved iteratively with Newton's method. The creation of a system of equations for two variables $C, \omega_{12}$ requires differentiation of Eqs. (76), (78) and (79) eliminating $\sigma_{12}$. Having differentiated (79), we therefore get

$$ \sin \left( {\sigma_{12} } \right) d\sigma_{12} = \rho_{1} \rho_{2} {\text{sin}}\left( {\omega_{12} } \right) d\omega_{12} $$

After substituting (77), we will get

$$ d\sigma_{12} = C d\omega_{12} $$

(80)

This equation allows for elimination of $d\sigma_{12}$ from the other equations.

We develop the fundamental Eq. (76) into a Taylor’s series, neglecting the terms of the second and higher orders:

$$ \begin{aligned} \delta \left( {C + dC,\omega _{{12}} + d\omega _{{12}} } \right) & \cong \delta \left( {C,\omega _{{12}} } \right) \\ & \quad + \frac{{\partial \delta }}{{\partial C}}dC + \frac{{\partial \delta }}{{\partial \omega _{{12}} }}d\omega _{{12}} = 0 \\ \end{aligned} $$

Hence

$$ d_{1} dC + d_{2} d\omega_{12} = \delta $$

(81)

And the details of the differentiation are presented below. Let start with variable C.

$$ d_{1} = - \frac{\partial \delta }{{\partial C}} = \left( {Ls \sigma_{12} + Le_{2} E_{2} - Le_{1} E_{1} } \right) + C\sigma_{12} \frac{\partial Ls}{{\partial C}} $$

From formula (26) we get

$$ Ls\left( {C^{2} } \right) \cong { }g_{1,1} + g_{1,2} \frac{{1 + C^{2} }}{2} + g_{1,3} \frac{{3 + 2C^{2} + 3C^{4} }}{8} $$

$$ \frac{\partial Ls}{{\partial C}} = Cg_{1,2} + \frac{{\left( {C + 3C^{3} } \right)}}{2}g_{1,3} $$

Finally

$$ \begin{aligned} d_{1} & = \left( {Ls~\sigma _{{12}} + Le_{2} ~E_{2} - Le_{1} ~E_{1} } \right) \\ &\quad\, + C^{2} \sigma _{{12}} \left( {g_{{1,2}} + \frac{{\left( {1 + 3C^{2} } \right)}}{2}g_{{1,3}} } \right) \\ \end{aligned} $$

(82)

In turn, differentiating with respect to $\omega_{12}$, yields

$$ d_{2} = - \frac{\partial \delta }{{\partial \omega_{12} }} = 1 + C Ls\frac{{d\sigma_{12} }}{{d\omega_{12} }} $$

And after substituting (80), we get

$$ d_{2} = - \frac{\partial \delta }{{\partial \omega_{12} }} = 1 + C^{2} Ls $$

(83)

Derivative $d_{2}$ is close to zero and greater than $d_{1}$ by two orders of magnitude. It gives clear preferences towards variable $\omega_{12}$ and, consequently, allowed the use of direct iterations in Vincenty's method. Unfortunately, Vincenty’s approach turned out to have some shortcomings in respect to near-antipodal arcs.

Differentiating Eq. (78) and substituting (80), we will get

$$ h_{1} dC + h_{2} d\omega_{12} = 0 $$

(84)

$$ h_{1} = {\text{sin}}\left( {\sigma_{12} } \right) $$

(85)

$$ \begin{aligned} h_{2} &= \;C~{\text{cos}}\left( {\sigma _{{12}} } \right)\frac{{d\sigma _{{12}} }}{{d\omega _{{12}} }} - \rho _{1} \rho _{2} {\text{cos}}\left( {\omega _{{12}} } \right) \\ & = C^{2} ~{\text{cos}}\left( {\sigma _{{12}} } \right)- \rho _{1} \rho _{2} {\text{cos}}\left( {\omega _{{12}} } \right) \\ \end{aligned} $$

(86)

Equations (81) and (84) form a Newton's method system of equations for unknowns $dC,d\omega_{12}$.

In order to investigate Eq. (86), we will use formula (3)

$$ \begin{aligned} h_{2} &= C^{2} ~{\text{cos}}\left( {\sigma _{{12}} } \right) - \rho _{1} \rho _{2} ~{\text{cos}}\left( {\omega _{{12}} } \right) \\ & = \rho _{1} \rho _{2} \left( {{\text{sin}}\left( {\alpha _{1} } \right)~{\text{sin}}\left( {\alpha _{2} } \right)~{\text{cos}}\left( {\sigma _{{12}} } \right) - {\text{cos}}\left( {\omega _{{12}} } \right)} \right) \\ \end{aligned} $$

and the cosine formula for the angles of polar spherical triangle (Fig. 7) (${\text{cos}}\left( {\omega_{12} } \right) = {\text{cos}}\left( {\alpha_{1} } \right)\, {\text{cos}}\left( {\alpha_{2} } \right) + {\text{sin}}\left( {\alpha_{1} } \right){\text{ sin}}\left( {\alpha_{2} } \right) \,{\text{cos}}\left( {\sigma_{12} } \right)$), then

$$ h_{2} = - \rho_{1} \,\rho_{2} \,{\text{cos}}\left( {\alpha_{1} } \right) \,{\text{cos}}\left( {\alpha_{2} } \right) $$

(87)

Equation (84) has the following specific characteristics:

• i. Coefficient $h_{1}$ equals zero in the case of bifurcation, resulting in $d\omega_{12} = 0$;

• ii. Coefficient $h_{2}$ equals zero in the case when an end point of geodesic arc coincides with its turning point or when the geodesic arc coincides with the equator, resulting in $dC = 0$;

• iii. Both coefficients $h_{1}$, $h_{2}$ are equal to zero, i.e. the equation is identically satisfied for all values of $d\omega_{12}$ and $dC$, when both end points of geodesic arc coincide with its turning points, but this is an extreme case of bifurcation, i.e. $d\omega_{12} = 0$, and Eq. (81) is sufficient.

The solution of the system of Eqs. (81), (84)

$$ \begin{aligned} ~D =\, & d_{1} ~h_{2} - d_{2} ~h_{1} ; \\ {\text{if}}~D = \,&0~{\text{then}}~dC = \,\frac{{ - \delta }}{{d_{1} }}~{\text{else}}~dC = \,h_{2} \frac{\delta }{D}; \\ ~d\omega _{{12}} = \,&- h_{1} \frac{\delta }{D} \\ \end{aligned} $$

(88)

reveals a direct impact of the coefficients of homogeneous Eq. (83) on the values of the unknowns. In practice, only the unknown with the greater effective impact on the free term of Eq. (80) is determined.

$$ \left| {h_{1} d_{2} } \right| < \left| {h_{2} d_{1} } \right| $$

(89)

If inequality (89) is true, $dC$ is determined (and Type = C assigned to the solution) or, otherwise, $d\omega_{12}$ is determined (Type = S). The choice of the dynamic (adaptative) parameter determines the way the polar triangle is defined, and results in fewer iterations.

After inserting the approximate values of $d_{1} ,d_{2}$, the inequality takes the following form

$$ {\text{sin}}\left( {\sigma_{12} } \right) < \left| {h_{2} } \right| f \sigma_{12} $$

It is true for $\sigma_{12} \approx \pi$, then the estimation turns into

$$ \pi - \sigma_{12} < \left| {h_{2} } \right|f \sigma_{12} $$

hence

$$ \sigma_{12} > \frac{\pi }{{1 + \left| {h_{2} } \right|f}} $$

(90)

This means that the Type C solution is more effective in near-antipodal cases of an inverse problem.

The adaptative parameter choice was verified numerically. The numerical derivative of $\frac{{d\lambda_{12} }}{dC}$ was found to be equal to $\frac{\delta }{dC}$, and the numerical derivative of $\frac{{d\lambda_{12} }}{{d\omega_{12} }}$ was found to be equal to $\frac{\delta }{{d\omega_{12} }}$.

The Newton's iterative process is stopped when $\left| \delta \right| \le \varepsilon$, where $\varepsilon$ is an arbitrary stop-criterion. Then, in the last stage of our extended inverse problem, the distance and area are calculated according to formulas (40) and (73).

The (81) and (84) system of equations is an indirect solution of the complex theoretical problem of the reduced length of the geodesic that was analysed by Gauss (1828) and Christoffel (1868). For example, Helmert’s (1880) solution, later used by Karney (2013), requires laborious calculation of a specific series twice.

7.2 Two types of polar spherical triangles

The one type of polar spherical triangle (Type = S) is defined by the determined value of $\omega_{12}$ and the constant values $\beta_{1} ,\beta_{2}$. This is a classical problem of spherical trigonometry, and its solution will be based on Napier’s rules that enable direct determination of $\alpha_{12}$. If in any spherical triangle, two sides $a,c$ are given and the angle $B$ between them, then the cotangential formulas define the angles as follows:

$$ {\text{cot}}\left( A \right) = \frac{{{\text{cos}}\left( a \right)\, {\text{sin}}\left( c \right) - {\text{sin}}\left( a \right) \,{\text{cos}}\left( c \right)\, {\text{cos}}\left( B \right)}}{{{\text{sin}}\left( a \right) \,{\text{sin}}\left( B \right)}} = \frac{{q_{1} }}{{q_{2} }} $$

(91)

$$ {\text{cot}}\left( C \right) = \frac{{{\text{cos}}\left( c \right) \,{\text{sin}}\left( a \right) - {\text{sin}}\left( c \right)\, {\text{cos}}\left( a \right) \,{\text{cos}}\left( B \right)}}{{{\text{sin}}\left( c \right) \,{\text{sin}}\left( B \right)}} = \frac{{q_{3} }}{{q_{4} }} $$

(92)

It results from the sine formula that $q_{2} = {\text{sin}}\left( b \right) {\text{sin}}\left( A \right)$. While from the 'five-part rule' of spherical trigonometry, it results that $q_{1} = {\text{sin}}\left( b \right) {\text{cos}}\left( A \right)$ hence, ${\text{sin}}\left( b \right) = \sqrt {q_{1}^{2} + q_{2}^{2} } = \sqrt {q_{3}^{2} + q_{4}^{2} } = q_{5}$. This allows the calculation of the missing triangle side:

$${\text{tan}}\left( b \right) = \frac{{q_{5} }}{{{\text{cos}}\left( a \right)\, {\text{cos}}\left( c \right) + {\text{sin}}\left( a \right) \,{\text{sin}}\left( c \right) \,{\text{cos}}\left( B \right)}} $$

(93)

Using these formulas for a polar triangle, we get

$$ {\text{tan}}\left( {\sigma_{12} } \right) = \frac{{q_{5} }}{{\zeta_{1} \zeta_{2} + \rho_{1} \rho_{2}\, {\text{cos}}\left( {\omega_{12} } \right)}} $$

(94)

$$ \begin{aligned} q_{1} = & \;\rho _{2} \,{\text{sin}}\left( {\omega _{{12}} } \right),~~~~~ \\ q_{2} = & \;\rho _{1} \zeta _{2} - \rho _{2} \zeta _{1} \,{\text{cos}}\left( {\omega _{{12}} } \right),~~ \\ q_{5} = & \;\sqrt {q_{1}^{2} + q_{2}^{2} } \\ \end{aligned} $$

(95)

We get the parameter of the geodesic from modified formula (77)

$$ C = \frac{{\rho_{1} \rho_{2} \,{\text{sin}}\left( {\omega_{12} } \right)}}{{{\text{sin}}\left( {\sigma_{12} } \right)}} = \frac{{\rho_{1} \rho_{2} \,{\text{sin}}\left( {\omega_{12} } \right)}}{{q_{5} }} $$

(96)

While, the azimuths and their difference will be derived directly from Napier's rules:

$${\text{tan}}\left( {\frac{{\alpha_{1} + \alpha_{2} }}{2}} \right) = \frac{{{\text{cos}}\left( {\frac{{\beta_{1} + \beta_{2} }}{2}} \right)\,{\text{sin}}\left( {\frac{{\omega_{12} }}{2}} \right)}}{{{\text{sin}}\left( {\frac{{\beta_{2} - \beta_{1} }}{2}} \right)\,{\text{cos}}\left( {\frac{{\omega_{12} }}{2}} \right)}} $$

(97)

$$ {\text{tan}}\left( {\frac{{\alpha_{2} - \alpha_{1} }}{2}} \right) = \frac{{{\text{sin}}\left( {\frac{{\beta_{1} + \beta_{2} }}{2}} \right){\text{sin}}\left( {\frac{{\omega_{12} }}{2}} \right)}}{{{\text{cos}}\left( {\frac{{\beta_{2} - \beta_{1} }}{2}} \right){\text{cos}}\left( {\frac{{\omega_{12} }}{2}} \right)}} $$

(98)

The second type of the polar spherical triangle is defined by two sides and the C parameter of the third side. It is used for near-antipodal regions and for bifurcation, where the triangle shape is better defined by the C parameter. The case is assigned Type = C since formula (3) allows the azimuths to be calculated using ${\text{sin}}\left( {\alpha_{i} } \right) = \frac{C}{{\rho_{i} }}$.

In order to obtain an unambiguous solution, we check whether the polar triangle is obtuse. We analyse the cosine formula that determines the triangle side opposite to the azimuth

$$ \begin{aligned} {\text{cos}}\left( {\frac{\pi }{2} - \beta _{2} } \right) = & \;\zeta _{2} = \zeta _{1} ~{\text{cos}}\left( {~\sigma _{{12}} } \right) \\ + & \;{\text{sin}}\left( {\frac{\pi }{2} - \beta _{1} } \right)~{\text{sin}}\left( {~\sigma _{{12}} } \right)~{\text{cos}}\left( {\alpha _{1} } \right) \\ \end{aligned} $$

The second segment of the equation is negative, if ${\alpha }_{1}>\frac{\pi }{2}$. Moving the first segment of the equation to the left side, we obtain the sign of the cosine of the missing angle:

$$ m_{1} = \sqrt {\rho_{1}^{2} - C^{2} } ; \,{\text{if}} \, \zeta_{2} - \zeta_{1} {\text{cos}}\left( { \sigma_{12} } \right) < 0 \; {\text{then}} \; m_{1} = - m_{1} $$

(99)

$$ m_{2} = \sqrt {\rho_{2}^{2} - C^{2} } ;\, {\text{if}} \,\zeta_{1} - \zeta_{2} {\text{cos}}\left( { \sigma_{12} } \right) \ge 0 \; {\text{ then}} \; m_{2} = - m_{2} $$

(100)

Finally, the azimuths are determined from the formulas

$$ {\text{tan}}\left( {\alpha_{i} } \right) = \frac{C}{{m_{i} }} $$

(101)

and their difference from the direct formula

$$ {\text{tan}}\left( {\alpha_{2} - \alpha_{1} } \right) = \frac{{C\left( {m_{1} - m_{2} } \right)}}{{m_{1} m_{2} + C^{2} }} $$

(102)

Let us determine $\omega_{1}$ from formula (33)

$$ {\text{sin}}\left( {\omega_{1} } \right) = \frac{{C \zeta_{1} }}{{\sqrt {1 - {\text{C}}^{2} } \rho_{1} }} $$

$$ {\text{cos}}\left( {\omega_{1} } \right) = \frac{{\sqrt {\left( {1 - {\text{C}}^{2} } \right)\rho_{1}^{2} - C^{2} \zeta_{1}^{2} } }}{{\sqrt {1 - {\text{C}}^{2} } \rho_{1} }} = \frac{{\sqrt {\rho_{1}^{2} - C^{2} } }}{{\sqrt {1 - {\text{C}}^{2} } \rho_{1} }} $$

After substituting formula (99)

$$ {\text{cos}}\left( {\omega_{1} } \right) = \frac{{m_{1} }}{{\sqrt {1 - {\text{C}}^{2} } \rho_{1} }} $$

Finally, we get

$$ {\text{tan}}\left( {\omega_{12} } \right) = \frac{{C \left( {m_{1} \zeta_{2} - m_{2} \zeta_{1} } \right)}}{{m_{1} m_{2} + C^{2} \zeta_{1} \zeta_{2} }} $$

(103)

Hence, using formula (34)

$$ {\text{sin}}\left( {\sigma_{1} } \right) = \frac{{\zeta_{1} }}{{\sqrt {1 - {\text{C}}^{2} } }} $$

we get

$$ {\text{tan}}\left( {\sigma_{12} } \right) = \frac{{m_{1} \zeta_{2} - m_{2} \zeta_{1} }}{{m_{1} m_{2} + \zeta_{1} \zeta_{2} }} $$

(104)

7.3 Initial approximation

In order to obtain the starting point for Newton’s method, we shall begin with $\lambda_{12} = \lambda_{2} - \lambda_{1}$ normalisation so that $\left| {\lambda_{12} } \right| \le \pi$. It determines the signs of the following values ${\text{sign}}\left( C \right) = {\text{sign}}\left( {\omega_{12} } \right) = {\text{sign}}\left( {\lambda_{12} } \right)$. We calculate $C_{{{\text{max}}}} = {\text{min}}\left( {\rho_{1} ,\rho_{2} } \right)$, $Ls_{{\text{m}}} = Ls\left( {C_{{{\text{max}}}} } \right)$, which enables the determination of the BIF variable value from formula (75). If bifurcation occurs, then $\sigma_{12} = \pi$; $\omega_{12} = \pi {\text{ sign}}\left( {\lambda_{12} } \right)$, otherwise $\sigma_{12}$ is calculated from formula (93) assuming that $\omega_{12} = \lambda_{12}$.

The components of the fundamental Eq. (76) are estimated based on Table 1. Ls is the function of unknown C but it varies in very narrow ranges, and the extreme value Ls = $Ls_{{\text{m}}}$ is taken for the initiation. Le takes low values, and Le = 0 is assumed for the initiation. Thus substituting the fundamental equation simplified this way:

$$ \omega_{12} = \lambda_{12} - { }Ls_{{\text{m}}}\, { }\sigma_{12} { }\,C $$

(105)

Table 2 Estimated values of constant coefficients defined in Eqs. (13) and (62)

into (77), we get

$$ C{\text{ sin}}\left( {\sigma_{12} } \right) = \rho_{1} \rho_{2} \,{\text{sin}}\left( {\lambda_{12} - Ls_{{\text{m}}} { }\,\sigma_{12} { }\,C} \right) $$

Expanding the sine into a Taylor’s series, neglecting the terms of the second and higher order one gets the relation

$$ C{\text{ sin}}\left( {\sigma_{12} } \right) \cong \rho_{1} \rho_{2}\, {\text{sin}}\left( {\lambda_{12} } \right) - Ls_{{\text{m}}} { }\,\sigma_{12} \,\rho_{1} \rho_{2} \,{\text{ cos}}\left( {\lambda_{{12{ }}} } \right){ }\,C $$

Hence, the formula for the initial value of the C parameter is as follows:

$$ \tilde{C} = \frac{{\rho_{1} { }\rho_{2} {\text{ sin}}\left( {\lambda_{12} } \right)}}{{{\text{sin}}\left( {\sigma_{12} } \right) + Ls_{{\text{m}}} { }\,\sigma_{12} { }\,\rho_{1} { }\rho_{2} \,{\text{ cos}}\left( {\lambda_{12} } \right)}} $$

(106)

Although the approximation obtained is weaker than the 4^th degree equation of Jacobi’s (1891) astroid, the tests revealed that it does not affect the number of iterations (see Sect. 7.5). As a result, the proposed algorithm is characterised by relatively low computational complexity.

Using the value of the variable $C=\tilde{C }$ in formula (76), the approximation of the value of the second variable can be improved.

$$ \tilde{\omega }_{12} = \lambda_{12} - \tilde{C} \left( {\sigma_{12} \,Ls\left( {\tilde{C}^{2} } \right) + E_{2} \,Le_{2} - E_{1} \,Le_{1} } \right) $$

(107)

7.4 Algorithm for calculating the surface area of a geodesic polygon

This pseudocode outlines how to calculate the area of a geodesic polygon, given a list of the geographic coordinates of its vertices. FUNCTION AreaPolygon uses FUNCTION AreaArc to calculate the area under the geodesic.

7.5 Comments on numerical application

The ellipsoid area amounts to 5.1⋅10¹⁴. This means that in order to record only the numerical values of the area with 1m² precision, we have to be able to use numbers with 15 decimal digits (49 binary digits). The common IEEE 754 Double Floating Point Format PC standard has a 52-bit mantissa, which is sufficient to record values with such high precision. Still, making such precise real calculations is difficult, which is illustrated by the following numerical relation in double precision arithmetic, which contradict mathematics.

Example. Inequality $1 + \varepsilon > 1$ is not met by $0 \le \varepsilon \le 2^{ - 52} \cong 2.22$ ⋅10^–16.

7.6 Numerical tests

The method was verified by comparing the results obtained in IEEE 754 Double Floating arithmetic (52 bits of mantissa) with the results obtained in IEEE 754 Extended Floating Point arithmetic (63 bits of mantissa). The tests data and results are listed in Table 3. In the first column the geodetic latitudes of the given points and the longitude difference are displayed. The distance, the forward and reverse azimuths, i.e. the solution of the inverse problem, are displayed in the second column. The last column lists the value of the surface area under the geodesic obtained in double and extended arithmetic (denoted as ‘Area double’ and ‘Area extended’, respectively), the difference between these results (denoted as ‘diff’), the type specifying the result of dynamic parameter choice, taking the value C or S (as discussed in ch.7.1, and denoted as ‘Type’), bifurcation occurrence indicator (‘BIF’) and the number of iterations (‘Iter’).

Table 3 Numerical results of the extended inverse problem solved using the authors’ algorithm. Ist column: the inverse problem input data; IInd column: the results, i.e. distance and the forward and reverse azimuths; IIIrd column: the extension, where (Area) is the area under geodesics obtained using double and extended arithmetic, (diff) is their difference, (Type) is the parametrisation type, (Iter) is the number of iterations and (BIF) is the bifurcation occurrence indicator. Adopted ellipsoid: GRS80, a = 6378137 m, f = 298.257222101 (IAG and IAU recommended). Authors of test data: Karney (2013) – tests no. 3, 4 and 10; Nowak & Nowak Da Costa – remaining tests

The relative accuracy of the area ranges from 5⋅10^–14 for the longest arcs to 5⋅10^–12 for the shortest ones. The tests revealed that 1m² accuracy is achieved for distances under 10,000 km. It needs to be emphasised that the presented innovative iterative algorithm is characterised by a low number of iterations (Iter, Table 3), compared to the number of iterations when Rapp's (1993) method was applied. Consequently, a quick convergence of the proposed algorithm was achieved and proven.

We covered the calculation of the area of a geodesic polygon using the Bermuda Triangle as an example. The extended inverse problems listed in the last three rows of Table 3 (lines 11–13) allow the calculation of the area of a geodesic triangle: Miami $25^{ \circ } 47^\prime 16''{\text{N}},{\text{80}}^{ \circ } {\text{13'27"W}}$, Bermuda $32^\circ {20}{\mathrm{^{\prime}}}\mathrm{N}, 64^\circ 45\mathrm{^{\prime}}\mathrm{W}$, Puerto Rico $18^\circ {15}{\mathrm{^{\prime}}}\mathrm{N},66^\circ 30\mathrm{^{\prime}}\mathrm{W}$. The sum of the areas under successive geodesics equals 1,154,292,256,682 m².

Additionally, for the sake of comparison with the published results of the calculation of the area under the geodesic, i.e. 84,275,623.42235 km² given ${\varphi }_{1}=40^\circ$, ${\alpha }_{1}=30^\circ$ and ${s}_{12}=10000\mathrm{km}$ (Karney 2013), the calculations were made for ellipsoid WGS-84: a = 6,378,137, 1/f = 298.257223563 (adopted by Karney as the reference surface) obtaining 84,275,623,422,354.4844 m² in double arithmetic and 84,275,623,422,354.4531 m² in extended arithmetic.

7.7 Comment on existing methods

The main distinctive feature of the methods described for solving the inverse geodetic problem on arc 1–2 (Fig. 8) is the employment of different polar spherical triangles for the same purpose. Namely, Karney’s or Tseng’s method uses two triangles, i.e. (0, P, 1) and (0, 2, P); Pittman’s method uses both the (1, P, V) and (2, P, V) triangles; while our method, like Vincenty's, use only one, i.e. the target triangle, i.e. (1, P, 2).

Fig. 8

Polar spherical triangles used in different approaches to solve the inverse geodetic problem

In those methods that use two triangles, in each iteration, the series to determine the longitude and distance for the two arcs are computed. Subtraction of the resulting large, often similar, numbers is then performed, resulting in a decrease in accuracy (number of significant digits). Additionally, in Pittman's method, the same series are also calculated for the geodesic arc from the equator to the turning point. In Karney's method, two dedicated series are calculated and then subtracted to determine the reduced length, which directly increases the number of calculations but also reduces the number of iterations by three times compared to Tseng's method.

The second distinguishing feature of our proposed approach is the use of recursively computed base functions common to both series, which greatly improves the computation of these series. In addition, our method uses strict derivatives that reduces the number of iterations.

Finally, it is worth emphasising that all methods of solving the inverse geodetic problem use Clairaut's formula (1733), most of them in its classical form, i.e.

$$ \cos \;\beta \;\sin \alpha = \sin \alpha _{0} = C $$

Pittman’s and Tseng's methods, on the contrary, use the following form

$$ \cos \beta \sin \alpha = \cos \beta_{{{\text{turningPoint}}}} = C \ge 0 $$

which implicitly limits the azimuths to the interval $0 \le \alpha \le \pi$ (note that in case this inequality is not satisfied, the right side of the Clairaut formula can be negative, while its left side is always positive). Consequently, such an approach complicates the algorithms as it requires additional operations, such as sign analysis of individual components.

In contrast to our work, the expansion of integrals into series was implemented by computer algebra systems in the work of Karney (2013) and Tseng (2014). Our approach, though laborious and time-consuming, resulted in the possibility of recursion discovery and brought us to an efficient (in a mathematical sense) and effective (in a computational sense) solution.

8 Conclusions

In this paper, we aimed to develop a coherent theory of the ellipsoid’s geometry and accomplished our goal elaborating original, rigorous and implement convenient formulas for a geodesic polygon area, the longitude and the arc length of a geodesic based directly on geographical coordinates, e.g. from GNSS. While there were a number of partial approaches focusing on practical solutions for a longitude and distance along a geodesic for which approximate solutions were sufficient, we are unaware of any previous work that simultaneously attains the above-mentioned threefold goal.

In deriving the elliptical integral formulas, a system of cylindrical coordinates was used, whereby integration was carried out in reference to the normalised radius of parallel ρ, and not in reference to the reduced latitude, as most authors have done since Bessel. When developing the integrals for the geodesic longitude and arc length into a series, their unique similarity was used, i.e. it helped to use the same basic functions. A direct consequence of this method, in addition to being convenient to apply, is the decomposition of a series into two parts: depending on the distance and on the azimuth. This decomposition makes it possible to solve only one spherical triangle instead of two, as used by Danielsen (1989) and Karney (2013). This contributes indirectly to solving the issue of function periodicity, with which some scientists, such as Bessel (1826) and Pittman (1986), struggled, because the distance-dependent part of the formula is interconnected with periodicity.

The other significant objective of this study was to provide the formulas implementation algorithm to effectively obtain results with an accuracy that is limited only by the computer arithmetic applied. This objective was achieved though (a) formulas, defined with the use of constant factors and not with finite series values, as in previous formulas (cf. Danielsen 1989), as well as (b) the use of Napier's rules for the direct determination of the difference between the forward and reverse azimuths, and consequently (c) the elimination of their subtraction that reduces the accuracy of the surface area calculation.

The provided formulas are of recursive form that facilitates and simplifies their coding, which further contributes to a greater efficiency of our algorithm seizing Newton's rapidly converging iterative process and a new simple method to determine the geodesic C parameter’s initial value for all cases, including those that previously required solving a 4^th degree equation describing Jacobi's (1891) astroid.

Newton's iterative process was constructed for two alternative parameters, C and ω. It is an indirect solution to a difficult theoretical problem of a reduced length of the geodesic investigated by Gauss (1828) and Christoffel (1868), and it replaces Helmert's solution (1880), which requires labour-consuming calculation of a special series. The proposed parameterisation with two alternative working unknowns includes all cases and, in particular, the cases when derivatives do not exist: near the turning point of a geodesic, where Karney's (2013) and Tseng’s (2014) algorithm lose their convergence, and in the near-antipodal case, where Vincenty's (1975) algorithm diverges.

The orthodromic characteristics of a geodesic were analysed, which helped to detect the bifurcation lines on which there are two orthodromically equivalent solutions of a geodetic inverse problem, offered by the algorithm presented. New solutions were obtained of a polar geodesic triangle for near-antipodal and bifurcation cases defined by the C parameter value.

The algorithm presented was tested and verified using numerical examples from known publications and our original ones. The numerical tests, performed in extended and double arithmetic, proved the accuracy of 1m² of the computed geodesic area using standard IEEE 754 Double Floating Point arithmetic, for geodesic lengths up to 10,000 km.

The method described converges rapidly, despite the fact that the accuracy of the solution of the inverse geodetic problem 1⋅10^–15, necessary to calculate the area, typically requires additional iteration. The number of iterations, like accuracy, depends on the length of the geodesic (cf. Tab.3): for short arcs, dominant in the applications, only one iteration is sufficient to obtain the final result, with the relative accuracy of the surface area being 5⋅10^–12; for long arcs, on the other hand, up to three iterations are sufficient and the relative accuracy is 5⋅10^–14.

Summing up, we have achieved our goal of developing strict formulas, which are convenient for use in applications and implementations, and an innovative algorithm, which is characterised by its relative simplicity (executable in only a few steps and computationally less complex), its rapid convergence, its high accuracy (depending only on the floating point implementation), and its universal character (applicable for any ellipsoid of revolution and all cases, including near-antipodal cases, cases near the turning point of a geodesic or bifurcation case).