International Journal
on Marine Navigation
and Safety of Sea Transportation
Volume 5
Number 4
December 2011
527
1 INTRODUCTION
The navigational problem is as follows (Fig. 1):
− we have known geographic coordinates of posi-
tion P
′
– φ
P′
, λ
P′
,
− we have ranges, bearings or courses from position
P
′
to position P,
− we search for geographic coordinates of position
P – φ
P
, λ
P
,
The most common solution of such a navigational
problem is a rather strange combination of flat and
ellipsoidal calculations:
− conversion ranges, bearings and courses, by sol-
ving flat triangles, to δx, δy increments in a flat
rectangular coordinate frame (with the y-axis
pointing north),
− conversion rectangular δx, δy flat increments to
geographic coordinates increments δφ, δλ, on the
reference ellipsoid, by the equations
)(R
y
PM
′
ϕ
δ
=δϕ
(1)
)cos()(R
x
PPN
′′
ϕϕ
δ
=δλ
(2)
where R
M
(φ
P′
) = the radius of curvature in meridian
for P′; and R
N
(φ
P′
)
= the radius of curvature in the
prime vertical for P′
given by the equations
3
P
22
2
0
M
)sine1(
)e1(a
R
′
ϕ−
−
=
(3)
P
22
0
N
sine1
a
R
′
ϕ−
=
(4)
where a
0
= the semi-major axis of the reference el-
lipsoid; and e = eccentricity
and finally
δϕ+ϕ=ϕ
′
PP
(5)
δλ+λ=λ
′
PP
(6)
if east longitude and north latitude are considered
positive and west longitude and south latitude are
considered negative.
Apart from the obvious errors of assuming the
part of the ellipsoid a plane there are also errors hid-
den in Equations 1 and 2 – although Equations 3 and
4 are accurate for the ellipsoid – the errors of as-
suming the main radii of curvature constant at points
P′ and P.
Solutions of Direct Geodetic Problem in
Navigational Applications
A.S. Lenart
Gdynia Maritime University, Gdynia, Poland
any simplifications for a plane or a sphere, by an application of solutions of direct geodetic problem are pre-
sented. The rigorous, rapid, non-iterative solution of the direct geodetic problem according to Sodano, for any
length of geodesics, is attached.
528
Figure 1.Definition of the problem and the position from range
and bearings
2 ERROR OF ASSUMING THE PART OF
ELLIPSOID TO BE A PLANE
The spherical excess in an equilateral spherical trian-
gle with sides d and spherical radius R is
2
2
R4
3d
=ε
(7)
For d = 22 km (≈ 12 n.m.) and R = 6370 km, we
get
rad 105
-6
×≈ε
This gives linear changes in the range of 11 cm.
3 ERRORS OF ASSUMING THE MAIN RADII
CONSTANT AT POINTS P′ AND P
A better approximation of Equation 1 should be
(Lenart 1985)
ϕ
δϕ+ϕ
δ
=
ϕ+ϕ
δ
=δϕ
′
′
d
dR
2
1
)(R
y
2
R
y
M
P
M
PP
M
*
(8)
Defining the resulting error as
*
δϕ−δϕ=δϕ∆
(9)
we get
ϕ
δϕ+ϕϕ
ϕ
δϕ
δ=δϕ∆
′′
d
dR
2
1
)(R)(R
d
dR
2
1
y
M
PMPM
M
(10)
which with
ϕ
δϕ+ϕδϕ=δ
′
d
dR
2
1
)(Ry
M
PM
(11)
yields
)(R
d
dR
)(
2
1
PM
M
2
′
ϕ
ϕ
δϕ
=δϕ∆
(12)
After substitution
M
2
M
R2sine
2
3
d
dR
ϕ≈
ϕ
(13)
we finally get
ϕδϕ≈δϕ∆ 2sin)(e
4
3
22
(14)
where, if δφ is in radians the result is also in radians.
For example, if
°=ϕ
×≈
′′
=δϕ
−
45
rad103521
4
we get
cm 396.012
rad 106.125 rad 101225
150
1
4
3
8-8-
≈
′′
≈
×≈×××≈δϕ∆
In the case of error in longitude we have, in ac-
cordance with Equation 2,
δϕ
+ϕ
ϕ
δϕ+ϕ
δ
=δλ
′′
2
cos
d
dR
2
1
)(R
x
P
N
PN
*
(15)
Then
*
δλ−δλ=δλ∆
(16)
and, after substitution
PPP
sin
2
cos
2
cos
′′′
ϕ
δϕ
−ϕ≈
δϕ
+ϕ
(17)
δϕ
+ϕ
ϕ
δϕ+ϕδλ=
δ
′′
2
cos
d
dR
2
1
)(Rx
P
N
PN
(18)
After simplification
PN
PNP
N
cosR
sinRcos
d
dR
2
1
′
′′
ϕ
ϕ−ϕ
ϕ
δϕδλ
≈δλ∆
(19)
Since
N
2
N
R2sine
2
1
d
dR
ϕ≈
ϕ
(20)
then finally
d
P'
P
N
529
ϕδϕδλ−≈δλ∆ tan
2
1
(21)
If δφ and δλ are in radians, the result is also in ra-
dians.
For example, if
)equatortheon
21torelateswhich(rad1020196
rad103521
80
4
4
′
×≈
′
=δλ
×≈
′
=δϕ
°=ϕ
−
−
then
63.41rad7.51020135
2
1
8
′′
−≈××××−≈δλ∆
−
At that latitude, this corresponds to –222 m!
4 THE DIRECT GEODETIC PROBLEM
It can be seen from the above, that the errors of sim-
plifications are neglectable or significant, depending
on the required accuracy and the values of φ, δφ and
δλ, but all of them are systematic and are integrated
in dead reckoning.
These simplifications have been necessary to re-
duce the number of calculations on the ellipsoid and
justified in times of manual mechanical or electronic
calculators, but are completely unnecessary and un-
justified in times of computer calculations. There-
fore we will directly apply the solution of the prob-
lem known in geodesy as direct geodetic problem.
In the solution of the direct geodetic problem
(Fig. 2) from the given coordinates φ
1
, λ
1
and azi-
muth α
1-2
at the start of geodesics P
1
and their length
S are calculated coordinates φ
2
, λ
2
of the endpoint P
2
and the reversed azimuth α
2-1
, on any reference el-
lipsoid.
E. M. Sodano (Sodano 1958, 1965, 1967) from
Helmert’s classical iterative formulae derived a ri-
gorous non-iterative procedure, for any length of ge-
odesics and for any required accuracy, which is at-
tached in Appendix A. This procedure will be used
in this paper in the formal form
φ
2
, λ
2
= SDGP (φ
1
, λ
1
, α
1-2
, S) (22)
α
2-1
= SDGP (φ
1
, λ
1
, α
1-2
, S) (23)
5 APPLICATION OF THE DIRECT GEODETIC
PROBLEM
5.1 Position from range and bearings
We search for the position P(φ
P
, λ
P
)
for which we
have the range d and the bearing α from, or the bear-
ing β to, known position P′(φ
P′
, λ
P′
) (Fig. 1).
The solution is
φ
P
, λ
P
= SDGP (φ
P′
, λ
P′
, α, d) (24)
or
φ
P
, λ
P
= SDGP (φ
P′
, λ
P′
, β – 180°, d) (25)
5.2 Dead reckoned position
We search for the position P(φ
P
, λ
P
) dead reckoned
from known position P′(φ
P′
, λ
P′
) with the speed over
ground V
g
and the course over ground C
g
during the
time interval Δt (Fig. 3).
The solution is
φ
P
, λ
P
= SDGP (φ
P′
, λ
P′
, C
g
, V
g
Δt) (26)
Figure 2. Direct and inverse geodetic problem
Figure 3. Dead reckoned position
N
1
2
1-2
P
2-1
P
S
P'
P
N
V
g
t
C
g
530
5.3 Position from two ranges
We search for the position P(φ
P
, λ
P
) for which we
have two ranges d
1
and d
2
to known positions
P′
1
(φ
P′1
, λ
P′1
) and P′
2
(φ
P′2
, λ
P′2
) (Fig. 4).
The solution is iterative
φ
P1
, λ
P1
= SDGP (φ
P′1
, λ
P′1,
α
1
= var, d
1
) (27)
φ
P2
, λ
P2
= SDGP (φ
P′2
, λ
P′2
, α
2
= var, d
2
) (28)
where α
1
and α
2
are adjusted by any small incre-
ments until e.g.
mincoscos)()(
1P2P
2
1P2P
2
1P2P
=ϕϕλ−λ+ϕ−ϕ
(29)
This iterative process, although looks as very
complicated, is very fast and simple with using e.g.
the Solver in Microsoft Excel.
5.4 Position from two bearings
We search for the position P(φ
P
, λ
P
) for which we
have the bearing α
1
from, or the bearing β
1
to,
known position P′
1
(φ
P′1
, λ
P′1
) and the bearing α
2
from, or the bearing β
2
to, known positions P′
2
(φ
P′2
,
λ
P′2
) (Fig. 4).
The solution is iterative
φ
P1
, λ
P1
= SDGP (φ
P′1
, λ
P′1,
α
1
, d
1
= var) (30)
φ
P2
, λ
P2
= SDGP (φ
P′2
, λ
P′2
, α
2
, d
2
= var) (31)
or
φ
P1
, λ
P1
= SDGP (φ
P′1
, λ
P′1,
β – 180°, d
1
= var) (32)
φ
P2
, λ
P2
= SDGP (φ
P′2
, λ
P′2
, β – 180°, d
2
= var) (33)
or any combination of the above, where d
1
and d
2
are
adjusted by any small increments until e.g. Equation
29 is fulfilled.
5.5 Position from range and bearing to different
positions
We search for position P(φ
P
, λ
P
) for which we have
the range d
1
to known position P′
1
(φ
P′1
, λ
P′1
) and the
bearing α
2
from, or the bearing β
2
to, known posi-
tions P′
2
(φ
P′2
, λ
P′2
) (Fig. 4).
The solution is iterative
φ
P1
, λ
P1
= SDGP (φ
P′1
, λ
P′1,
α
1
= var, d
1
) (34)
φ
P2
, λ
P2
= SDGP (φ
P′2
, λ
P′2
, α
2
, d
2
= var) (35)
or
φ
P2
, λ
P2
= SDGP (φ
P′2
, λ
P′2
, β – 180°, d
2
= var) (36)
where α
1
and d
2
are adjusted by any small incre-
ments until e.g. Equation 29 is fulfilled.
5.6 Position from any combination of ranges and
bearings
The above can be easily extended to any number of
combination of ranges and bearings - we search for
position P(φ
P
, λ
P
) for which we have n ranges d or
bearings α from, or bearings β to, n known positions
P′(φ
P′
, λ
P′
) (Fig. 5).
The solution is iterative
φ
P1
, λ
P1
= SDGP (φ
P′1
, λ
P′1,
α
1
= var, d
1
) (37)
φ
P2
, λ
P2
= SDGP (φ
P′2
, λ
P′2
, α
2
, d
2
= var) (38)
……………………………………….
Figure 4. Positions from two ranges or bearings
Figure 5. Position from n ranges and bearings
d
P'
P
N
P'
d
2
2
1
1
1
2
2
1
d
P'
P
N
P'
d
i
n
1
1
1
n
n
i
1
d
i
n
P'
i
531
φ
Pi
, λ
Pi
= SDGP (φ
P′i
, λ
P′i
, β
i
– 180°, d
i
= var) (39)
……………………………………….
φ
Pn
, λ
Pn
= SDGP (φ
P′n
, λ
P′n
, β
n
– 180°, d
n
= var) (40)
where d or α or β are adjusted by any small incre-
ments until e.g.
∑
=
=ϕλ−λ+ϕ−ϕ
n
1i
Pi
22
Pi
2
Pi
mincos)()(
(41)
where
∑
=
ϕ=ϕ
n
1i
Pi
n
1
(42)
∑
=
λ=λ
n
1i
Pi
n
1
(43)
It is worth mentioning, that we will achieve the
least square errors position in the case of excessive
number of position lines.
5.7 Position lines of different accuracies
In the case of position lines of different accuracies
we can extend Equation 29 or 41 with weights – e.g.
reciprocal of mean square errors m
i
∑
=
=
ϕλ−λ+ϕ−ϕ
n
1i
2
i
Pi
22
Pi
2
Pi
min
m
cos)()(
(44)
5.8 Bearings for long ranges
For long ranges
α ≠ β – 180° (45)
If this difference is significant, for Section 5.1,
we at first iteratively search for α from the equation
α
2-1
= SDGP (φ
P′
, λ
P′
, α = var, d) (46)
until
α
2-1
= β (47)
and then
φ
P
, λ
P
= SDGP (φ
P′
, λ
P′
, α, d) (48)
For Section 5.4 and 5.5 Equations 32, 33, 36, 39
and 40 becomes respectively to
φ
Pi
, λ
Pi
= SDGP (φ
P′i
, λ
P′i
, α
i
= var, d
i
= var) (49)
and additionally
α
2-1i
= SDGP (φ
P′i
, λ
P′i
, α
i
= var, d
i
= var) (50)
as well as Equations 29 and 41 should be supple-
mented by the component e.g.
ϕ
β
+
ϕϕ
β
α−β
−
2
PiM
ii
2
PiPiN
ii
2
i12i
)(R
cosd
cos)(R
sind
)(
(51)
for each β
i
.
6 ACCURACY OF THE SOLUTION OF THE
DIRECT GEODETIC PROBLEM
“The accuracy of geodetic distances computed
through the e
2
, e
4
, e
6
order for very long geodesics is
within a few meters, centimeters and tenth of milli-
meters respectively. Azimuths are good to tenth,
thousandths and hundreds thousandths of a second.
Further improvement of results occurs for shorter
lines” (Sodano 1958).
7 DIRECT COMPUTATION FORM SIMPLIFIED
For shorter distances (the abovementioned “very
long geodesics” means even 20 000 km) or lower
required accuracies we can use equations from Ap-
pendix A reduced to e
2
and f order. Therefore Equa-
tion A 9 becomes to
)cossin(em
4
1
sinea
2
1
SSS
2'
1
S
2'
1S0
ΦΦ+Φ−+
Φ−Φ=Φ
(52)
and Equation A 12 becomes to
γ+βΦ−=
0S
cosfL
(53)
8 CIRCULAR FUNCTIONS
The angles α
2-1
and γ from Equations A 10 and A 11
have to be calculated with the circular function
tan
-1
(), but this function gives solutions in the range
(-90°, 90°). For full range (0°, 360°) retrieving tables
of quadrants are used in Sodano 1965.
For computer calculations a special procedure
should be used to retrieve the full range (0°, 360°)
from the signs of the numerator N and the denomi-
nator D and to detect and support a division by zero
case e.g.:
For
D
N
tanangle
1−
=
IF D ≠ 0 THEN ANGLE = ATN(N/D)
IF D < 0 THEN ANGLE = ANGLE + 180°: END IF
532
ELSE
ANGLE = (2 - SIGN(N))*ABS(SIGN(N))*90°
END IF
IF ANGLE<0 THEN ANGLE=ANGLE+360°: END IF
9 CONCLUSIONS
Presented procedures are quite general and univer-
sal. They can be used for any number of ranges and
bearings, any combinations of ranges and bearings,
any ranges – from meters up to 20 000 km, with al-
most any required accuracy, on any reference ellip-
soid and can calculate the optimal position according
to any objective function.
REFERENCES
Lenart A.S. 1985. Errors of algorithms for position determina-
tion from hyperbolic navigation systems. Marine Geodesy
9(1): 93-111.
Sodano E.M. 1958. A rigorous non-iterative procedure for ra-
pid inverse solution of very long geodesics. Bulletin Géo-
désique 47/48: 13-25.
Sodano E.M. 1965. General non-iterative solution of the in-
verse and direct geodetic problems. Bulletin Géodésique
75: 69-89.
Sodano E.M. 1967. Supplement to inverse solution of long ge-
odesics. Bulletin Géodésique 85: 233-236.
APPENDIX A
Direct computation form (Sodano 1965)
Given: φ
1
, λ
1
, α
1-2
, S
Required: φ
2
, λ
2
, α
1-2
Reference ellipsoid: a
0
, b
0
= semi-major and semi-
minor axes
Flattening
0
0
a
b
1f −=
(A 1)
Second eccentricity squared
2
0
2
0
2
0
2
b
ba
e
−
=
′
(A 2)
11
tan)f1(tan ϕ−=β
(A 3)
2110
sincoscos
−
αβ=β
(A 4)
211
coscosg
−
αβ=
(A 5)
0
S
b
S
=
Φ
(A 6)
)cos1)(sin
2
e
1(m
0
2
1
2
2'
1
β−β+=
(A 7)
)sinsingcos)(sinsin
2
e
1(a
S1S1
2
1
2
2'
1
Φβ+Φββ+=
(A 8)
)cossin
8
5
cos
4
1
sin
8
3
(ema
)cossin
32
5
cos
8
1
cossin
64
13
64
11
(em
cossinea
8
5
)cossin(em
4
1
sinea
2
1
S
2
S
SSS
4'
11
S
3
SS
2
S
SSS
4'
2
1
SS
4'
2
1
SSS
2'
1
S
2'
1S0
ΦΦ−
ΦΦ+Φ+
ΦΦ+ΦΦ−
ΦΦ−Φ+
ΦΦ+
ΦΦ+Φ−+
Φ−Φ=Φ
(A 9)
010
0
12
sinsin
cosg
cos
tan
Φβ−Φ
β
=α
−
(A 10)
210101
210
1
cossinsincoscos
sinsin
tan
−
−
αΦβ−Φβ
αΦ
=γ
(A 11)
10SSS1
2
S1
2
S
cos)]cossin(mf
4
3
sinaf
2
3
f[L
γ+βΦΦ−Φ+
Φ+Φ
−=
(A 12)
L
12
+λ=λ
(A 13)
0012
singcossinsin Φ+Φβ=β
(A 14)
)f1(
tan
tan
2
2
−
β
=ϕ
(A 15)