GEODESICS ON AN ELLIPSOID – PITTMAN'S METHOD

R.E. Deakin and M.N. Hunter

School of Mathematical and Geospatial Sciences, RMIT University
GPO Box 2476V, MELBOURNE, VIC 3001
Email: [email protected]
Presented at the Spatial Sciences Institute Biennial International Conference (SSC2007),
Hobart, Tasmania, 14-18 May 2007.

ABSTRACT
The direct and inverse problems of the geodesic on an ellipsoid are fundamental
geodetic operations. This paper presents a detailed derivation of a set of recurrence
relationships that can be used to obtain solutions to the direct and inverse problems with
sub-millimetre accuracies for any length of line anywhere on an ellipsoid. These
recurrence relationships were first described by Pittman (1986), but since then, little or
nothing about them has appeared in the geodetic literature. This is unusual for such an
elegant technique and it is hoped that this paper can redress this situation. Pittman's
method has much to recommend it.

BIOGRAPHIES OF PRESENTERS
Rod Deakin and Max Hunter are lecturers in the School of Mathematical and Geospatial
Sciences, RMIT University; Rod is a surveyor and Max is a mathematician, and both
have extensive experience teaching undergraduate students.

INTRODUCTION
Twenty-one years ago (March 1986), Michael E. Pittman, an assistant professor
of mathematical physics with the Department of Physics, University of New Orleans,
Louisiana USA, published a paper titled 'Precision Direct and Inverse Solutions of the
Geodesic' in Surveying and Mapping (the journal of the American Congress on
Surveying & Mapping, now called Surveying and Land Information Systems). It was
probably an unusual event – a physicist writing a technical article on geodetic
computation – but even more unusual was Pittman's method; or as he put it in his paper,
"The following method is rather different." And it certainly is.
Usual approaches could be roughly divided into two groups: (i) numerical
integration schemes and (ii) series expansion of elliptic integrals. The first group could
be further divided into integration schemes based on simple differential relationships of
the ellipsoid (e.g., Kivioja 1971, Jank & Kivioja 1980, Thomas & Featherstone 2005),
or numerical integration of elliptic integrals that are usually functions of elements of the
ellipsoid and an auxiliary sphere (e.g., Saito 1970, 1979 and Sjöberg 2006). The second
group includes the original method of F. W. Bessel (1826) that used an auxiliary sphere
and various modifications to his method (e.g., Rainsford 1955, Vincenty 1975, 1976 and
Bowring 1983, 1984).

1
 Pittman developed simple recurrence relationships for the evaluation of elliptic
integrals that yield distance and longitude difference between a point on a geodesic and
the geodesic vertex. These equations can then be used to solve the direct and inverse
problems. Pittman's technique is not limited by distance, does not involve any auxiliary
surfaces, does not use arbitrarily truncated series and its accuracy is limited only by
capacity of the computer used.
Pittman's paper was eight pages long and five of those contained a FORTRAN
computer program. In the remaining three pages he presented a very concise
development of two recurrence relationships and how they can be used to solve the
direct and inverse problems of the geodesic on an ellipsoid (more about this later). His
paper, a masterpiece of brevity, contained a single reference and an acknowledgement
to Clifford J. Mugnier – then a lecturer in the Department of Civil Engineering,
University of New Orleans – for numerous discussions. Unlike other published
methods which have been discussed and developed in detail over the years, Pittman's
method seems to have received no further treatment to our knowledge in the academic
literature, excepting brief mentions in bibliographies and reference lists. Our purpose,
in this paper, is to explain Pittman's elegant method as well as provide some useful
information about the properties of the geodesic on an ellipsoid.

The Direct and Inverse problems of the geodesic on an ellipsoid

In geodesy, the geodesic is a unique curve on the surface of an ellipsoid defining
the shortest distance between two points. A geodesic will cut meridians of an ellipsoid
at angles α , known as azimuths and measured clockwise from north 0º to 360 . Figure
1 shows a geodesic curve C between two points A (φ A , λA ) and B (φB , λB ) on an
ellipsoid. φ , λ are geodetic latitude and longitude respectively and an ellipsoid is taken
to mean a surface of revolution created by rotating an ellipse about its minor axis, NS.

n vertex φmax
dia

C
e ri

• B
•α
hm

αAB de s ic BA
geo
Greenwic

al

A
• b
rm

no
rm
no

a
φA l O φB
•
a HA • H a
B
equator
λB
λA

Fig. 1: Geodesic curve on an ellipsoid

2
 The geodesic curve C of length s from A to B has a forward azimuth α AB
measured at A and a reverse azimuth α BA measured at B and α AB ≠ α BA . The direct
problem on an ellipsoid is: given latitude and longitude of A and azimuth α AB and
geodesic distance s, compute the latitude and longitude of B and the reverse azimuth
α BA . The inverse problem is: given the latitudes and longitudes of A and B, compute
the forward and reverse azimuths α AB , α BA and the geodesic distance s.
The geodesic is one of several curves of interest in geodesy. Other curves are: (i)
normal section curves that are plane curves containing the normal at one of the terminal
points; in Figure 1 there would be two normal section curves joining A and B and both
would be of different lengths and also, both longer than the geodesic; (ii) curve of
alignment that is the locus of all points Pk where the normal section plane through Pk
contains the terminal points of the line; and (iii) great elliptic arcs that are plane curves
containing the terminal points of the line and the centre of the ellipsoid. Normal section
curves, curves of alignment and great elliptic arcs are all longer than the geodesic and
Bowring (1972) gives equations for the differences in length between these curves and
the geodesic.

Some ellipsoid relationships

The size and shape of an ellipsoid is defined by one of three pairs of parameters:
(i) a, b where a and b are the semi-major and semi-minor axes lengths of an ellipsoid
respectively, or (ii) a, f where f is the flattening of an ellipsoid, or (iii) a, e 2 where e 2
is the square of the first eccentricity of an ellipsoid. The ellipsoid parameters a, b, f , e 2
are related by the following equations
a −b b a 2 − b2 b2
f = = 1 − ; b = a (1 − f ) ; e 2 = = 1 − = f (2 − f ) (1)
a a a2 a2
The second eccentricity e′ of an ellipsoid is also of use and
a 2 − b2 e2 f (2 − f )
( e′) = = =
2
(2)
1− e (1 − f )
2 2 2
b
In Figure 1, the normals to the surface at A and B intersect the rotational axis of
the ellipsoid (NS line) at H A and H B making angles φ A , φB with the equatorial plane of
the ellipsoid. These are the latitudes of A and B respectively. The longitudes λA , λB are
the angles between the Greenwich meridian plane and the meridian planes ONAH A and
ONBH B containing the normals through A and B. φ and λ are curvilinear coordinates
and meridians of longitude (curves of constant λ ) and parallels of latitude (curves of
constant φ ) are parametric curves on the ellipsoidal surface. Planes containing the
normal to the ellipsoid intersect the surface creating elliptical sections known as normal
sections. Amongst the infinite number of possible normal sections at a point, each
having a certain radius of curvature, two are of interest: (i) the meridian section,
containing the axis of revolution of the ellipsoid and having the least radius of
curvature, denoted by ρ (rho), and (ii) the prime vertical section, perpendicular to the
meridian plane and having the greatest radius of curvature, denoted by ν (nu).

3
 a (1 − e 2 ) a
ρ= and ν= (3)
(1 − e sin φ ) (1 − e sin φ )
3 1
2 2 2 2 2 2

In the development that follows, use

will be made of relationships that can be α+dα
obtained from the differential rectangle on Q
• φ+dφ
the ellipsoid shown in Figure 2. Here P

ρ dφ
α
and Q are two points on the surface ds
P
connected by a curve of length ds with •
azimuth α at P. The meridians λ and ν cos φ dλ φ
λ + d λ , and parallels φ and φ + dφ form
λ
a differential rectangle on the surface of λ+dλ
the ellipsoid.
From Figure 2 the following Fig. 2: Differential rectangle on ellipsoid
relationships can be obtained
ds sin α = ν cos φ d λ and ds cos α = ρ dφ (4)

Mathematical definition of a geodesic

A curve drawn on a surface so that its osculating plane at any point on the surface
contains the normal to the surface is a geodesic (Lauf 1983). This definition, including
a definition of the osculating plane, can be explained briefly by the following.
A point P on a curve (on a surface) has a position vector
r ( t ) = x ( t ) i + y ( t ) j + z ( t ) k where i,j,k are unit vectors in the directions of the x,y,z
Cartesian coordinate axes and t is some scalar parameter. As t varies then the vector r
sweeps out the curve C on the surface, hence the distance s along the curve is a function
ds d
of t, given via = r ( t ) . Differentiating the vector r with respect to s gives a unit
dt dt
tangent vector t and differentiating t with respect to s gives the curvature vector κ n ,
perpendicular to t. n is the principal normal vector, κ (kappa) is the curvature and
1
ρ = is the radius of curvature and also the radius of the osculating (kissing) circle
κ
touching P.
The osculating plane at P contains both t and n (and the osculating circle), and
when this plane also contains the normal to the surface then the curvature κ is least and
ρ is a maximum; this is Meunier's theorem (Lauf 1983), a fundamental theorem of
surfaces. Therefore, if P and Q are very close and both lie on the surface and in the
osculating plane, then the distance ds between them is the shortest possible distance on
the suface.

4
The characteristic equation of a geodesic
The mathematical definition of a geodesic does little to help us develop solutions
to the problem of computing distances of geodesics on an ellipsoid. It does lead to the
characteristic equation of a geodesic, and this equation is the basis of all solutions to
computing geodesic distances. This equation
ν cos φ sin α = constant (5)
is known as Clairaut's equation in honour of the French mathematical physicist Alexis-
Claude Clairaut (1713-1765). In a paper in 1733 titled Determination géométric de la
perpendicular à la méridienne tracée par M. Cassini, ... Clairaut made an elegant study
of the geodesics of surfaces of revolution and stated his theorem embodied in the
equation above (Struik 1933). His paper also included the property already pointed out
by Johann Bernoulli (1667-1748): the osculating plane of the geodesic is normal to the
surface (DSB 1971)
The characteristic equation of a geodesic shows that the geodesic on the ellipsoid
has the intrinsic property that at any point, the product of the radius r = ν cos φ of the
parallel of latitude and the sine of the azimuth, sin α , of the geodesic at that point is a
constant. This means that as r decreases in higher latitudes, in both the northern and
southern hemispheres, sin α changes until it reaches a maximum or minimum of ±1 .
Such a point is known as a vertex and the latitude φ will take maximum value φ0 .
P0 vertex
φ max •
s
P
φ •
∆λ
A B A'
• equator •
•
node λ λ0 node node A
∆ λ4

φ min •
vertex

Fig. 3: Schematic diagram of the oscillation of a geodesic on an ellipsoid

Thus the geodesic oscillates over the surface of the ellipsoid between two
parallels of latitude having a maximum in the Northern and Southern Hemispheres and
crossing the equator at nodes. As we will demonstrate later, due to the eccentricity of
the ellipsoid, the geodesic will not repeat after a complete revolution.
Figure 3 shows a schematic diagram of the oscillation of a geodesic on an
ellipsoid. P is a point on a geodesic that crosses the equator at A, heading in a north-
easterly direction reaching a maximum northerly latitude φmax at the vertex P0 (north),
then descends in a south-easterly direction crossing the equator at B, reaching a
maximum southerly latitude φmin at P0 (south), then ascends in a north-easterly
direction crossing the equator again at A′ . This is one complete revolution of the
geodesic, but λA′ does not equal λA due to the eccentricity of the ellipsoid. Hence we
say that the geodesic curve does not repeat after a complete revolution.

5
EQUATIONS FOR COMPUTATION ALONG GEODESICS
Using Clairaut's equation and simple differential relationships, expressions for
distances s and longitude differences ∆λ (see Figure 3) between P on a geodesic and
the vertex P0 can be obtained. These expressions are in the form of elliptic integrals,
which by their nature do not have exact (or closed) solutions.
Expanding the integrands into infinite series, integrating term-by-term and then
truncating to a finite number of terms is the usual technique to obtain working solutions
for s and ∆λ (e.g., Thomas 1970). In this section, we show how this method can be
simplified by using recurrence relationships to generate solutions to the integrals in the
series. Our relationships are slightly different from Pittman (1986) and our notation is a
little different but in all other respects, we have followed his elegant approach.

Relationships between parametric latitude ψ and geodetic latitude Φ

Development of formulae is simplified if parametric latitude ψ is used rather
than geodetic latitude φ . The connections between the two latitudes can be obtained
from the following relationships.
Figure 4 shows a portion of a meridian
NPE of an ellipsoid having semi-major axis z
auxiliary circle
OE = a and semi-minor axis ON = b . P is a N'
point on the ellipsoid and Q is a point on an Q
N
auxiliary circle centred on O of radius a. P and •
Q have the same perpendicular distance from the • P
ta n
axis of revolution ON. The normal to the b ge
nt
ellipsoid at P cuts the major axis at an angle φ
l
ma
a

(the geodetic latitude) and intersects the

n or

rotational axis at H and the distance PH = ν .

The angle QOE = ψ is the parametric latitude. O ψ φ M E
a w
The Cartesian equation of the ellipse is
w2 z 2 H•
+ = 1 and the Cartesian equation of the
a 2 b2
auxiliary circle is w2 + z 2 = a 2 . We may re- Fig. 4: Meridian section of ellipsoid
arrange both equations so that w2 is on the left-hand side of the equals sign giving
a2
w2 = a 2 − 2 z 2 (ellipse) and w2 = a 2 − z 2 (circle). Now, since the w-coordinates of P
b
a2 2 b
and Q are the same then a − 2 z P = a 2 − zQ2 which leads to z P = zQ .
2

b a
Using this relationship
w = OM = a cosψ and z = MP = b sin ψ (6)
dw dz
and differentiating equations (6) with respect to ψ gives = − a sinψ , = b cosψ
dψ dψ
dz dz dψ b
and the chain rule gives = = − cotψ .
dw dψ dw a

6
 dz
Now by definition, is the gradient of the tangent and from Figure 4 we may write
dw
= − tan ( 90 − φ ) = − cot φ . Equating the two expressions for dz dw gives a
dz
dw
relationship between ψ and φ as
b
tanψ = tan φ = (1 − f ) tan φ (7)
a
From equation (6) and Figure 4, w = a cosψ = ν cos φ and using equation (3) gives
cos φ
cosψ = (8)
(1 − e sin 2 φ )
2 12

Alternatively, using the trigonometric identity sin 2 A + cos 2 A = 1 , equation (8) can be
written as
sinψ
sin φ = (9)
(1 − e cos 2 ψ )
2 12

The latitudes Φ0 and ψ0 of the geodesic vertex

Denoting the latitude of the vertex as φ0 (a maximum), Clairaut's equation (5) gives

ν 0 cos φ0 = constant = ν cos φ sin α (10)

Denoting the parametric latitude of the vertex as ψ 0 and using a cosψ = ν cos φ from
before, equation (10) becomes a cosψ 0 = a cosψ sin α and ψ 0 is defined as

cosψ 0 = cosψ sin α (11)

Squaring both sides of equation (11) and using again the identity sin A + cos A = 1 we
2 2

can obtain the azimuth α of a geodesic as

cos 2 ψ − cos 2 ψ 0
cos α = (12)
cosψ
From equation (11) we see that if the azimuth α of a geodesic is known at P having
parametric latitude ψ , the parametric latitude ψ 0 of the vertex P0 can be computed.
Conversely, given ψ and ψ 0 of points P and P0 the azimuth of the geodesic between
them may be computed from equation (12).
ds dλ
In the following sections, two differential equations; one for and the other for ,
dψ dψ
will be developed that will enable solutions for the geodesic distance s and the
longitude difference ∆λ between P and the vertex P0 .

7
 ds dλ
Differential equations for distance and longitude difference
dψ dψ
From equation (9) we may write sin 2 ψ = (1 − e 2 cos 2 ψ ) sin 2 φ and differentiating
implicitly and re-arranging gives

dφ (1 − e sin φ ) sinψ cosψ

2 2

= (13)
dψ (1 − e 2 cos 2 ψ ) sin φ cos φ
ds
Using the chain rule and equation (4) gives an expression for the derivative as
dψ

ds ds dφ ρ (1 − e sin φ ) sinψ cosψ

2 2

= = (14)
dψ dφ dψ cos α (1 − e 2 cos 2 ψ ) sin φ cos φ
b2
Using equations (7), (8), (9) and the fact that 1 − e 2 = , we may write
a2

(1 − e2 cos2 ψ )
12
ds
= a cosψ (15)
dψ ( cos2 ψ − cos2 ψ 0 )
12

Similarly, the chain rule and equations (4) and (15) gives

(1 − e2 cos2 ψ )
12
d λ d λ ds sin α
= = a cosψ (16)
dψ ds dψ ν cos φ ( cos2 ψ − cos2 ψ 0 )
12

Using equation (10) and the relationship a cosψ = ν cos φ , we may write

d λ cosψ 0 (1 − e cos ψ )
2 2 12

= (17)
dψ cosψ ( cos 2 ψ − cos 2 ψ )1 2
0

Equations (15) and (17) are the basic differential equations that will yield solutions for
distance s and longitude difference ∆λ along the geodesic curve between P and the
vertex P0 .

Formula for computing geodesic distance s between P and the vertex P0

Equation (15) can be simplified by letting u = sinψ and u0 = sinψ 0 , so that
du
= cosψ and cos 2 ψ − cos 2 ψ 0 = u02 − u 2 , hence
dψ

du (1 − e cos ψ )
2 2 12
ds
=a (18)
dψ dψ ( u 2 − u 2 )1 2
0

8
 du a (1 − e cos ψ )
2 2 12
ds ds
The chain rule gives = = but using cos 2 ψ = 1 − sin 2 ψ
du dψ dψ (u 2 − u 2 )
12
0

and equations (1) and (2) we are able to obtain, after some manipulation

ds b (1 + ε u )
2 12

= (19)
du ( u 2 − u 2 )1 2
0

where ε = ( e′ ) . The geodesic distance s between P and the vertex P0 is given by

2

p = u0
(1 + ε p ) 2 12

s=b ∫ u −p dp (20)
p =u ( )
2
0
2 12

where sinψ ≤ p ≤ sinψ 0 . Equation (20) can be simplified by use of the binomial series
and the numerator of the integrand is given by
∞
(1 + ε p )
2 12
= ∑ Bn2 ( ε p 2 )
1 n
(21)
n=0

1
where Bn2 are binomial coefficients computed from the recurrence relationship
3 − 2n 12
Bn2 =
Bn −1 , n ≥ 1 and B02 = 1
1
1
(22)
2n
Equation (20) can now be written as
u0 ∞ ∞ u0 ∞
1 p 2n
s = b∫ ∑B ε p dp = b∑ Bn ε ∫ dp = b∑ ε n Bn2 I n
1 1 1
n 2n n
2 2
(23)
(u ) (u )
n
2 12 2 12
u
2
0 −p n =0 n =0 u
2
0 −p n =0

u0
p 2n
where In = ∫ dp , for n ≥ 0 (24)
(u − p2 )
2 12
u 0

The solution of the integral I n is fundamental to the computation of the distance s

along the geodesic between P and P0 , and the usual technique is to find solutions for
each integral I n and expand equation (23) into a finite series; e.g. Thomas (1970, pp.
33-34). Pittman's (1986) approach, outlined below, was to developed the integral I n as
a recurrence equation having the general form I n = an −1 + bn −1 I n −1 where the coefficients
an −1 and bn −1 are functions of n, ψ and ψ 0 and an initial value of I 0 is a function of
ψ and ψ 0 only.
−p
u0 u0 u0
p 2n
∫ dp = − ∫ p 2 n −1 dp = − ∫ p 2 n −1
d 2
( u0 − p 2 ) dp
12
Now In =
u (u 2
0 −p )
2 12
u (u 2
0 −p )
2 12
u
dp
and using integration by parts (e.g., Ayres 1972) the integral I n becomes

9
 p = u0
I n = − ⎡⎢ p 2 n −1 ( u02 − p 2 ) − ∫ ( u02 − p 2 ) ( 2n − 1) p 2 n − 2 dp ⎤⎥
12 12

⎣ ⎦ p =u
p 2n−2
u0

= u 2 n −1 ( u02 − u 2 ) + ( 2n − 1) ∫ ( u02 − p 2 )
12
dp
u (u 2
0 −p 2 12
)
= u 2 n −1 ( u02 − u 2 ) + ( 2n − 1) ⎡⎣u02 I n −1 − I n ⎤⎦
12
(25)
and

2n I n = u 2 n −1 ( u02 − u 2 ) + ( 2n − 1) u02 I n −1
12
for n = 1, 2,3,… (26)

so that u = Uu0 , u02 − u 2 = u02 (1 − U 2 ) giving

u
Let U =
u0

2n I n = (Uu0 ) u0 (1 − U 2 ) + ( 2n − 1) u02 I n−1

2 n −1 12
for n = 1, 2,3,… (27)
2n I n 2 ( n − 1) u02
Let J n = so that J n −1 = I n −1 and the recurrence formula for I n becomes
u02 n u02 n
a simpler recurrence formula for J n
2n − 1
J n = U 2 n −1 1 − U 2 + J n −1 for n = 2, 3,… (28)
2 ( n − 1)
with initial condition
2 I1
J1 =
= U 1 − U 2 + I0 (29)
u02
I 0 has a simple result derived from equation (24) as follows:

( )
u0
2 −1 2
I 0 = (1 u0 ) ∫ 1 − [ p u0 ] dp (30)
u

and with the transformation p = u0 cos θ , dp dθ = −u0 sin θ and 1 − [ p u0 ] = 1 − cos 2 θ

2

0
⎛u⎞
I0 = ∫ ( −1) dθ = arccos ⎜ ⎟ = arccos U (31)
⎛u ⎞
θ = arccos⎜ ⎟
⎝ u0 ⎠
⎝ u0 ⎠

Using these results, the distance s along the geodesic between P and the vertex P0 is

⎧ 1 ∞ 1 ⎫
s = b ⎨ I 0 + ∑ ε nu02 n Bn2 J n ⎬
1

⎩ 2 n =1 n ⎭
b b b
= bI 0 + ε u02 B12 J1 + ε 2u04 B22 J 2 + ε 3u06 B32 J 3 +
1 1 1

2 4 6
= D0 + D1 + D2 + D3 + (32)

10
Formula for computing difference in longitude ∆λ between P and P0
Using the binomial series we may write equation (17) as
dλ ∞
cos 2 n −1 ψ
= cosψ 0 ∑ ( −1) e 2 n Bn2
n 1
(33)
dψ ( cos2 ψ − cos2 ψ 0 )
12
n=0

and the difference in longitude between P and the vertex P0 is

ψ0
dλ ∞
∆λ = ∫ dθ = cosψ 0 ∑ ( −1) e 2 n Bn2 Ln
n 1
(34)
θ =ψ dθ n=0

where the integral Ln is

ψ0
cos 2 n θ
Ln = ∫
θ =ψ cos θ ( cos 2 θ − cos 2 ψ 0 )
12
dθ , n≥0 (35)

Again, let u = sinψ , u0 = sinψ 0 and put p = sin θ . Then dθ dp = sec θ ,

cos 2 θ = 1 − p 2 , and with

( cos 2 θ ) (1 − p2 )
n n
cos 2 n θ
dp = (1 − p 2 ) dp
n −1
dθ = cos θ dθ =
cos θ cos θ2
1− p 2

and

( cos θ − cos ψ ) (
= 1 − sin 2 θ − (1 − sin 2 ψ 0 ) ) = ( u02 − p 2 )
12 12 12
2 2
0

giving

Ln =∫
(1 − p )
u0 2 n −1

dp, n ≥ 1 (36)
(u − p )
u
2
0
2 12

Using the binomial series, the numerator of the integrand can be expanded into a
n −1
polynomial (1 − p 2 )
n −1
= ∑ ( −1) Bmn −1 p 2m , where the binomial coefficients Bmn −1 are
m

m =0

given by
n − m n −1
Bmn −1 =
Bm −1 for m = 2,3, 4,… (37)
m
with an initial value B1n −1 = n − 1 and noting that B0n −1 = 1 .
Using these results, equation (36) becomes
n −1 u0 n −1
p 2m
Ln = ∑ ( −1) Bmn −1 ∫ dp = ∑ ( −1) Bmn −1 I m
m m
(38)
m=0 u ( u02 − p )
2 12
m =0

u0
p 2m
where Im = ∫ dp , for m ≥ 0 (39)
( u02 − p 2 )
12
u

11
and equation (39) is the same as equation (24) except for a change of index variable.
Using this similarity and the expressions above, the longitude difference given by
equation (34) can be expressed as
⎧ ∞ n −1
⎫
∆λ = cosψ 0 ⎨ L0 + ∑ ( −1) e 2 n Bn2 ∑ ( −1) Bmn −1 I m ⎬
n 1 m
(40)
⎩ n =1 m =0 ⎭
Equation (40) can expanded as
⎧ ⎡ ∞
1 ⎤
∆λ = cosψ 0 ⎨ L0 + ⎢ −e2 B12 + ∑ ( −1) e2 n Bn2 ⎥ I 0
1 n

⎩ ⎣ n=2 ⎦
(41)
∞ n −1
⎫
+ ∑ ( −1) e Bn ∑ ( −1) n −1
n 1 m
2n 2
B I ⎬
m m
n=2 m =1 ⎭
and then simplified by use of the binomial series, where
∞ ∞ ∞
(1 − e2 ) = ∑ ( −1) e 2 n Bn2 = 1 + ∑ ( −1) e 2 n Bn2 = 1 − e 2 B12 + ∑ ( −1) e 2 n Bn2
12 n 1 n 1 1 n 1
(42)
n=0 n =1 n=2

The terms in [ ] of equation (41) are the last two terms on the right-hand side of
equation (42) and using this equivalence gives
⎧
( ) ⎫
∞ n −1
∆λ = cosψ 0 ⎨ L0 + 1 − e 2 − 1 I 0 + ∑ ( −1) e 2 n Bn2 ∑ ( −1) Bmn −1 I m ⎬
n 1 m

⎩ n=2 m =1 ⎭
⎧⎪ ( −1) ⎫⎪
( )
m
∞ n −1
= cosψ 0 ⎨ L0 + 1 − e −1 I0 + ∑ ( −1) e Bn ∑ u02 m Bmn −1 J m ⎬
n 1
2 1 2n
2
2
(43)
⎪⎩ n=2 m =1 m ⎪⎭
where I 0 is obtained from equation (31) and J m are given by equation (28), noting that
2m
as before J m = 2 m I m .
u0
A simple expression for L0 is obtained from equation (35) as follows
ψ0 ψ0
1 sec 2 θ
L0 = ∫ cos θ ( cos 2 θ − cos 2 ψ 0 )
dθ = ∫ dθ (44)
( sin ψ 0 − tan θ cos ψ 0 )
12 2 2 2 12
θ =ψ θ =ψ

Putting x = cotψ 0 tan θ then dθ dx = tanψ 0 cos 2 θ and

⎛ cos 2 ψ 0 ⎞
sin ψ 0 − tan θ cos ψ 0 = sin ψ 0 ⎜ 1 − tan θ
2 2 2 2
⎟
2

⎝ sin 2 ψ 0 ⎠
= sin 2 ψ 0 (1 − tan 2 θ cot 2 ψ 0 )
= sin 2 ψ 0 (1 − x 2 )
so that
tanψ 0
1
dx
L0 =
sinψ 0 ∫ψ
tan 1 − x2
(45)
x=
tanψ 0

12
 dx ⎧ arcsin x
since ∫ 1− x
= ⎨π
2
⎩ 2 − arccos x
, then using the second result gives

1
dx ⎛ tanψ ⎞
L0 = secψ 0 ∫ψ
tan 1− x
= secψ 0 arccos ⎜
2
⎟
⎝ tanψ 0 ⎠
(46)
x=
tanψ 0

Equation (40) can be simplified further to give the longitude difference ∆λ between P
and the vertex P0 as

∆λ = cosψ 0 {M 0 + M 1 + M 2 + M 3 + } (47)

⎧L for n = 0
⎪ 0
where
⎪
⎪ 1
(
M n = ⎨ 1 − e2 − 1 I 0 ) for n = 1 (48)

⎩⎪ 2 Bn ( −1) e K n for n ≥ 2
1 n 2n
2

( −1)
m
n −1
and Kn = ∑ u02 m Bmn −1 J m for n = 2,3, 4,… (49)
m =1 m

A GEODESIC ON AN ELLIPSOID DOES NOT REPEAT AFTER A SINGLE

REVOLUTION
Earlier, it was mentioned that due to the eccentricity of the ellipsoid, the geodesic will
not repeat after a complete revolution. Here is a demonstration of that fact.
When P is at the node A of Figure 3 then ∆λ = ∆λ4 and using equation (17) we have

(1 − e cos θ )
12
ψ0 2 2

4 ( ∆λ4 ) = 4 cosψ 0 ∫ cos θ cos θ − cos ψ dθ (50)

( )
2 2 12
θ =0 0

Since this integral is difficult to evaluate, we instead determine upper and lower bounds
for the quantity 4 ( ∆λ4 ) by using the bounds of the integration variable θ . This allows
certain terms within the integral to be disposed of and a simplified integral evaluated.

For 0 ≤ θ ≤ ψ 0 , the bounds on the numerator of the integrand are

(1 − e )2 12
≤ (1 − e2 cos 2 θ ) ≤ (1 − e2 cos 2 ψ 0 )
12 12
so that on the one hand

(1 − e cos ψ )
12
ψ0 2 2

4 ( ∆λ ) ≤ 4 cosψ ∫ dθ
0

cos θ ( cos θ − cos ψ )

4 0 12
2 2
θ =0 0

= 4 cosψ (1 − e cos ψ ) L 2 2 12
0 0 0 ψ =0

= 4 cosψ (1 − e cos ψ ) π secψ 2 2 12

1
0 0 2 0

= 2π (1 − e cos ψ ) 2 2 12
0 (51)

13
while on the other hand

4 ( ∆λ4 ) ≥ 4 cosψ 0
ψ0
(1 − e )
2 12

∫
θ =0 cos θ ( cos 2 θ − cos 2 ψ 0 )
12
dθ

= 2π (1 − e 2 )
12
(52)
Combining these inequalities gives the bounds for the quantity 4 ( ∆λ4 ) as

2π (1 − e2 ) ≤ 4 ( ∆λ4 ) ≤ 2π (1 − e 2 cos 2 ψ 0 )
12 12
(53)
Therefore, after a single revolution, 4 ( ∆λ4 ) < 2π when 0 < ψ 0 < 90 . Note that when
ψ 0 = 0 the geodesic is an arc of the equator (a circle) and when ψ 0 = 90 the geodesic
is an arc of the meridian (an ellipse).

NUMERICAL RESULTS FOR DISTANCE AND LONGITUDE EQUATIONS

Equations (32) and (47) for computing distance s and longitude difference ∆λ between
P and the vertex P0 are relatively simple summations of terms. To test the number of
terms required for accurate answers, a geodesic was chosen with an azimuth
α = 43 12′ 36′′ at P having latitude φ = 9 35′ 24′′ on the ellipsoid of the Geodetic
Reference System 1980 (GRS80) (Moritz 1980), defined by a = 6378137 metres and
f = 1 298.257 222101 .

Numerical constants for GRS80 ellipsoid and geodesic

b = a (1 − f ) = 6356752.314140356 metres
ψ = arctan ⎡⎣(1 − f ) tan φ ⎤⎦ = 0.166826262923 radians
ψ 0 = arccos [ cosψ sin α ] = 0.829602797993 radians
u = sinψ = 0.166053515348; u0 = sinψ 0 = 0.737663250899
u sinψ
U = = = 0.225107479796; I 0 = arccos U = 1.343742980976 radians
u0 sinψ 0
tanψ
V = = 0.154125311675; L0 = secψ 0 arccos V = 2.097333540996 radians
tanψ 0

14
 Table 1: Ellipsoid and geodesic constants and binomial coefficients for
equations (32) and (47)
1
n e2 n εn u02 n Bn2
1 6.694380022901e-003 6.739496775479e-003 0.544147071727 0.500000000000
2 4.481472389101e-005 4.542081678669e-005 0.296096035669 -0.125000000000
3 3.000067923478e-007 3.061134482735e-007 0.161119790759 0.062500000000
4 2.008359477428e-009 2.063050597570e-009 0.087672862339 -0.039062500000
5 1.344472156450e-011 1.390392284997e-011 0.047706931312 0.027343750000
6 9.000407545482e-014 9.370544321391e-014 0.025959586974 -0.020507812500
7 6.025214847044e-016 6.315275323850e-016 0.014125833235 0.016113281250
8 4.033507790574e-018 4.256177768135e-018 0.007686530791 -0.013092041016

Table 2: Recurrence formula values and distance components for equation (32)

n Jn Dn
1 1.563072838216 8.541841303930e+006 8541841.303930 m
2 2.355723441968 9.109578467516e+003 9109.5784675
3 2.945217495733 -6.293571169346e+000 -6.2935712
4 3.436115617261 9.618619108010e-003 0.0096186
5 3.865631515581 -1.929070816523e-005 -0.0000193
6 4.252194740421 4.456897529564e-008 0.0000000
7 4.606544305836 -1.123696751599e-010 -0.0000000
8 4.935583185013 3.006580650377e-013 0.0000000
sum 8.550944598425e+006 s = 8550944.598425 m

Table 3: Recurrence formula values and longitude components for equation (47)
n Mn
Jn
0 2.097333540996e+000
1 1.563072838216 -4.505315819380e-003
2 2.355723441968 2.382298926901e-006
3 2.945217495733 1.267831357153e-008
4 3.436115617261 6.525291638252e-011
5 3.865631515581 3.431821056093e-013
6 4.252194740421 1.852429353592e-015 ∆λ = cosψ 0 ( sum ) ≅ 1.413013969112 radians
7 4.606544305836 1.023576994037e-017
= 80.959736823113 degrees
8 4.935583185013 5.769507252421e-020
sum 2.092830620219e+000 = 80 57′ 35.052563′′

Inspection of these numerical values indicates than an upper limit of N = 8 in the

summations is more than sufficient for accuracies of 0.000001 metre in distances and
0.000001 second of arc for longitude differences. [Results for s and ∆λ can be
confirmed using Vincenty's equations (Vincenty 1975) that have been programmed in a
Microsoft™ Excel workbook that can be downloaded from the website of Geoscience
Australia at http://www.ga.gov.au/]

15
 It should be noted here that the distance and longitude equations [equations (32)
and (47)] are not themselves, solutions to the direct or inverse problems. Instead, they
are the basic tools, which if used in certain ways, enable the solution to those problems.
In a computer program, equations (32) and (47) would be embedded in a function
that returned s and ∆λ given the ellipsoid parameters ( a, f ) , parametric latitudes
(ψ ,ψ 0 ) and the upper limit of summations ( N ) . A brief explanation of how such a
function might be used is given below.

USING THE DISTANCE AND LONGITUDE EQUATIONS TO COMPUTE THE

DIRECT AND INVERSE PROBLEM

α2
P2 •
α1 = α12 •
∆λ 2 α21
s2 P0 vertex P1 •
φmax P2 •
•
s s1
P1 s4
φ1 •
s3
∆λ1
A B equator A'
• • •
node λ1 λ0 node node A

∆λ4

φmin •
vertex

Fig. 5: Schematic diagram of a geodesic between P1 and P2 on an ellipsoid

Direct solution
The key here is to use the distance equation in an iterative computation of sinψ 2 . Once
this is known, then φ2 , λ2 and α 21 follow. The steps in the computation are:

1. Test the azimuth to determine whether the geodesic is heading towards or away
from the nearest vertex P0 , noting that P0 will be in the same hemisphere as P1 .
2. Compute ψ 1 and ψ 0 ; then use the distance and longitude equations to compute s1
and ∆λ1 between P1 and P0 , as well as λ0 . (see Fig. 5).
3. With u = sinψ = 0 , compute s4 and ∆λ4 between the node and P0 .

16
 ⎧ s − s1 if geodesic is heading towards P0
4. Compute s2 = ⎨ . If s2 > 0 then P2 is
⎩ s + s1 if geodesic is heading away from P0
after P0 and closer to another vertex P0′ in which case s2 is reduced by multiples of
2s4 until s2 < s4 and the number of vertices n determined (vertices are 2s4 apart).
If s2 < 0 then P2 is before P0 . (Note that in Fig. 5, s2 < 0 and P2 is before P0 )
5. Compute ψ 2 by iteration. An approximate value ψ 2′ is found from equations (32)
s ⎛ sinψ ⎞
by taking the first term only; hence = I 0 = arccos ⎜ ⎟
b ⎝ sinψ 0 ⎠
⎛s ⎞
and sinψ 2′ = sinψ 0 cos ⎜ 2 ⎟ .
⎝b⎠
ds u02 − u 2
Now a re-arrangement of the differential equation (19) gives du =
b 1+ ε u2
where u = sinψ 2′ , ds = s2′ − s2 and s2′ is computed from the distance equation with
the approximate parametric latitude ψ 2′ . Equation (19), linking ds and du, is the
basis of the iterative solution for sinψ 2 (and hence φ2 ).
6. After computing ψ 2 the longitude difference ∆λ2 is computed and depending on
the number of vertices and the direction of the geodesic, λ2 is determined. The
azimuth α 2 follows from Clairaut's equation and the reverse azimuth α 21 obtained.

Inverse solution
This is the more difficult of the two solutions since ψ 0 is unknown and must be
determined by iteration, using approximations for s, α1 and α 2 obtained by
approximating the ellipsoid with a sphere and using spherical trigonometry. The steps
in the computation are:

1. Convert longitudes of P1 and P2 to east longitudes in the range 0 < λ1 , λ2 < 360
and determine a longitude difference ∆λ in the range −180 ≤ ∆λ ≤ 180 . ±∆λ
corresponding to east/west direction of the geodesic from P1 .
2. Compute parametric latitudes ψ 1 and ψ 2 then use these and ∆λ as latitudes and
longitude difference on a sphere to compute spherical distance σ and spherical
angles β1 and β 2 . These can be used to determine approximations of s and α12 .
3. Compute ψ 0 by iteration. Approximations ∆λ1′ and ∆λ2′ can be obtained from
⎛ tanψ ⎞
equation (47) noting that M 0 = secψ 0 arccos ⎜ ⎟ and ignoring terms
⎝ tanψ 0 ⎠
M 1 , M 2 , M 3 ,…
⎛ tanψ 1 ⎞ ⎛ tanψ 2 ⎞
This gives ∆λ1′ = arccos ⎜ ⎟ and ∆λ2′ = arccos ⎜ ⎟ , and
⎝ tanψ 0 ⎠ ⎝ tanψ 0 ⎠

17
 ⎧⎪ ⎛ tanψ 1 ⎞ ⎛ tanψ 2 ⎞ ⎫⎪
f (ψ 0 ) = ∆λ ′ − ∆λ = ⎨± arccos ⎜ ⎟ ± arccos ⎜ ⎟ ± ∆λ4′ ⎬ − ∆λ where the ±
⎪⎩ ⎝ tanψ 0 ⎠ ⎝ tanψ 0 ⎠ ⎪⎭
signs are associated with the east/west direction of the geodesic.
ψ 0 can be found using Newton's iterative method (Williams 1972)
f (ψ 0 )
(ψ 0 )n +1 = (ψ 0 )n − (54)
f ′ (ψ 0 )
where f ′ (ψ 0 ) is the derivative of f (ψ 0 ) . An initial value of ψ 0 can be computed
from equation (11).
4. Once ψ 0 is known then s1 , ∆λ1 ; s2 , ∆λ2 and s4 , ∆λ4 can be computed from the
distance and longitude equations and s obtained. The forward and reverse azimuths
can be found from Clairaut's equation (5).

CONCLUSION
Pittman's (1986) recurrence relationships for evaluating integrals allow beautifully
compact equations for distance s and longitude difference ∆λ along a geodesic between
P and the vertex P0 . These equations can be easily translated into a computer program
function returning s and ∆λ given a, f, u and u0 . Using such a function, algorithms (as
outlined above), can be constructed to solve the direct and inverse problems on the
ellipsoid. Pittman's (1986) paper (which included FORTRAN computer code) has a
concise development of the necessary equations and algorithms. The paper here has a
more detailed development of the recurrence relationships (with a slightly different
formulation) as well as additional information on the definition and properties of a
geodesic.
Interestingly, Pittman's (1986) method is entirely different to other approaches
that fall (roughly) into two groups: (i) numerical integration techniques and (ii) series
expansion of integrals; the latter of these with a history of development extending back
to Bessel's (1826) method. Numerical integration, a technique made practical with the
arrival of computers in the mid to late 20th century, is relatively modern. So too is
Pittman's method.
To our knowledge, this is the first paper (since the original) discussing his elegant
method; a method that has much to recommend it, and one that we hope might become
the object of study in undergraduate surveying courses and discussion in the geodetic
literature.

REFERENCES
Ayres, F., 1972. Calculus, Schaum's Outline Series, Theory and problems of
Differential and Integral Calculus, 2nd edn, McGraw-Hill Book Company, New
York.
Bessel, F. W., 1826, 'Uber die Berechnung der Geographischen Langen und Breiten aus
geodatischen Vermessungen. (On the computation of geographical longitude and
latitude grom geodetic measurements)', Astronomische Nachrichten (Astronomical
Notes), Band 4 (Vol. 4), No. 86, Spalten 241-254 (Columns 241-254).

18
Bowring, B. R., 1972, 'Correspondence: Distance and the spheroid', Survey Review,
Vol. 21, No. 164, pp. 281-284.
Bowring, B. R., 1983, 'The geodesic inverse problem', Bulletin Geodesique, Vol. 57,
No. 2, pp. 109-120.
Bowring, B. R., 1984, 'Note on the geodesic inverse problem', Bulletin Geodesique,
Vol. 58, p. 543.
DSB, 1971. Dictionary of Scientific Biography, C.C. Gillispie (Editor in Chief),
Charles Scribener's Sons, New York.
Jank, W., Kivioja, L.A., 1980, 'Solution of the direct and inverse problems on reference
ellipsoids by point-by-point integration using programmable pocket calculators',
Surveying and Mapping, Vol. 15, No. 3, pp. 325-337.
Kivioja, L. A., 1971, 'Computation of geodetic direct and indirect problems by
computers accumulating increments from geodetic line elements', Bulletin
Geodesique, No. 99, pp. 55-63.
Lauf, G.B., 1983. Geodesy and Map Projections, TAFE Publications Unit,
Collingwood, Australia
Moritz, H., 1980, 'Geodetic reference system 1980', The Geodesists Handbook 1980,
Bulletin Geodesique, Vol. 54, No. 3, pp. 395-407.
Pittman, M.E., 1986. 'Precision direct and inverse solutions of the geodesic', Surveying
and Mapping, Vol. 46, No. 1, pp. 47-54, March 1986.
Rainsford, H. F., 1955, 'Long geodesics on the ellipsoid', Bulletin Geodesique, No. 37,
pp. 12-22.
Saito, T., 1970, 'The computation of long geodesics on the ellipsoid by non-series
expanding procedure', Bulletin Geodesique, No. 98, pp. 341-374.
Saito, T., 1979, 'The computation of long geodesics on the ellipsoid through Gaussian
quadrature', Bulletin Geodesique, Vol. 53, No. 2, pp. 165-177.
Sjöberg, Lars E., 2006, 'New solutions to the direct and indirect geodetic problems on
the ellipsoid', Zeitschrift für Geodäsie, Geoinformation und Landmanagement (zfv),
2006(1):36 pp. 1-5.
Struik, D.J., 1933. 'Outline of a history of differential geometry', Isis, Vol. 19, No.1, pp.
92-120, April 1933. (Isis is an official publication of the History of Science Society and has
been in print since 1912. It is published by the University of Chicago Press - Journals Division:
http://www.journals.uchicargo.edu/)
Thomas, P.D., 1970. Spheroidal Geodesics, Reference Systems, & Local Geometry,
Special Publication No. 138 (SP-138), United States Naval Oceanographic office,
Washington.
Thomas, C. M. and Featherstone, W. E., 2005, 'Validation of Vincenty's formulas for
the geodesic using a new fourth-order extension of Kivioja's formual', Journal of
Surveying Engineering, Vol. 131, No. 1, pp. 20-26.
Vincenty, T., 1975, 'Direct and inverse solutions on the ellipsoid with application of
nested equations', Survey Review, Vol. 22, No. 176, pp. 88-93.
Vincenty, T., 1976, 'Correspondence: solutions of geodesics', Survey Review, Vol. 23,
No. 180, p. 294.
Williams, P. W., 1972, Numerical Computation, Thomas Nelson and Sons Ltd,
London.

19