1. INTRODUCTION
In traditional navigation, the calculation of the elements of the shortest navigation path between two points on the surface of the Earth is usually conducted by the use of a spherical model of the Earth and the assumption that one minute of arc of any great circle is equal to one international nautical mile. It is well known that more accurate results can be obtained by the adoption of an ellipsoidal model of the Earth and the calculation of geodesic distances and azimuths.
The discrepancies between the results of shortest navigational path calculations on the spherical model of the Earth as great circle arcs and on the ellipsoidal model as geodesics, are in the order of 0·27% according to Tobler (Reference Tobler1964), and in the order of 0·5% according to Earle (Reference Earle2006). In reality these discrepancies can exceed 15 nautical miles (about 28·5 km). An example of such a discrepancy is shown through the calculation of the shortest navigational distance from a departure location on the east coast of Australia, such as the entrance to Sydney Harbour (ϕ: 33° 46·21′ S, λ: 151° 31·964′ E) to a destination point on the west coast of South America such as the approaches to Valparaiso in Chile (ϕ: 32° 59·998′ S, λ: 71° 36·675′ W). The calculation of this distance on the spherical Earth model with the assumption that one minute of a great circle arc is equal to the international nautical mile (1852 metres) yields a distance of 6113 nautical milesFootnote 1. The calculation of same distance on the WGS-84 ellipsoid, with very accurate geodetic methods of sub-metre accuracy as the method of Vicenty (1975), yields 6128·4 nautical miles. For this example the difference in calculated distances on the spherical model from those on the ellipsoid is more than 15 nautical miles (~28·5 km).
In practice very accurate results can be obtained by the execution of the calculations on the great ellipse rather than on the geodesic. The great ellipse is the line of intersection of the surface of the ellipsoid with the plain passing through its geometric centre O and the departure and destination points P1 and P2 (Figure 1). For surface navigation applications the great elliptic arc P1P2 approximates very closely the geodesic line and even for the longest possible navigational paths the discrepancies between geodesic and great elliptic arc are practically negligible (see Section 5 of this paper presenting the results of numerical tests and comparisons). Seeking the higher accuracies of the geodesic for sailing calculations does not have any practical value for marine navigation and simply adds more complexity to the calculations. Great elliptic calculations are simpler than those on the geodesic and much more accurate than the traditional navigation calculations on the great circle with a spherical earth model.
In the span of the last 25 years many interesting methods and formulas for great elliptic sailing computations have been proposed from the direct and inverse solutions of Bowring (Reference Bowring1984) up to the vector solutions for azimuths of Earle (Reference Earle2008). Most of these methods have offered valuable contributions for the complete, straightforward and accurate solution of the great elliptic sailing problem, without the need to use advanced numerical methods or commercial mathematical software.
Our proposed analytical method and algorithm solves the complete problem including not only the great elliptic arc distance, as in other proposed methods, but also other elements of the sailing such as the course at the departure and destination points and the geodetic coordinates of the vertex and the intermediate points along the great elliptic arc.
The second section of this paper presents the great elliptic sailing (GES) problem and the basic parameters of the great ellipse that are used for the subsequent calculations. The third section presents an overview of various methods and formulas which have been proposed in the span of the last 25 years, from the direct and inverse solutions of Bowring (Reference Bowring1984) up to the vector solutions for azimuths of Earle (Reference Earle2008). The fourth section presents the proposed equations and the algorithm for the complete solution of the great elliptic arc sailing, including the coordinates of intermediate points along the great elliptic arc. The fifth section presents the results of numerical tests for the evaluation of the proposed algorithm. The sixth section concludes the paper.
2. THE GREAT ELLIPTIC SAILING PROBLEM AND THE BASIC PARAMETERS OF THE GREAT ELLIPSE
As in other sailing methods, in Great Elliptic Sailing (GES) the calculations are distinguished in the direct and inverse problem. In the direct problem the known parameters are the geodetic coordinates of the departure point, the distance and azimuth of the destination point, the semi major axis and the eccentricity of the ellipsoid; the computed parameters are the geodetic coordinates of the destination point. In the inverse problem the known parameters are the geodetic coordinates of the departure and destination points, the semi major axis and the eccentricity of the ellipsoid; the computed parameters are the great elliptic arc length (great elliptic sailing distance) and the course at the departure and destination point.
The navigator is interested only in the solution of the inverse sailing problem. Nevertheless the solution of the direct problem can be exploited for the derivation of calculation methods of the coordinates of intermediate points of the great elliptic arc. (see Section 4.4). These intermediate points are required in order to approximate the GES by a series of successive rhumb line sailings.
2.1. Calculation of the length of the great elliptic arc
The calculation of the length of the great elliptic arc can be conducted by the use of standard geodetic Formula (1) for the length of the arc of the meridian M 0φ (Figure 2) after the proper modification of the parameters of the meridian ellipse with those of the great ellipse. This is better understood if we consider the great ellipse as an inclined version of the meridian ellipse (Figure 3). More specifically the two basic parameters that are required for the calculation of the length of the great elliptic arc are:
• The eccentricity e ge of the great ellipse. This is equivalent to the eccentricity e of the reference ellipsoid used in Formula (1).
• The geodetic great elliptic angle ϕge. This is equivalent to the geodetic latitude ϕ used in Formula (1).
Another basic element of the great ellipse is the geocentric great elliptic angle θge. This parameter is a prerequisite for the calculation of the geodetic great elliptic angle ϕge. Specific formulas for the calculation of these parameters of the great ellipse are presented in Section 4 of this paper.
Adopting the parameters ege and ϕ′ of the great ellipse into the meridian arc Formula (1), we obtain Formula (2) that provides the distance along the great elliptic arc. Many methods have been used for the computation of the integral of Formula (1). All these methods and formulas can be used for the calculation of the distance along the great elliptic arc by Formula (2).
Equations (1) and (2) can be transformed to an elliptic integral of the second type, which cannot be evaluated in a closed form. The calculation can be performed either by numerical integration methods, such as Simpson's rule, or by the binomial expansion of the denominator to rapidly converging series, retention of a few terms of these series and further integration by parts. This process yields results like Formula (3).
Equation (3) is the standard geodetic formula for the accurate calculation of the meridian arc length, which is proposed in a number of textbooks such as in Torge's Geodesy using up to sin(2ϕ) terms.
According to Snyder (1987) and Torge (Reference Torge2001), Simpson's numerical integration of Formula (1) does not provide satisfactory results and consequently the standard computation methods for the length of the meridian arc are based on the use of series expansion formulas, such as Formula (3).
2.2. Calculation of the forward and inverse azimuth (course at the departure and destination point)
Most of the methods and formulas that have been proposed for the calculation of the forward and inverse azimuth of the great elliptic arc require advanced numerical methods and thus they suggest either the use of commercial mathematical software (Earle Reference Earle2000, Reference Earle2008), or the use of the simple formulas for the calculation on the sphere (Walwyn Reference Walwyn1999). A straightforward method for the calculation of the forward and inverse azimuth, that can be easily used in the compilation of algorithms, without recourse to advanced numerical methods, has been proposed by Bowring (Reference Bowring1984). The methodology proposed by Bowring consists of two steps. Initially the forward and backward azimuths are calculated on the auxiliary unit sphere as in the classical spherical Earth model of the traditional navigation. The second and final step is the reduction of the calculated spherical azimuths to their ellipsoidal values for the great elliptic arc sailing.
3. EXISTING METHODS FOR THE SOLUTION OF THE GREAT ELLIPTIC SAILING
3.1. Bowring method for the direct and inverse solutions for the great elliptic line
Bowring (Reference Bowring1984) provides formulas for the solution of the direct and inverse great elliptic sailing problem. Bowring's formulas can be used for the calculation of the great elliptic arc length and the forward and backward azimuths, but no solution is proposed for the calculation of the coordinates of the vertex and the intermediate points.
The first set of formulas proposed by Bowring concerns the computation of some prerequisite parameters, such as the spherical azimuth, the minor eccentricity of the great ellipse and the parametric latitudes of the departure and destination points. The proposed parameters are subsequently used for the solution of the direct and inverse problem.
The method of Bowring for the calculation of the great elliptic arc length employs the use of an auxiliary geodetic sphere and various types of coordinates, such as, geodetic, geocentric, Cartesian and polar. These formulas for the great elliptic distance have been tested and it was proved that they provide very satisfactory results in terms of obtained accuracy. Nevertheless other simpler computation methods of the length of the great elliptic arc can be used by the employment of standard geodetic formulas for the length of the arc of the meridian, after the proper modification of the parameters of the meridian ellipse with those of the great ellipse, such as Formulas (1) and (2). The formulas used by Bowring for the calculation of the forward and backward azimuths, unlike those for the distance, are very much simpler than other methods of the same accuracy.
3.2. William's method for the computation of the distance along the great elliptic arc
Williams (Reference Williams1996) provides formulas for the computation of the sailing distance along the arc of the great ellipse. These formulas have the general form of the integral of Formula (2). For the computation of the eccentricity e ge and the geodetic great elliptic angle ϕge of Formula (2), Williams provides simple and compact formulas. For the evaluation of this integral Williams employs the cubic spline integration method of Plythian and Williams (Reference Plythian and Williams1985). In his paper Williams (Reference Williams1996) presents the results of comparative calculations for eight test lines that have been permanently used by Haiwara (Reference Hairawa1987) in comparison tests of sailing calculation methods. The method of Williams does not provide formulas for the computation of azimuths and coordinates of intermediate points along the great elliptic arc.
3.3. Earle's method for Vector Solutions (2000 and 2008)
Earle (Reference Earle2000) has proposed a method of computing distance along a great ellipse that allows the integral for distance to be computed directly using the built-in capabilities of commercial mathematical software. This obviates the need to write code in arcane computer languages. According to Earle, his method has been prepared with the syntax of a particular commercial mathematics package in mind. This work of Earle (Reference Earle2000) has been recently updated for the calculation of azimuth (Earle Reference Earle2008). The first paper of Earle (Reference Earle2000) also presents the results of comparative calculations for the great elliptic arc distance using the same eight test lines used by Williams (Reference Williams1996) and Haiwara (Reference Hairawa1987). Earle does not provide numerical results for the computation of azimuth (course) and the coordinates of the intermediate points along the great elliptic angle.
3.4. Walwyn's great ellipse algorithm
Walwyn (Reference Walwyn1999) presented an algorithm for the computation of the arc length along the great ellipse and the initial heading to steer. The algorithm uses various formulas for the calculation of distance and azimuths (courses). In some cases, probably for the sake of simplicity, these formulas are not the right ones used in standard geodetic computations, as the formulas for the transformation of the geodetic latitudes to geocentric. Walwyn does not provide numerical results for the computation of azimuth (course) and the coordinates of the intermediate points along the great elliptic angle.
4. THE PROPOSED NEW ALGORITHM FOR THE GREAT ELLIPTIC SAILING
Our proposed algorithm was initially developed as a supporting tool in another research work of the first author on the implementation of sailing calculations in GIS navigational systems (ECDIS and ECS). The complete great elliptic sailing problem is solved including, in addition to the great elliptic arc distance, the geodetic coordinates of an unlimited number of intermediate points along the great elliptic arc. The algorithm has been developed having a mind to avoid the use of advanced numerical methods, in order to allow for the convenient implementation even in programmable pocket calculators. Numerical tests and results showed that the accuracy achieved is a little higher than other methods (see Section 5).
The algorithm starts with the calculation of the eccentricity of the great ellipse and the geocentric and geodetic great elliptic angles (Figure 3) of the points of departure and destination. For this part of the algorithm we used the formulas proposed by Williams (Reference Williams1996) because they are simple, straightforward and provide accurate results. For the calculation of the length of the great elliptic arc we used the standard geodetic series expansion formulas for the meridian arc length that are presented in basic geodesy textbooks like Torge (Reference Torge2001) after their proper modification for the great ellipse.
For the calculation of the initial and final course we adopted the methodology proposed by Bowring (Reference Bowring1984). Initially the forward and backward azimuths are calculated on the auxiliary unit sphere as in the classical spherical earth model of the traditional navigation. Then the calculated spherical azimuths are reduced to their ellipsoidal values for the great elliptic arc. These parameters are required for the subsequent calculation of the geodetic coordinates of the intermediate points along the great elliptic arc.
The calculation of the geodetic coordinates of the intermediate points along the great elliptic arc is conducted by successive solutions of the direct great elliptic problem using the formulas proposed by Bowring (Reference Bowring1984). In these successive calculations, in order to avoid propagation of errors, the initial point is always the point of departure and the destination point is the intermediate point concerned. The known parameters in these direct problem solutions are: the geodetic coordinates of the point of departure, the calculated initial course at this point and the distance of the intermediate point from the point of departure. The method can be implemented for the calculation of the coordinates of the intermediate points either at a given successive distance along the great elliptic arc, or for any desired number of waypoints at equally spaced distances. Examples of calculations are presented in Section 5 of this paper. The complete set of the proposed new algorithm is:
4.1. Part I: Calculation of basic parameters on the great ellipse
Geocentric latitudes θ1, θ2:
Where:
ϕ1 and ϕ2 are the geodetic latitudes of the departure and destination points P1 and P2.
θ1 and θ2 are the corresponding geocentric latitudes.
Transformation of geodetic to Cartesian coordinates:
Where:
Parameters λw, μw that define the equation z=λwx+μwy of the plain of the great ellipse containing points O, P1 and P2:
Geocentric coordinates of the vertex V:
Eccentricity e ge of the great ellipse:
Longitude of the point E where the great ellipse crosses the equator (see Figure 3):
Geocentric great elliptic angles θ1, θ2:
Geodetic great elliptic angles ϕge1, ϕge2:
4.2. Part II: Calculation of the great elliptic distance
Length of the great elliptic arc:
4.3. Part III: Calculation of the initial and final course (forward and backward azimuths)
Spherical azimuths:
Where Δλ=λ2−λ1
Reduction of spherical azimuths to elliptic:
4.4. Part IV: Calculation of the coordinates of intermediate points
Number of intermediate points:
New Cartesian coordinates of the point of departure:
Auxiliary parameters of the great ellipse:
Parametric latitude of the point of departure:
Auxiliary angle Θ:
Where:
Distance of the point of departure from the major axis of the great ellipse:
For the nth intermediate point (n=1, 2, …, N):
Distance from the point of departure:
Distance from major axis and auxiliary angle Θ:
Parametric latitude of the nth intermediate point:
Where:
Angle between great ellipse's semi-major axis and X axis:
Cartesian coordinates of the nth intermediate point:
Geodetic coordinates (ϕn, λn) of the nth intermediate point:
5. NUMERICAL TESTS AND COMPARISONS
For our numerical tests and comparisons we used as a basis older numerical tests and comparisons conducted by Williams (Reference Williams1996) and Earle (Reference Earle2000). These tests and comparisons use eight specific routes (lines) that were initially used by Haiwara (Reference Hairawa1987) for the comparison of great circle and geodesic sailing calculations. These eight lines lie on the successive parallels 10°, 20° … 80° and all of them have the same difference of longitude equal to 100°. For better visual perception, these eight routes are shown over the geographical area of the Pacific Ocean in Figure 4.
The numerical tests and comparisons of Williams (Reference Williams1996) and Earle (Reference Earle2000) are restricted to the calculations of the sailing distance only and do not include other computed elements of the great elliptic sailing, such as the coordinates of intermediate points. The assessment of the accuracy in the calculated elliptic arc distance in these older numerical tests and comparisons is conducted by their comparison to the results of the Lambert correction method for the computation of the corresponding distances on the geodesics.
The weakness of these older tests employing the Haiwara dataset (Figure 3) is not only the lack of computations for the coordinates of intermediate points but also that they have not employed the most accurate methods for the calculation of the length of geodesics as a basis for their comparison tests, since they use the old Lambert correction method and not newer more accurate methods and algorithms such as those of Sodano (Reference Sodano1965) and Vicenty (Reference Vicenty1975). This weakness was also reported by Zukas (Reference Zukas1994).
For the above mentioned reasons, in the first stage of our numerical tests we recalculated the lengths of the Haiwara's eight lines on the Bessel ellipsoid using Vicenty's algorithm for the precise calculation of the length of the corresponding geodesic with sub-metre accuracy and compared these results with those obtained by the Lambert correction method for the length of the geodesic, as well as Williams' method, Earle's method and our proposed algorithm for the length of the great elliptic arc (see Figure 5 and Table 1).
Sv, Se, Sw, Sl, Sp are computed distances. Δe, Δw, Δl, Δp are errors computed as discrepancies from the results of Vicenty's algorithm.
The results of these comparison tests showed that:
• The average discrepancy between the calculated geodesic length by the Lambert Correction Method and Vicenty's algorithm is 12·6 metres, which is twice the discrepancy of the results of Williams and Earle. In other words the calculation of the sailing distance on the great elliptic arc by Williams and Earle is more accurate than the calculation of corresponding sailing distance along the geodesic conduced by the Lambert correction method. In view of this finding the comparison tests have to be re-evaluated on the basis of their comparison to Vicenty's algorithm rather than to the Lambert correction method.
• The average error of our proposed method in the calculation of the great elliptic arc distances is 4·38 metres which is smaller than the 6·54 metres average errors of the methods of Williams and Earle.
The maximum error of our proposed method in the calculation of the great elliptic arc distances is 10·85 metres which is smaller than the 15·72 metres maximum error of the methods of Williams and Earle. This maximum error occurs in the line on the 10° parallel. This line has a length of 5526·95 geographical miles.
The second stage of our numerical tests was the solution (in WGS-84) of the complete great elliptic sailing (GES) problem including the computation of the coordinates of intermediate points along the great ellipse for two very long navigational routes.
The first numerical example of very long navigational route with difference of longitude greater than 130° is the sailing from the approaches of Sydney Harbour-Australia (33° 46·21′ S, 151° 31·964′ E) to the approaches of Valparaiso-Chile (32° 59·998′ S, 71° 36·675′ W). The results of these calculations are shown in Table 2. In this calculation we selected a distance of 50 nautical miles between successive intermediate points along the great ellipse (the distance between intermediate points is selected by the user and can be as short as desired). Comparing the calculated with our algorithm value of great elliptic sailing distance between Sydney and Valparaiso (6129·12 nautical miles) with the corresponding geodesic distance calculated by Vicenty's algorithm (6128·41 nautical miles) we see that even for this extremely long distance with difference of longitude of some 137°, the bigger discrepancy (0·71 nautical miles) is still negligible for the practical purposes of marine navigation.
The second numerical example of very long navigational route with difference of longitude greater than 145° is the sailing from Valparaiso-Chile (32° 59·998′ S, 71° 36·675′ W) to Hokohama-Japan (34° 26·178′ N, 139° 51·39′ E). In this example the discrepancy between the great elliptic distance computed with our algorithm (9242·80 nautical miles) and the geodesic distance computed by Vicenty's algorithm (9241·92 nautical miles) is a little bigger (0·88 nautical miles) but still acceptable for the practical purposes of marine navigation.
It is noted that even for these two extreme cases where the difference of longitude between departure and destination points is greater than 130° and 145° the resulting discrepancies, that are still less than one nautical mile, are practically diminished in the process of the computation of the coordinates of the intermediate points. Our algorithm computes very easily these coordinates for as many intermediate points as desired (see Table 2).
6. CONCLUSIONS
Numerical tests and comparisons showed that great elliptic sailing distances are practically the same with the equivalent geodesic distances. Even for extremely long navigational paths (≈9·000 nautical miles) with difference in longitude greater than 135° the bigger discrepancy in the computed distances between the geodesic and the great elliptic arc is still negligible for the practical purposes of marine navigation.
The extremely high accuracy of the proposed new method and algorithm for the complete solution of the great elliptic sailing problem has been verified with Vicenty's algorithm which is one of the most accurate geodesic methods for the computation of long geodesics with sub-metre accuracy. Comparison tests showed that the proposed algorithm compared to the methods of Williams and Earle for the computation of the great elliptic distance, yields slightly better accuracies.
The main advantages of the proposed algorithm for the great elliptic sailing calculations, is not the above mentioned slightly better accuracy but the capacity to solve the complete sailing problem and not only the great elliptic arc distance as in other proposed methods. Other elements of the great elliptic sailing are also calculated, such as the geodetic coordinates of the vertex and those of an unlimited number of intermediate points along the great elliptic arc. The proposed algorithm does not require the use of advanced numerical integration methods as is the case in other methods that have been proposed in the past. The simplicity of the method also permits its implementation on programmable pocket calculators for the complete solution of the great elliptic arc sailing.