Hostname: page-component-745bb68f8f-v2bm5 Total loading time: 0 Render date: 2025-02-06T10:02:12.508Z Has data issue: false hasContentIssue false
  • English
  • Français

Accurate Harmonic Series for Inverse and Direct Solutions for the Great Ellipse

Published online by Cambridge University Press:  07 June 2011

Michael A. Earle
Rights & Permissions [Opens in a new window]

Abstract

In this paper, formulae drawn from the discipline of map projections are applied to provide simple and accurate solutions for the Inverse and Direct problems on the Great Ellipse. Distance along the meridional arc of the spheroid as a function of geodetic latitude is defined in terms of an elliptic integral which will be replaced here with a compact harmonic series approximation possessing simplicity and high accuracy. Latitude as a function of distance along the meridional arc will also be obtained via another equally simple inversion series that also possesses high accuracy. When these two series are applied with their constants modified to suit the section ellipse, they will be shown to provide accurate solutions to the inverse and direct navigation problems pertaining to the Great Ellipse, and thereby provide a complete solution that is also simple to implement.

Keywords

NavigationGreat EllipseGreat CircleInversion
Type
Research Article
Information
The Journal of Navigation , Volume 64 , Issue 3 , July 2011 , pp. 557 - 570
Copyright
Copyright © The Royal Institute of Navigation 2011

1. INTRODUCTION

The modern Global Reference Surface (GRS) is a spheroid generated by the WGS 84 reference ellipse rotated about its minor axis. It enables navigational track solutions on the Earth's surface to be determined that are more precise than those that would be obtained if the Earth was considered as a true sphere. The intersection of the surface of the spheroid with a plane passing through its poles produces a meridian ellipse. At any other inclination away from the poles, the intersection through the centre produces a section ellipse and the curve at the intersection is a Great Ellipse (GE). Similarly, the intersection of a sphere and a plane through its centre creates a Great Circle (GC). Mathematics for GC navigation on the standard navigation sphere have worked well as the basis for navigation; it has the advantage of being well established and practiced throughout the maritime community and is also valuable in an instructional setting. Formulae for GC computation are based on the properties of spherical triangles that enable solutions to both the inverse problem and the direct or forward problem to be obtained with ease. A limitation of GC navigation is that distances on the sphere when compared to those obtained for the GE on the spheroid can differ by 0·5% at worst but are typically less (~0·2%) depending upon position (Earle, Reference Earle2006). Differences between azimuth directions on the sphere and the spheroid are known to be quite small and will be further remarked upon later.

Navigation is now aided by satellite positioning systems (SPS) and facilitated by Electronic Chart Display and Information Systems (EDS and ECDIS) which, taken together, are the new norm. These systems provide high positioning accuracy and have become critical to the planning and execution of voyages. Uncertainty of position derived from these systems is now so small as to be on the order of, and frequently less than, the dimensions of the vessel being steered and so computational tools commensurate with this level of positioning accuracy are required. These requirements are met with the mathematics of the GE as it applies to the WGS 84 reference datum, which is a common default datum for GPS receivers, but are not adequately met by use of GC computation.

Transoceanic route planning usually entails establishing an intended GE track from a set of geographic endpoint coordinates together with coordinates for a number of intermediate geographic waypoints along the GE. Overall distance is then computed, along with distance and direction between waypoints. Steering between waypoints is done as a series of rhumb line legs chosen to approximate the desired GE track. In practice, waypoint selection is an operational matter that may at times require track adjustments to maintain overall operational efficiency. This may be achieved by adjusting the aspect of the vessel to the prevailing weather and sea state conditions, taking into account the hull form and freeboard of the vessel. These perturbations to the ideal track are outside the scope of this paper though some sense of their significance to track optimization has been discussed (Martinez de Osés, Reference Martinez2005).

Modern ship bridge automation and technology has simplified much of the effort in route planning and execution. The process of defining the end points of an intended track and then determining the overall distance is known as the inverse problem. On the other hand, the direct or forward problem specifies a point of departure and an initial azimuth and determines arrival coordinates after travelling a specified distance along a GE track. At first sight, even though the GE can be established from a single geographic position, such a solution might be considered to be too open-ended for the safe conduct of a vessel at sea and may be of little interest to the navigator. The direct solution benefits the navigator if the practice of route planning requires waypoints to be set at predetermined separations in distance along an intended GE line between a departure point and a pre-established destination point. If however, waypoints are to be separated by divisions in longitude, then having developed the GE from two geographic positions, the inverse solution provides all the required data. Mathematical schemes that provide both the inverse and the forward solution for the GE are known as complete solutions, for which there is a new algorithm (Pallikaris and Latsas, 2009), that provides a consolidated approach that draws upon the work of Bowring (Bowring, Reference Bowring1984).

Calculating distance between two points on a great ellipse can be done via a solution to a line integral using numerical methods. This requires a continuous mathematical relationship between latitude ψ (or ϕ) and longitude θ to be defined via what is called the Equation of the Great Ellipse (EGE) that will be described later. Although precise numerical evaluation employing an adaptive algorithm is quite straightforward, it is also possible to calculate distance by an approximate series formula that possesses very good accuracy. An appeal of the series solution is that software preparation is simplified somewhat. In the following sections, one such harmonic series approximation for distance along the meridian ellipse as a function of geodetic latitude is discussed and results obtained from it compared to results from the integral method. After that, the properties of a companion harmonic series inversion formula that allows distance along the meridian ellipse to be converted back to geodetic latitude will be discussed. In the sections that follow on, these two series expressions for the meridian ellipse will be modified for application to the GE on the section ellipse. This will be done by changing the constant terms in the two series to account for the reduced ellipticity of the section ellipse. The section ellipse is then the GE on which distance is to be calculated as a function of an elliptic angle. In the inversion process, distance along the ellipse will be converted to an elliptic angle and then into geodetic coordinates. The combination of these two schemes provides the complete solution.

2. FUNDAMENTAL SOLUTION FOR GE DISTANCE

The square of a line element of distance ds on the surface of a spheroid is established from a vector of position X(φ, θ) and written as:

(1)
ds^{\setnum{2}} \equals \left\vert {{\bi X}_\phi } \right\vert^{\setnum{2}} \lpar d\phi \rpar ^{\setnum{2}}\plus \left\vert {{\bi X}_\theta} \right\vert^{\setnum{2}} \lpar d\theta \rpar ^{\setnum{2}}

This is known as the first fundamental form of X=X(φ, θ) for which, Xφ=∂X/∂φ and Xθ=∂X/∂θ whose magnitudes are known as the fundamental coefficients of the first order (Lipshutz, Reference Lipshutz1969), (Massey and Kestelman, Reference Massey and Kestelman1958). Equation 1 is applicable to any regular curve on the surface of the ellipsoid; it is definitive and contains no assumptions; consequently any solution to Equation 1 is therefore a fundamental solution and is exact. One such line integral solution to Equation 1 (Earle, Reference Earle2000) based upon the EGE has been provided for the GE and is recast here as:

(2)
s \equals a\alpha \int\limits_{\theta _{\setnum{1}} }^{\theta _{\setnum{1}} } \left[ {{1 \over {1 \plus \alpha \tan \psi ^{\setnum{2}} }}} \left( {1 \plus {{\lpar 1 \plus \tan \psi ^{\setnum{2}} \rpar \lpar \lambda \sin \theta \minus \mu \cos \theta \rpar ^{\setnum{2}} } \over {\lpar 1 \plus \alpha \tan \psi ^{\setnum{2}} \rpar ^{\setnum{2}} }}} \right)\right]^{1 \sol 2} d\theta

In this equation, a=10800/π is the equatorial radius in geodetic miles for which one geodetic mile is equal to 1855·324847 metres or 1·001795 nautical miles, |s|=distance along the GE in the same units, while α=1−ε2 where ε is the ellipticity of the meridian ellipse. Latitude ψ is a function of longitude defined by the EGE of Equation 10, shown later. Although the above expression is exact, it must be solved using numerical methods preferably with an adaptive algorithm. Results obtained later from Equation 2 will be used to establish a baseline against which other results will be compared. When distance has been specified, forward solutions for longitude θ can be determined from Equation 2 by means of iteration methods, after which, latitude can be determined on substituting the solved value of θ into the EGE. In the developments that follow, iteration methods are avoided. For the determination of forward and backward azimuths associated with each position on the spheroid, some of the solutions previously presented (Earle, Reference Earle2008) can be readily applied.

3. THE MERIDIAN ELLIPSE

The WGS 84 meridional ellipse has an ellipticity ε=0·081819191 which when rotated defines the GRS. A quarter section through the meridian ellipse is illustrated in Figure 1 where OE=equatorial radius a, ON is the polar axis and r=OP is the radius at P. The position of a point P on the surface of a spheroid is located by latitude and longitude coordinates. Latitude is specified north or south of the equator as geodetic latitude ψ or geocentric latitude ϕ, while longitude θ is specified east or west of the prime meridian. The geodetic latitude ψ, sometimes called geographic latitude, is the angle between the inward projected normal to the surface tangent at a point P and the equatorial plane.

Figure 1.

Angles φ and ψ are related by tan φ=α tan ψ. Distance on the meridian ellipse from the equator at E to a point at P at latitude ψ is given in geodetic miles by:

(3)
F\lpar \psi \rpar \equals a\alpha \int\limits_{\setnum{0}}^{\psi } {{{d\xi } \over {\sqrt {\lpar 1 \minus \varepsilon ^{\setnum{2}} \sin \xi ^{\setnum{2}} \rpar ^{\setnum{3}} } }}}

Equation 3 is easily evaluated numerically and even elementary methods such as Simpson's rule will work but may not have sufficient precision, although an algorithm described in (Williams, Reference Williams1998) is known to work well. It is preferable however, to use an adaptive algorithm that adjusts the intervals of the integrand according to the slope of the function. Overall distance M 12 between two latitudes on the meridional arc in the same hemisphere is that from:

M_{\setnum{12}} \equals F\lpar \psi _{\setnum{2}} \rpar \minus F\lpar \psi _{\setnum{1}} \rpar

3.1. Series Approximations for the Meridian Ellipse

The function f 1(ψ) below is a compact harmonic series approximation to Equation 3 for meridional distance (Snyder, Reference Snyder1987).

(4)
f_{\setnum{1}} \lpar \psi \rpar \equals a\left[ {a_{\setnum{0}} \psi \plus \sum\limits_{n \equals \setnum{1}}^{\setnum{3}} {a_{n} \ \sin \lpar 2n\psi \rpar } } \right]

The coefficients are:

\eqalign{ a_{\setnum{0}} \equals \tab 1 \minus {1 \over 4}\varepsilon ^{\setnum{2}} \minus {3 \over {64}}\varepsilon ^{\setnum{4}} \minus {5 \over {256}}\varepsilon ^{\setnum{6}} \cr a_{\setnum{1}} \equals \tab \minus \left( {{3 \over 8}\varepsilon ^{\setnum{2}} \plus {3 \over {32}}\varepsilon ^{\setnum{4}} \plus {{45} \over {1024}}\varepsilon ^{\setnum{6}} } \right) \cr a_{\setnum{2}} \equals \tab \left( {{{15} \over {256}}\varepsilon ^{\setnum{4}} \plus {{45} \over {1024}}\varepsilon ^{\setnum{6}} } \right) \cr a_{\setnum{3}} \equals \tab \minus \left( {{{35} \over {3072}}\varepsilon ^{\setnum{6}} } \right) \cr}

Distance M 12 between two latitudes on the meridional arc can be determined using Equation 4 i.e.

M_{\setnum{12}} \equals f_{\setnum{1}} \lpar \psi _{\setnum{2}} \rpar \minus f_{\setnum{1}} \lpar \psi _{\setnum{1}} \rpar

Loss of significant digits is reduced for small angular separations if differencing is applied to Equation 4 resulting in:

(5)
M_{\setnum{12}} \equals a\left[ {a_{\setnum{0}} \lpar \psi _{\setnum{2}} \minus \psi _{\setnum{1}} \rpar \plus \sum\limits_{n \equals \setnum{1}}^{\setnum{3}} {2a_{n} \cos \lpar n\lpar \psi _{\setnum{2}} \plus \psi _{\setnum{1}} \rpar \rpar \sin \lpar n\lpar \psi _{\setnum{2}} \minus \psi _{\setnum{1}} \rpar \rpar } } \right]

which will be adapted later to give distance on the great ellipse. There is also a companion harmonic inversion series to Equation 4, described by Snyder and attributed to an earlier work (Adams, Reference Adams1921) that used the Lagrange Inversion Theorem to construct the inversion series. It provides geodetic latitude as a function of normalized meridional distance. The condensed form of this harmonic inversion series is:

(6)
f_{\setnum{2}} \lpar u\rpar \equals b_{\setnum{0}} u \plus \sum\limits_{n \equals \setnum{1}}^{\setnum{4}} {b_{n} \ \sin \lpar 2nu\rpar }

the constants for which are:

\eqalign{\tab b_{\setnum{0}} \equals 1 \cr \tab b_{\setnum{1}} \equals {3 \over 2}\varepsilon _{\setnum{1}} \minus {{27} \over {32}}\varepsilon _{\setnum{1}}^{\setnum{3}} \cr \tab b_{\setnum{2}} \equals {{21} \over {16}}\varepsilon _{\setnum{1}}^{\setnum{2}} \minus {{55} \over {332}}\varepsilon _{\setnum{1}}^{\setnum{4}} \cr \tab b_{\setnum{3}} \equals {{151} \over {96}}\varepsilon _{\setnum{1}}^{\setnum{3}} \cr \tab b_{\setnum{4}} \equals {{1097} \over {512}}\varepsilon _{\setnum{1}}^{\setnum{4}} \cr}

and

\varepsilon _{\setnum{1}} \equals \lpar 1 \minus \surd \alpha \rpar \sol \lpar 1 \plus \surd \alpha \rpar

For each value of the normalized distance u \equals {\pi \over 2}{M \over {M_{\setnum{0}} }}, the function f 2(u) returns a value of geodetic latitude ψ corresponding to the given meridional distance M. The constant M 0 is the meridional distance from the equator to the pole i.e. M_{\setnum{0}} \equals f_{\setnum{1}} \left( {{\pi \over 2}} \right) or, equivalently, M 0=a(a 0 π/2). Both of these series are periodic and can be used over arcs spanning any interval in the range 0<ψ<2π.

3.2. Numerical Results for the Meridional Ellipse

Table 1 shows a comparison of results for meridional arc distance acquired from Equation 3 and Equation 4 for angles L ψ1=180 ψ/π degrees in the range 0<ψ<π, together with the magnitude of their maximum difference.

Table 1.

A second evaluation involving Equations 4 (or 5) and 6 has been performed using angles L ψ1=180 ψ/π degrees for the range 0<ψ<π, which were then applied to Equation 5 to obtain distance. Results from Equation 5 in geodetic miles were first normalized as described then applied to Equation 6 to determine corresponding meridional angles. After converting to degrees, the results from the inversion provided by Equation 6 then closely replicated the initial angles as L ψ2 as shown in Table 2 where the magnitude of the maximum difference found among them is also shown in the bottom row.

Table 2.

From the results in Table 1, meridional distance computed from Equation 4 is seen to closely match those results from Equation 3 at a level of precision quite adequate for GE distance to be computed on the section ellipse. Latitude differences in Table 2 show that Equation 4 and thus Equation 5, together with inversion via the series of Equation 6, return the original values of latitude to within <7×10−9 degrees. After applying an adjustment to their constant terms, the two meridional equations, Equations 5 and 6, can be expected to provide inverse and forward solutions for the GE on the section ellipse with comparable accuracy as will be demonstrated.

4. APPLICATION TO THE GREAT ELLIPSE

A GE as illustrated in Figure 2 is the locus of all points at the intersection of the surface of a spheroid and an inclined plane that passes through its centre. From a node such as N1 where the GE crosses the equator, the path of the GE rises in latitude up to a turning point called the vertex (V). From the vertex it then turns back to the equator crossing it at the next node N2 180° from the first. It then proceeds towards the vertex in the opposite hemisphere, then on towards the starting node. The GE is a section ellipse and is seen to resemble an inclined version of the meridional ellipse but passes through the vertexes instead of the poles and has reduced ellipticity. Figure 2 shows a section of the GE in the northern hemisphere upon which P is a point defined by geocentric coordinates ϕ,θ. Geodetic latitude corresponding to ϕ is ψ and geocentric latitude of the vertex is ϕv.

Figure 2.

Figure 2 also shows the angle N1OP=φ′ that lies in the plane of the ellipse and is known as the geocentric great elliptic angle (GEA). Associated with the geocentric GEA is the geodetic GEA ψ′. At the vertex V, the line OV makes a right angle with the equatorial line N1 O N2.

4.1. Equation of the Great Ellipse

In Figure 3, P1 and P2 correspond to two position vectors X1 and X2 defined in terms of their geocentric latitude and longitude; each has the form:

{\bi X} \equals r \lpar \phi \rpar x

for which

(7)
x \equals \lsqb \hat{\bimiota }\ cos \ \phi \ cos \ \theta \plus \widehat {\bi j} \ cos \ \phi \ sin \ \theta \plus \widehat {\bi k}\ sin \ \phi \rsqb
r \lpar \phi \rpar \equals \sqrt {{\alpha \over {\lpar 1 \minus \varepsilon ^{\setnum{2}} \ cos \phi ^{\setnum{2}} \rpar }}}

and

\alpha \equals \lpar 1 \minus \varepsilon ^{\setnum{2}} \rpar

Geocentric latitudes and longitudes for P1 and P2 are φ1, θ1 and φ2, θ2 which are then inserted as appropriate into Eq.7 for each vector. Since the two vectors X1 and X2 are constants, then their cross product also contains constants, that is:

(8)
x_{\setnum{12}} \equals {{{\bi X}_{\setnum{1}} \times {\bi X}_{\setnum{2}} } \over {\left\vert {{\bi X}_{\setnum{1}} \times {\bi X}_{\setnum{2}} } \right\vert}}

which is perpendicular to the plane of the GE containing P1, P2 and the centre. The three constant components are:

(9)
\matrix{ {\lsqb \hat{\bimiota }m} \tab {\widehat {\bi j} n} \tab {\widehat{{\bi k}{p} }\rsqb \equals x_{\setnum{12}} } \cr}

An arbitrary vector X=X(φ, θ) in the plane of the section ellipse defined by any other position along the GE must satisfy the condition that X⋅x 12=0. The solution determined for this condition is the equation of the great ellipse (EGE), written here with \lambda \equals m \sol p and \mu \equals n \sol p as:

(10)
\tan \psi \equals \minus \lpar \lambda \ cos \theta \plus {\rm \mu} \ sin \theta \rpar \sol \alpha

The EGE relates latitude to longitude for all points along the ellipse and is required for Equation 2. The zeros of Equation 10 where φ=0, determine the longitudes of the two nodes as \theta _{n} \equals \tan ^{ \minus \setnum{1}} \left( {{\raise0.7ex\hbox{${\minus \lambda }$} \!\mathord{\left/ {\vphantom {{ \minus \lambda } \mu }}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$\mu $}}} \right) \pm \pi , while longitudes of the two vertexes occur at \theta _{v} \equals \tan ^{ \minus \setnum{1}} \left( {{\raise0.7ex\hbox{${ \mu }$} \!\mathord{\left/ {\vphantom {{ \minus \mu } \lambda }}\right.\kern-\nulldelimiterspace}\!\lower0.7ex\hbox{$\lambda $}}} \right) \pm \pi where the geocentric vertex latitudes are \phi _{v} \equals \pm {\rm arctan} \sqrt {\lambda ^{\setnum{2}} \plus \mu ^{\setnum{2}} } . The azimuth angle γ0 between the GE and the meridian on crossing the equator at a node is related to the vertex latitude φv by \gamma _{\setnum{0}} \equals {\pi \over 2} \minus \phi _{v} .

4.2. The Great Elliptic Angle (GEA)

For the section ellipse shown in Figure 3 to be a replica of the meridian ellipse, its own ellipticity must be replaced with a smaller value ε′ (0⩽ε′⩽ε) given by:

(11)
\varepsilon \prime \equals {{\varepsilon {\kern 1pt} {\rm sin}\phi _{v} } \over {\sqrt {\lpar 1 \minus \varepsilon ^{\setnum{2}} \cos \phi _{v}\hskip0pt^{\setnum{2}} \rpar } }}

A seen in Figure 3, the geocentric GEA φ′1 of position P1 is the angle between Xn and X1 defined by:

(12)
\cos \phi \prime_{\setnum{1}} \equals {{{\bi X}_{n} \cdot {\bi X}_{\rm \setnum{1}} } \over {\left\vert {{\bi X}_{n} } \right\vert\left\vert {{\bi X}_{\rm \setnum{1}} } \right\vert}}

Similarly for P2:

(13)
\cos \phi \prime_{\setnum{2}} \equals {{{\bi X}_{n} \cdot {\bi X}_{\rm \setnum{2}} } \over {\left\vert {{\bi X}_{n} } \right\vert\left\vert {{\bi X}_{\rm \setnum{2}} } \right\vert}}

The corresponding geodetic GEA for P1 is ψ′1 given by:

(14)
\psi \prime_{\setnum{1}} \equals {\rm arcos} \left[ {{{\alpha \prime{\kern 1pt} \cos \phi \prime_{\setnum{1}} } \over {\sqrt {\lpar 1 \minus \lpar 1 \minus \alpha \prime^{\setnum{2}} \rpar \cos \phi \prime_{\setnum{1}}\hskip0pt^{\setnum{2}} \rpar } }}} \right]

Similarly for P2

(15)
\psi \prime_{\setnum{2}} \equals {\rm arcos} \left[ {{{\alpha \prime{\kern 1pt} \cos \phi \prime_{\setnum{2}} } \over {\sqrt {\lpar 1 \minus \lpar 1 \minus \alpha \prime^{\setnum{2}} \rpar \cos \phi \prime_{\setnum{2}}\hskip0pt^{\setnum{2}} \rpar } }}} \right]

where α′=1−ε′.

For an arbitrary vector X=X(φ, θ) defined by any position P on the GE, the geocentric GEA is:

(16)
\cos \phi \prime \equals {{{\bi X}_{n} \cdot X} \over {\left\vert {{\bi X}_{n} } \right\vert\left\vert {\bi X} \right\vert}}

and the corresponding geodetic GEA is:

(17)
\psi \prime \equals {\rm arcos} \left[ {{{\alpha \prime{\kern 1pt} \cos \phi \prime} \over {\sqrt {\lpar 1 \minus \lpar 1 \minus \alpha \prime^{\setnum{2}} \rpar \cos \phi \prime^{\setnum{2}} \rpar } }}} \right]

Figure 3.

4.3. Great Elliptic Distance (GED)

The GED measured along the GE from the nearest node can be calculated after replacing the constants α, ε and ψ in Eq.3 with α′ ε′ and ψ′ i.e.

(18)
F \prime\lpar \psi \prime\rpar \equals a\alpha \prime\int\limits_{\setnum{0}}^{\psi \prime} {{{d\xi } \over {\sqrt {\lpar 1 \minus \varepsilon \prime ^{\setnum{2}} \sin \xi ^{\setnum{2}} \rpar ^{\setnum{3}} } }}}

The distance D′12 between two positions P1 at ψ′1 and P2 at ψ′2 then becomes:

(19)
D \prime_{\setnum{12}} \equals F \prime\lpar \psi \prime_{\setnum{2}} \rpar \minus F \prime\lpar \psi \prime_{\setnum{1}} \rpar

It is however a simpler matter to also make use of Equation 5 by suitably modifying the constant terms. This entails redefining the constants a n by replacing ε therein with ε′ so that the a n become a′n. With the constants so adjusted, the series for distance D′12 along the GE arc written as a function becomes:

(20)
D \ \lpar \psi \prime_{\setnum{1}}\comma \psi \prime_{\setnum{2}} \rpar \equals a\left[ {a \prime_{\setnum{0}} \lpar \psi \prime_{\setnum{2}} \minus \psi \prime_{\setnum{1}} \rpar \plus \sum\limits_{n \equals \setnum{1}}^{\setnum{3}} {2a \prime_{n} \cos \lpar n\lpar \psi \prime_{\setnum{2}} \plus \psi \prime_{\setnum{1}} \rpar \rpar \sin \lpar n\lpar \psi \prime_{\setnum{2}} \minus \psi \prime_{\setnum{1}} \rpar \rpar } } \right]

4.4. Numerical Evaluation (GED)

Numerical results obtained from Equation 20 for GED in geodetic miles have been compared with results from Equation 2 and are shown in Table 3. The paths chosen for the comparison begin and end at latitudes shown in column 1, for which the chosen span of longitude in each case is 1350. Column 2 contains the calculated results from the line integral of Equation 2, while column 3 shows results from Equation 20. In the worst case, the difference between the two schemes amounted to a maximum of 7.7×10−7 miles. From these results, it is apparent that Equation 20 is a useful alternative for calculating distance along a great ellipse.

Table 3.

5. GE FORWARD SOLUTION BY INVERSION FORMULA

The forward solution entails setting a distance to go from a given geographic position in a specified direction and although the EGE can be determined from a single position if the azimuth is also specified, by itself such a procedure is least likely to appeal to the navigator. If a series of waypoint distances are specified for a solution to the EGE established from two geographic end points, then the latitudes and longitudes of each successive waypoint can be determined from the forward solution set out as follows.

Application of the inversion formula of Equation 6 to the GE requires that the constants b n therein be replaced by b′ n which are obtained on replacing ε1 with ε′1 in the series for b n where now:

\varepsilon \prime_{\hskip-1pt \setnum{1}}\equals \lpar 1 \minus \surd {\alpha \prime} \rpar \sol \lpar 1 \plus \surd {\alpha \prime} \rpar

So then:

(21)
f \prime_{\setnum{2}} \lpar u\rpar \equals b \prime_{\setnum{0}} u \plus \sum\limits_{n \equals \setnum{1}}^{\setnum{4}} {b \prime_{n} } \sin \lpar 2nu\rpar

When distances between waypoints have been set, then the geocentric GEA to each waypoint from the node can be determined from Equation 21. But first, it is necessary to add the distance D ge1 between the nearest node and ψ′1 to the given distance D ge between the departure point at ψ′1 and the point or waypoint whose overall angular distance ψ′ is to be found.

So with

(22)
D_{\setnum{0}} \equals D\left( {0\comma {\pi \over 2}} \right)
(23)
D_{ge\setnum{1}} \equals D\lpar 0\comma \psi \prime_{\setnum{1}} \rpar

The required normalized GED is then:

(24)
D \equals {\pi \over 2}\left( {{{D_{ge\setnum{1}} \plus D_{ge} } \over {D_{\setnum{0}} }}} \right)

From Equation 21 the corresponding geodetic GEA is obtained as:

(25)
\psi \prime \equals f \prime_{\setnum{2}} \lpar D\rpar

and the corresponding geocentric GEA is:

(26)
\phi \prime \equals arctan \lpar \alpha \prime\tan \psi \prime\rpar

It is then a simple matter using spherical formulae, to recover the geocentric latitude ɸe corresponding to each GEA i.e.

(27)
\phi _{e} \equals arcsin \lpar \sin \phi \prime\sin \phi _{\nu } \rpar

followed by the geodetic latitude ψe from:

(28)
\psi _{e} \equals arctan \left( {{{\tan \phi _{e} } \over \alpha }} \right)

Forward azimuth γe on the GE at each waypoint is also recovered on using:

(29)
\gamma _{e} \equals arctan \lpar \tan \gamma _{s} \cos \lpar \psi _{e} \minus \phi _{e} \rpar \rpar

in which γs is the forward azimuth on the sphere given by:

(30)
\gamma _{s} \equals arcsin \left( {{{\cos \phi _{v} } \over {\cos \phi _{e} }}} \right)

When required, the reverse or back azimuth γr at each waypoint is readily determined as it is the supplement of the forward azimuth i.e. γr+γe=π. As remarked in the introduction section, the azimuth difference between the GC and the GE is small; analysis shows that if the term cos(ψe−φe) is ignored, then it can be demonstrated that the difference between these azimuths is no greater than {\varepsilon ^{\setnum{4}}\sol 16} or 1·6×10−4 degrees.

Intervals of longitude Δ, between successive waypoints and the longitude of the nearest node θn, are recovered upon using:

(31)
\rmDelta \equals \arcsin \left( {{{\sin \phi _{e} \sin \gamma _{s} } \over {\sin \phi _{v} }}} \right)

Or

(32)
\rmDelta \equals \arcsin \left( {{{\tan \phi _{e} } \over {\tan \phi _{v} }}} \right)

Adding these differences to the longitude of the node at θn gives the longitude of each waypoint i.e.

(33)
\theta \equals \theta _{n} \plus \rmDelta

Equations 27 through 32 are simple solutions to the spherical triangle located in the northern hemisphere. A variety of alternative expressions also exists (American Practical Navigator, 1981) and can be found in several other standard references as well.

In the following subsections, a baseline solution to the GE is first established that gives distance and azimuth over a wide span of longitude that contains a number of waypoints equi-spaced in longitude. Distance solutions obtained from the baseline are then applied to the equations of this section so as to recover as closely as possible the baseline data and thereby assess the accuracy and versatility of the inverse and forward solution methods described.

5.1. Baseline of Numerical Results

In this sub-section, a baseline solution for a GE track is established, against which results from a forward solution obtained using the expressions in this section have been compared. The baseline was established from the line integral solution using Equation 2 for distance plus Equation 17 of Earle, (Reference Earle2008) for azimuth (alternatively, Equation 35 therein could also have been used). A number of intermediate waypoints equi-spaced in longitude were placed between two widely separated end points, followed by the computation of their distances and azimuths. These results are shown in Table 4 and set out as follows. Column 1 is the range of longitude from a departure at 0°E up to a destination at 135°E with eight intermediate waypoints stepped in 15° intervals for a track beginning and ending at 10°N. Column 2 is the geodetic latitude on the GE computed from the EGE of Equation 10 converted to degrees. Column 3 is the distance to each point along the track in geodetic miles from Equation 2, while column 4 contains the forward azimuth at each waypoint in degrees.

Table 4.

5.2. The Forward Solution-Numerical Results

The same geographic endpoint pairs used for Table 4 were next used to re-establish the GE track on the section ellipse using the formulae of sub-sections 4.1 and 4.2. The geodetic GEAs corresponding to each waypoint position were also recalculated together with new values for distance from Equation 20. After using Equation 24 to normalize the results for distance from Equation 20, these normalized distance values were then inserted into the inversion formula Equation 21 via Equation 25, which allowed new geodetic GEAs to be recovered. New values of geocentric GEA were then recovered via Equation 26, after which, latitudes, azimuths and longitudes, corresponding to the baseline data set were determined anew from Equations 27 through 33.

In Table 5, the results from the forgoing steps are set out as follows. Longitude in column 1 is the value recovered using Equations 32 and 33 after using Equations 21 through 28 to obtain the latitudes shown in column 2. Column 3 shows the given distances derived from Equation 20, while azimuth values in column 4 are the values recovered from Equation 29. These results summarized in Table 5 are seen to closely replicate the baseline data of Table 4.

Table 5.

The final row of Table 5 shows maximum differences found in each column for the recovered data when compared with its baseline counterpart in Table 4. Other tests over several variants of the exercise showed similar results. For example, for baseline latitude beginning and ending at 80°N, the resulting maximum differences for columns 1, 2, 3, and 4 of Table 5 became 1·54×10−8, 3·7×10−11, 5·4×10−7 and 1·5×10−8 respectively.

6. CONCLUDING REMARKS

A closed form harmonic series approximation for determining distance between points on a great ellipse has been described and demonstrated. It is found to be a useful alternative to an established integral expression for distance and also found to give accurate results. A companion inversion harmonic series has also been demonstrated that converts elliptic distance into a geodetic GEA, from which the geodetic position coordinates for waypoints on the GE were constructed along with azimuth directions. These two series are seen to be simple, direct and computationally efficient. They provide accuracies for distance in the sub-metre range that are commensurate with the current levels of accuracy achieved by satellite positioning systems, while at the same time, they provide a complete solution to the Great Ellipse.

References

REFERENCES

Adams, . (1921). Latitude Developments Connected With Geodesy and Cartography. Special Publication 67, Department of Commerce, 123–125.Google Scholar
Bowditch, , (1981) American Practical Navigator, Volume II. Defense Mapping Agency, Pub No 9, 604–607.Google Scholar
Bowring, . (1984). The Direct and Inverse Solutions for the Great Elliptic Line on the Spheroid. Bulletin Geodesique 58, 101–108.CrossRefGoogle Scholar
Earle, . (2000). A Vector Solution for Navigation on the Great Ellipse. The Journal of Navigation, 53, 473–481.CrossRefGoogle Scholar
Earle, . (2006). Sphere to Spheroid Comparisons. The Journal of Navigation, 59, 491–496.CrossRefGoogle Scholar
Earle, . (2008). Vector Solutions for Azimuth. The Journal of Navigation, 61, 537–545.CrossRefGoogle Scholar
Lipshutz, . (1969) Differential Geometry. Outline Series, McGraw Hill, 171–173.Google Scholar
Martinez, de Osés. (2005). Graphics Theory to Optimize the Navigation. European Journal of Navigation. ISSN 1503-3953, 1–6.Google Scholar
Massey, and Kestelman, . (1958) Ancillary Mathematics. Sir Isaac Pitman & Sons, London, 874.Google Scholar
Pallikaris, and Latsas, . New Algorithm for Great Elliptic Sailing (GES). The Journal of Navigation, 62, 493–507.CrossRefGoogle Scholar
Snyder, . (1987). Map Projections-A Working Manual. Professional Paper No. 1395. U.S. Geological Survey.CrossRefGoogle Scholar
Williams, . (1998). Geometry of Navigation. Horwood Publishing. ISBN 1-898563-46-2.Google Scholar
Figure 0

Figure 1.

Figure 1

Table 1.

Figure 2

Table 2.

Figure 3

Figure 2.

Figure 4

Figure 3.

Figure 5

Table 3.

Figure 6

Table 4.

Figure 7

Table 5.