# Note on the direct computation of geodetic distances and azimuths on an ellipsoid of revolution.

1. INTRODUCTION

A classical problem of geodesy and many branches of geophysics is the computation of the distance and azimuths between two points on the reference earth ellipsoid when only the latitude and longitude of each point are known. This problem is referred to as the second main geodetic problem. The principal difficulty is that the general geodesic on the ellipsoid of revolution is a space curve, and hence the distance and azimuths cannot be expressed directly in terms of the geographical coordinates of the endpoints. The problem has been solved by applying iterative techniques (Helmert, 1880), or expansions with respect to the flattening.

In seismology we frequently encounter this problem in determining epicentral distances and the corresponding azimuths. The problem also acquires increasing importance in the controlled-source seismology, since many portable seismic instruments are now equipped with receivers of the GPS signals for determining their position. Consequently, the traditional applications of topographic maps for determining distances and azimuths are being replaced by computations from the geographical coordinates of the endpoints. Similar problems are also encountered in air navigation, etc. In these practical applications, the geodetic accuracy is not necessary, but a greater attention is rather paid to a simplification of the mathematical solution and a reduction of the computational time. Expansions with respect to the Earth's flattening seem to be sufficient and convenient solutions in these cases.

The method of expansions with respect to the flattening has its origins in the 19th century. Its development was accompanied with many discussions, disputes about the authorship, corrections of some errors and gradual improvements (Pick, 1998). Some expansions to first order terms in flattening are contained in the Andoyer-Lambert formulae (Andoyer, 1932; Lambert, 1942), but an exact theory of their derivation is not known. Thomas (1965) extended these formulae to second-order terms in flattening. Extensive tests of his formulae were performed by Pick (1998).

Thomas's formulae are convenient and sufficiently accurate in most geophysical problems. We have used them for many years in determining epicentral distances and azimuths in various studies of seismic surface waves (e.g., Proskuryakova et al., 1981). Recently, we have used them in the processing of data from the Central European Lithospheric Experiment Based on Refraction 2000 (CELEBRATION 2000 experiment); see Malek et al. (2001). However, we have found some specific situations when the accuracy in computing the azimuths was reduced. Consequently, we propose little changes in some formulae for the azimuths.

Consider an oblate ellipsoid of revolution. Denote its flattening f = (a - b)/a where a, b are the semi-major and semi-minor axes, respectively. Let A and B be two points on the ellipsoid. Denote [[phi]sub.A] and [[lambda].sub.A] the geodetic latitude and longitude of point A, and [[phi].sub.B], [[lambda].sub.B] the coordinates of point B. As opposed to Thomas (1965), east longitudes will be considered positive.

2. DISTANCES AND ANGLES ON A SPHERE

First, outline the expressions for a sphere, as they represent the zero-order terms in the expansions. For the moment, identify latitudes [[phi].sub.A], [[phi].sub.B] with geocentric latitudes. Consider spherical triangle [DELTA]ABC with vertex C at the North Pole (Fig. 1). We know the sides of the triangle a = 90[degrees] - [[phi].sub.B], b = 90[degrees] - [[phi].sup.A] and angle [gamma] = [DELTA][lambda] = [[lambda].sub.B] - [[lambda].sub.A].

Side c, in principle, can be determined from the cosine theorem,

cos c = sin [[phi].sub.A] sin [[phi].sub.B] + cos [[phi].sub.A] + cos [phi].sub.B] cos [DELTA][lambda]. (1)

[FIGURE 1 OMITTED]

However, if the angular distance c is close to zero or 180[degrees], the theorem leads to the computation of arccos x in the neighborhood of x = [+ or -] 1. Since the graph of this function is very steep there, a small error in its argument may lead to a large error in the result. To avoid this problem, write the first and last cosines in (1) as

cos c = 1 - 2 [sin.sup.2] c/2, cos [DELTA][lambda] = 1 - 2 [sin.sup.2] [DELTA][lambda]/2 - (2)

Using simple algebra and trigonometric identities we arrive at

[sin.sup.2] c/2 = [sin.sup.2] [[phi].sub.m] + H [sin.sup.2] [DELTA][lambda]/2 (3)

where

[[phi].sub.p] = 1/2 ([[phi].sub.A] + [[phi].sub.B]) [[phi].sub.m] = 1/2 ([[phi].sub.A] - [[phi].sub.B]), (4)

H = [cos.sup.2] [[phi].sub.p] - [sin.sup.2] [[phi].sub.m];

see (Thomas, 1965). Formula (3) already removes the above-mentioned disadvantages of the cosine theorem (1).

In computing angles [alpha] and [beta] we now use the following equations (Rektorys 1994):

cos [[[alpha] + [beta]]/2] cos c/2 = sin [[phi].sub.p] sin [DELTA][lambda]/2, (5)

sin [[[alpha] + [beta]]/2] cos c/2 = cos [[phi].sub.m] cos [DELTA][lambda]/2, (6)

cos [[[alpha] + [beta]]/2] sin c/2 = cos [[phi].sub.p] sin [DELTA][lambda]/2, (7)

sin [[[alpha] + [beta]]/2] sin c/2 = sin [[phi].sub.m] sin [DELTA][lambda]/2, (8)

[[phi].sub.p] and [[phi].sup.m] being given by (4). To determine ([alpha] + [beta])/2 we use that of formulae (5) and (6) whose right-hand side is smaller in the absolute value. The reason is again reason is again to avoid the computation of inverse trigonometric functions in the regions of their rapid changes. In a similar way we use formulae (7) or (8) in computing ([alpha] - [beta])/2.

By multiplying both sides of (5), (6), and both sides of (7), (8), we arrive at

sin ([alpha] + [beta]) = k/M sin ([DELTA][lambda]), (9)

sin ([alpha] - [beta]) = K/L sin ([DELTA]([lambda]), (10)

where

L = [sin.sup.2] c/2, M = 1 - L, (11)

k = sin [[phi].sub.p] cos [[phi].sub.m], K = cos [phi].sub.p] sin [[phi].sup.m].

Formulae (9) and (10) have been used by Thomas (1965). They seem to be simpler than formulae (5) to (8). However, it again follows from the graph of arcsine that, for example, the accuracy of formula (9) will be reduced if [alpha] + [beta] is close to 90[degrees]. Consequently, the more complicated formulae (5) to (8) should be preferred in writing a computer code.

3. A NUMERICAL EXAMPLE FOR A SPHERE

Let us construct a special simple case where the application of formulae (9) and (10) will not be convenient for computing the angles. Consider again a spherical triangle with vertex C at the North Pole (Fig. 1), but with vertices A and B at the same geocentric latitude of [[phi].sub.A] = [[phi].sub.B] 45[degrees]. Then it follows from (4), (3), (11) and (9) that H = 0.5, k [square root of 2]/2, K = 0,

L = [sin.sup.2] c/2 = 1/2 [sin.sup.2] = [DELTA][lambda]/2 = [[1 - cos [DELTA][lambda]]/4],

sin([alpha] + [beta]) = [2[square root of 2] sin [DELTA][lambda]/[3 + cos [DELTA][lambda]]

Now, let us specify a suitable value of the longitude difference [DELTA][lambda]. Choose cos [DELTA][[lambda] = -(1/3), i.e. approximately [DELTA][lambda] = 109.47[degrees]. Equations (9) and (10) then take the forms sin([alpha] + [beta]) = 1, sin([alpha] - [beta]) = 0. Although these equations yield the desired analytical solution, [alpha] = [beta] = 45[degrees], they indicate possible computational problems. Namely, formula (6) takes the form sin[([alpha] + [beta])/2] = [square root of 2]/2. We see that the solution of Eq. (9) requires the evaluation of arcsin1, whereas formula (6) operates with a smaller argument of arcsine, arcsin([square root of 2]/2). This should result in a higher accuracy of formula (6) for [DELTA][lambda] close to the value mentioned above.

We have verified the expected loss of precision of formulae (9) and (10) in the vicinity of [DELTA][lambda] = 109.47[degrees] by numerical tests. Using Fortran programs with single and double precision, we found that formulae (9) and (10) in single precision yielded errors of several tens of arc seconds. A part of the results are given in Table 1. These errors exceeded considerably the truncation errors of the abovementioned expansions for an ellipsoid, which, according to Thomas (1965), are only one second of arc. In applying (5) to (8) to the same problem, the accuracy was higher by about three orders. Moreover, the application of formula (9) in single precision even failed around [DELTA][lambda] = 109.47[degrees] since, as a consequence of rounding errors, the right-hand side of (9) became greater than one (see the missing values in Tab. 1). Consequently, we prefer formulae (5) to (8) to Thomas's formulae (9) and (10).

4. SPECIFICATION OF THE AZIMUTH

We shall measure the azimuth from the north in the clockwise direction. Denote [[alpha].sub.AB] the azimuth of B as seen from A, and analogously [[alpha].sub.AB] the azimuth of A from B. Thus,

[[alpha].sub.AB] = [alpha],[[alpha].sub.BA] = 360[degrees] - [beta] for sin [DELTA][lambda] [less than or equal to] 0, (12a)

which is the situation shown in Fig. 1, and

[[alpha].sub.AB] = 360[degrees] - [[alpha], [[alpha].sub.BA] = [beta] for sin [DELTA][lambda] < 0. (12b)

Note that Thomas (1965) has also measured the azimuth in the clockwise direction, but from the south.

5. CORRECTIONS FOR THE FLATTENING

We adopt the corrections for the flattening from Thomas (1965) without significant changes. In addition to (4) and (11), let us introduce the following notations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (13)

The geodetic arc length on the ellipsoid to second-order terms in flattening is then

S = a (c + f/4 W + [f.sup.2]/128 Z], (14)

c being the spherical angular distance given by (3).

The ellipsoidal corrections to azimuths, [delta][[alpha].sub.AB] and [[delta][alpha].sub.BA], satisfy the equations

[[delta][alpha].sub.BA] + [[delta][alpha].sub.AB] = -fH (T + 1) K/L sin([DELTA][lambda]), (15)

[[delta][alpha].sub.BA] - [[delta][alpha].sub.AB] = -fH (T - 1) k/M sin ([DELTA][lambda])

if our notation is used, i.e. [DELTA][lambda] = [[lambda].sub.B] - [[lambda].sub.A] and [[lambda] is measured positive toward the east. The azimuths on the ellipsoid then read

[([[alpha].sub.AB]).sub.el], = [[alpha].sub.AB] + [[delta][alpha].sub.AB], [([[alpha].sub.BA]).sub.el] = [[alpha].sub.BA] + [delta][[alpha].sub.BA], (16)

where [[alpha].sub.AB] and [[alpha].sub.BA] are the azimuths on a sphere; see (12).

Finally, let us give a numerical example taken from Thomas (1965). Consider the Clarke Ellipsoid 1866 with a = 6 378 206.4 m and f = 0.0033900753. Let the coordinates of points A and B on the ellipsoid be [[phi].sub.A] = 8[degrees] 58'25.0", [[lambda].sub.A] = -79[degrees]34'24.0" (Panama), and [[phi].sub.B] = 21[degrees]26'06.0", [[lambda].sub.B] = -158[degrees] 01'33.0" (Hawaii). Using double precision, the above-mentioned formulae have yielded the distance between the points S = 8466621.02m, and the azimuths [([[alpha].sub.AB]).sub.el] = 289.95470[degrees], [([[alpha].sub.BA]).sub.el] = 85.61963[degrees]. Note only that the appropriate longitudes in Thomas (1965) are positive, and the azimuths differ by 180[degrees] from the values just mentioned.

6. CONCLUSIONS

For geodetic distances we have adopted the method of Thomas (1965). The algorithm is described by formulae (4), (3), (11), (13) and (14). However, we have recognized some situations when Thomas's formulae for azimuths have a reduced accuracy. Consequently, for computing the azimuths we recommend a modified method based on formulae (4) to (8), (11) to (13), (15) and (16).

ACKNOWLEDGMENTS

We are indebted to Prof. M. Pick for his valuable comments and discussions. This work was supported by the Ministry of Education, Youth and Sports, Czech Republic, within the Research Project MSM113200004, and by the Grant Agency of the Czech Republic, Project 205/01/0481.

REFERENCES

Andoyer, H.: 1932, Formule donnant la longueur de la geodesique joignant 2 points de l'ellipsoide donnes par leurs coordonees geographiques. Bull. Geodesique, 34, 77-81.

Helmert, F.R.: 1880, Die mathematischen und physikalischen Theorieen der hoheren Geodasie. Teubner, Leipzig.

Lambert, W.D.: 1942, The distance between two widely separated points on the surface of the earth. J. Wash. Acad. Sci., 32, 125-130.

Malek, J., Broz, M., Fischer, T., Horalek, J., Hrubcova, P., Jansky, J., Novotny, O., Ruzek, B. and the CELEBRATION working group: 2001, Seismic measurements along short profiles in western Bohemia during the Celebration 2000 experiment. Acta Montana IRSM AS CR, ser. A, No. 18 (121), 15-28.

Pick, M.: 1998, On the exactness in geodesy. Voj. tech. informace, c. 58. VZU Praha, pp. 5-14, (in Czech).

Proskuryakova, T.A., Novotny, O. and Voronina, E.V.: 1981, Studies of the Earth's Structure by the Surface-wave Method (Central Europe). Nauka, Moscow, (in Russian).

Rektorys, K.: 1994, Survey of Applicable Mathematics. Kluwer Academic Publishers, Dordrecht.

Thomas, P.D.: 1965, Geodesic arc length on the reference ellipsoid to second-order terms in the flattening. J. Geophys. Res., 70, 3331-3340.

Oldrich NOVOTNY (1) and Jiri MALEK (2)

(1) Department of Geophysics, Faculty of Mathematics and Physics, Charles University, V Holesovickach 2, 180 00 Prague 8, Czech Republic, (Fax: +420-221912555; e-mail: oldrich.novotny@mff.cuni.cz)

(2) Institute of Rock Structure and Mechanics, Academy of Sciences of the Czech Republic, V Holesovickach 41, 182 09 Prague 8, Czech Republic, (Fax: +420-26886645; e-mail: malek@irsm.cas.cz)
```Table 1 The angles in a spherical triangle with vertex C at
the North Pole and vertices A and B at a geocentric latitude
of 45[degrees]. For a given longitude difference,
[DELTA][lambda], the angle at vertex A was computed by
several methods: from formula (9) with single precision
(value [[alpha].sub.1]); from formula (6) with single
precision ([[alpha].sub.2]); using (6) or (9) with double
precision ([alpha]). All angles are given in degrees. For
missing values, formula (9) failed.

[[alpha]                [[alpha]
[DELTA]                [[alpha].    .sub.1]    [[alpha]    .sub.2]
[lambda]    [alpha]      sub.1]    [alpha]      .sub.2]    [alpha]

109.40     45.037762   45.035664   -0.002099   45.037766   0.000003
109.42     45.027160   45.026169   -0.000991   45.027161   0.000001
109.44     45.016556   45.013988   -0.002567   45.016560   0.000004
109.46     45.005950         --          --    45.005955   0.000004
109.48     44.995344         --          --    44.995346   0.000002
109.50     44.984736   44.986012   0.001275    44.984737   0.000001
```