253
1 INTRODUCTION
The traditional method of computing the shortest
pathbetweentwopointsontheEarthknownasGreat
CirclemethodapproximatetheEarthasthesphere.
These simplifications have been necessary to re‐
ducethenumberofcalculationsandjustifiedintimes
of manual mechanical or electronic calculators, but
are completely unnecessary and unjustified in times
of computer calculations. Therefore we will directly
applythesolutionoftheproblemknowningeodesy
astheinversegeodeticproblem.
In the solution of the inverse geodetic problem
(Fig.1)fromthegivencoordinatesφ
1,λ1 atthestart
ofgeodesicP
1andcoordinatesφ2,λ2oftheendpoint
P
2arecalculatedthelengthS,theazimuthα1‐2andthe
reversedazimuthα
2‐1,onanyreferenceellipsoid.
E. M. Sodano (Sodano 1958, 1965, 1967) from
Helmert’s classical iterative formulas derived a
rigorous non‐iterative procedure, for any length of
geodesic and for any required accuracy, which is
attached (as exemplary) in Appendix A. This
procedure (or any other solution of the inverse
geodetic
problem) will be used in this paper in the
formalnotation
S=IGP(φ
1,λ1,φ2,λ2) (1)
α
1‐2,α2‐1=IGP(φ1,λ1,φ2,λ2) (2)
N
1
2
1-2
P
2-1
P
S
Figure1.Thedirectandtheinversegeodeticproblem.
Solutions of Inverse Geodetic Problem in Navigational
Applications
A.S.Lenart
GdyniaMaritimeUniversity,Gdynia,Poland
ABSTRACT: Solutions of such navigational problems as an orthodromic navigation (courses, distances and
intermediate points), maximum latitude and a composite navigation with limited latitude as well as, for
comparison, a loxodromic navigation (courses, distances) without any simplifications for a sphere, by an
application of solutions
of the inverse geodetic problem are presented. An exemplary rigorous, rapid, non‐
iterativesolutionoftheinversegeodeticproblemaccordingtoSodano,foranylengthofgeodesic,isattached.
http://www.transnav.eu
the International Journal
on Marine Navigation
and Safety of Sea Transportation
Volume 7
Number 2
June 2013
DOI:10.12716/1001.07.02.13
254
2 ORTHODROMICDISTANCEANDCOURSES
Wedefinetheorthodromeastheshortestpathonany
surfaceandnotonlytheGreatCircledistanceonthe
sphereascommonlyisused.
Thegeodesicis(locally‐notlongwayround)the
shortest path between two points on an ellipsoid of
revolution. Therefore
we can obtain orthodromic
distanceandcoursesdirectlyfromEquations1and2
withnavigationalsubstitutions
C
gs=α1‐2 (3)
C
ge=α2‐1‐180° (4)
whereC
gs=thecourseovergroundatthestartofthe
orthodrome;andC
ge=thecourseovergroundatthe
endoftheorthodrome.
Eastlongitudesandnorthlatitudesareconsidered
positiveandwestlongitudesandsouthlatitudesare
considerednegative.
3 ACCURACYOFTHESOLUTIONOFTHE
INVERSEGEODETICPROBLEM
“The accuracy of geodetic distances computed
throughthee
2
,e
4
,e
6
orderforverylonggeodesicsis
within a few meters, centimeters and tenth of milli‐
meters respectively. Azimuths are good to tenth,
thousandthsandhundredsthousandths of a second.
Further improvement of results occurs for shorter
lines”(Sodano1958).
This accuracy can be easy tested in the case of
equatorial
orthodrome. Substitutionφ1 =φ2 = 0 to
EquationsA2toA10yields
L)ff1(bS
2
0
(5)
whereasthecorrectvalueisgivenbytheequation
L
f1
b
LaS
0
0
*
(6)
thereforetherelativeerroris
*
9
*
38 10 38cm /10 000 km
SS
S
S

(7)
4 ERRORSOFCALCULATIONSONTHESPHERE
According to Euler’s theorem for an ellipsoid of
revolutiontheradiusofcurvatureinmeridianisthe
smallest and the radius of curvature in the prime
verticalisthelargestatapoint.Theseradiiaregiven
respectivelybytheequations
322
2
0
M
)sine1(
)e1(a
R
(8)
22
0
N
sine1
a
R
(9)
The widest span has the radius of curvature in
meridiansince
0
2
0
MminM
a
b
)0(RR
(10)
0
2
0
NMmaxM
b
a
)°90(R)°90(RR
(11)
The substitute radius of curvature of any
orthodromewillbewithintheselimits.Theminimum
absolutevalueofdeviationgivesanassumptionthat
the substitute radius of sphere is given by the
equation(foraglobalrangeoflatitudes)
2
RR
R
minMmaxM
S
(12)
Thenthemaximal relative errorof calculation on
such a sphere, instead of an ellipsoid, gives the
equation
km00010/km50%5.0
RR
)RR(
S
minMmaxM
minMmaxM
S
(13)
These results are similar to obtained by Earle
(2006)withmuchmorecomplicatedmethods.
5 INTERMEDIATEPOINTSONTHE
ORTHODROME
For calculating intermediate points on the
orthodromewecanuse,asexemplary,thesolutionof
the direct geodetic problem presented in Lenart
(2011), also according to Sodano, having similar
accuracy.
In the solution of the direct geodetic problem
(Fig.1)fromthegivencoordinatesφ
1,λ1andazimuth
α
1‐2atthe start of geodesic P1 and their length S are
calculated coordinatesφ
2,λ2 of the endpoint P2 and
thereversedazimuthα
2‐1,onanyreferenceellipsoid.
Thisprocedure(oranyothersolutionofthedirect
geodetic problem) will be used in this paper in the
formalnotation
φ
2,λ2 =DGP(φ1,λ1,α1‐2,S) (14)
α
2‐1=DGP(φ1,λ1,α1‐2,S) (15)
255
TheorthodromeofthelengthSwewilldividefor
nsuborthodromes(Fig.2)ofanylengthS
isuchas
n
1i
i
SS (fromEquation1) (16)
andintermediatepointsarecalculatedinniterations:
FORi=1ton
IFi=1THEN
(φ
2,λ2)i‐1 = φ1,λ1 (17)
(
2‐1)i‐1‐180°=2‐1(fromEquation2) (18)
ENDIF
(φ
2,λ2)i =DGP((φ2,λ2)i‐1,(α2‐1)i‐1‐180°,Si) (19)
(α
2‐1)i=DGP((φ2,λ2)i‐1,(α2‐1)i‐1‐180°,Si) (20)
NEXTi
Figure2.Intermediatepointsontheorthodrome.
orevenn‐1iterationsbecausethelastiterationisfor
verificationonlythat
(φ
2,λ2)n =P2(φ2,λ2) (21)
In traditional navigation intermediate points are
calculated during planning the voyage to navigate
between them along a loxodrome. If we have
programmed procedure on the bridge during the
voyage then the situation can be quite different. In
thiscasetheintermediatepointsareneedede.g.for
theverification
ofthepathonthemaponly.During
thevoyageifweenterasP
1thecurrentpositionand
P
2isconstantlytheendpointthentheproceduregives
the course over ground for the orthodrome to the
endpoint,evenifwee.g.duesetanddriftareoutof
track, and not to an intermediate point. We can
calculateanewcoursefortheorthodromeaftereach
position fixing
and to navigate always along the
current orthodrome without these calculated
intermediatepoints.
6 MAXIMUMLATITUDE
According to Clairaut’s relation for a geodesic on a
surfaceofrevolution
rsin
=const.=C (22)
wherer=theradiusofparallel.
Foranellipsoidofrevolution
r=R
Ncosφ (23)
Atφ
max
1sin
(24)
therefore
C
sine1
cosa
max
2
max0
(25)
andfinally
222
0
22
0
max
Cea
Ca
sin
(26)
whereCe.g.is
C=r(φ
1)sinα1‐2 (27)
Theaboveequationsarevalidif
1sin
(28)
existsonourorthodromei.ewhen
(C
gs‐90°)(Cge‐90°)<0 (29)
or
(C
gs‐270°)(Cge‐270°)<0 (30)
orelse
),max(
21max
(31)
7 COMPOSITENAVIGATION
Ifforanyreasonφ
maxislimitedtoφlimthen(Fig.3)we
have:
1 The orthodome I (Ort I) – fromφ1,λ1 toφlim,
λ2OI.
2 Theloxodrome(Lx)–atφlimfromλ2OItoλ1OII.
3 TheorthodomeII(OrtII)–fromφlim,λ1OIItoφ2,
λ2.
We will obtainλ
2OI andλ1OII in iterative
procedures:
256
C
ge=IGP(φ1,λ1,φlim,λ2OI=var) (32)
whereλ
2OIisadjustedbyanysmallincrements until
C
ge=90°ifL>0or270°ifL<0
Figure3.Compositenavigation.
and
C
gs=IGP(φlim,λ1OII=var,φ2,λ2) (33)
whereλ
1OIIisadjustedbyanysmallincrementsuntil
C
gs=90°ifL>0or270°ifL<0
This iterative process, although looks as very
complicated, is very fast and simple with using e.g.
theSolverinMicrosoftExcel.
8 LOXODROMICDISTANCEANDCOURSE
The ortodromic navigation is for shorter distances
then in the loxodromic navigation therefore
to
calculate this difference we will calculate, for
comparison, the loxodromic distance and course on
anellipsoid.
lx
21M
lx
cos
),(S
S
(34)
whereS
M(φ1,φ2)=themeridiandistancebetweenφ1
andφ
2; and lx = the course over ground for
loxodrome.
L
tan
lx
(35)
where
)
sine1
sine1
ln
sine1
sine1
(ln
2
e
))2/4/(tgln())2/4/(tgln(
1
1
2
2
12
(36)
Inourcase
S
M(φ1,φ2)=IGP(φ1,λ1,φ2,λ1) (37)
Equation33isvalidifφ
1≠φ2orelse
L)(r = S
21lx
(38)
9 INVERSECOMPUTATIONFORMSIMPLIFIED
For shorter distances (the very long geodesic in
paragraph3meanseven20000km)orlowerrequired
accuracies we can use equations from Appendix A
reduced to f order (having the accurate solution for
referenceinerrorscalculations).ThereforeEquationA
10becomes
to
]2/sin)cosma2)(ff(
2/m)ff()ff1[(bS
2
22
0
(39)
andEquationA11becomesto
Lc)ff(
2
(40)
or
]sin)cosma2(c[BAS
2
SS
(41)
where
0
2
S
b)]ff(
2
1
1[A
(42)
0
2
S
b)]ff(
2
1
[B
(43)
It is evident that f + f
2
are from reduced higher
orderelementsfromseries
1
f1
1
fff
32
(44)
Notingthat
0
0
a
f1
b
(45)
Equations42and43arethus
2/)ba(A
00S
(46)
2/)ba(B
00S
(47)
This simplified computation form gives errors in
therangeofmeters(andhasnoerrorsforequatorial
orthodromes.
10 CIRCULARFUNCTIONS
Theanglesα
1‐2,α2‐1fromEquationsA12,A13andlx
from Equation 35 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).
257
For computer calculations a special procedure
should be used to retrieve the full range (0°, 360°)
fromthe signsof the numerator Nand the denomi‐
natorDandtodetectandsupportadivisionbyzero
casee.g.:
For
D
N
TANANGLE
1
IFD≠0THEN
ANGLE=ATN(N/D)
IFD<0THENANGLE=ANGLE+180°:ENDIF
ELSE
ANGLE=(2‐SGN(N))*ABS(SGN(N))*90°
ENDIF
IFANGLE<0THENANGLE=ANGLE+360°:ENDIF
11 CONCLUSIONS
Thesetofpresentedproceduresarequitegeneraland
universal.Theycanbeusedwithanysolutionsofthe
inverse(anddirect)geodeticproblemsas
wellduring
the voyage planning as during the voyage in real
time,for“full”orthodromenavigation.
REFERENCES
Earle M.A. 2006. Sphere to Spheroid Comparisons. The
JournalofNavigation59:491‐496.
Lenart A.S. 2011. Solutions of Direct Geodetic Problem in
Navigational Applications. Transnav – International
Journal on Marine Navigation and Safety of Sea
Transportation5(4):527‐532.
Sodano E.M. 1958. A rigorous non‐iterative procedure for
rapid inverse
solution of very long geodesics. Bulletin
Géodésique47/48:13‐25.
SodanoE.M.1965.Generalnon‐iterativesolutionofthein‐
verse and direct geodetic problems. Bulletin Géodésique
75:69‐89.
Sodano E.M.1967.Supplementtoinverse solution of long
geodesics.BulletinGéodésique85:233‐236.
APPENDIXA
Inversecomputationform(Sodano1965,1967)
Given:φ
1,λ1,φ2,λ2
Required:α
1‐2,α2‐1,S
‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
Reference ellipsoid: a
0, b0 = semi‐major and semi‐
minoraxes
Flattening
0
0
a
b
1f
(A1)
‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
11
tan)f1(tan
(A2)
22
tan)f1(tan
(A3)
21
sinsina
(A4)
21
coscosb
(A5)
12
L
(A6)
Lcosbacos
(A7)
sin
Lsinb
c
(A8)
2
c1m
(A9)
]
2
csc)cosma)(m1(f
8
cossin)cosma2(f
16
)cossin(mf
2
sin)cosma2)(ff(
2
m)ff(
)ff1[(bS
22
22
22
2
2
2
0
(A10)
Lcm]cotf
4
cossinf
4
f5
[
a]cscf
2
sinf
[)ff(
22
22
22
2
2
(A11)
coscossincossin
sincos
tan
2112
2
21
(A12)
)coscossincos(sin
sincos
tan
1221
1
12
(A13)