(15C) Haversine Navigation

[Thomas Klemm]

Motivation

Quote:Further solutions using the Spherical Law of Cosines fall apart when the distances are very small; short distances are most important to the vast majority of weekend sailors. Haversines resolved these issues.

Haversine Function

\(
\begin{align}
\operatorname {hav} (\theta) = \sin^2\left(\frac{\theta}{2}\right) = \frac{1 - \cos(\theta)}{2}
\end{align}
\)

To avoid cancellation for small angles only the expression using the sine function should be used.

Haversine Formula

This formula can be used to calculate the great-circle distance between two points on a sphere:

\(
\begin{align}
\operatorname {hav} \left(\theta \right) = \operatorname {hav} \left(\varphi_2 - \varphi_1\right) + \cos \left(\varphi_1\right)\cos \left(\varphi_2\right)\operatorname {hav} \left(\lambda_2-\lambda_1\right)
\end{align}
\)

We can see that this avoids cancellation since only \( \operatorname {hav} \) is used for small differences.
This happens if the two points are close which appears to be common in sailing.

Compare this to the common formula:

\(
\begin{align}
\cos(\theta) = \sin(\varphi_1) \sin(\varphi_2) + \cos(\varphi_1) \cos(\varphi_2) \cos(\lambda_{2}-\lambda_{1})
\end{align}
\)

Here cancellation happens if \(\cos(\Delta \lambda) \) is calculated for small \( \Delta \lambda = \lambda_2-\lambda_1 \).
A similar problem arises if \( \theta \) is small.

Program

This is the original program:

Code:

001 - 42,21,11  LBL A
002 -    45  5  RCL 5
003 -    32  2  GSB 2
004 -    45  8  RCL 8
005 -    32  2  GSB 2
006 -       30  -
007 -    32  0  GSB 0
008 -    45  5  RCL 5
009 -    32  2  GSB 2
010 -       23  SIN
011 -    45  8  RCL 8
012 -    32  2  GSB 2
013 -       23  SIN
014 -       20  x
015 -    45  6  RCL 6
016 -    45  9  RCL 9
017 -       30  -
018 - 43, 5, 0  CF 0
019 - 43,30, 2  TEST 2
020 - 43, 4, 0  SF 0
021 -    32  0  GSB 0
022 -       20  x
023 -       40  +
024 -    32  1  GSB 1
025 -    44  3  STO 3
026 -    45  8  RCL 8
027 -    32  2  GSB 2
028 -    32  0  GSB 0
029 -    45  5  RCL 5
030 -    32  2  GSB 2
031 -    45  3  RCL 3
032 -       30  -
033 -    32  0  GSB 0
034 -       30  -
035 -    45  5  RCL 5
036 -    32  2  GSB 2
037 -       23  SIN
038 -    45  3  RCL 3
039 -       23  SIN
040 -       20  x
041 -       10  /
042 -    32  1  GSB 1
043 - 43, 6, 0  F? 0
044 -    22  3  GTO 3
045 -        3  3
046 -        6  6
047 -        0  0
048 -       34  x<>y
049 -       30  -
050 - 42,21, 3  LBL 3
051 -    44  2  STO 2
052 -        6  6
053 -        0  0
054 - 44,20, 3  STOx 3
055 -       34  x<>y
056 -    45  3  RCL 3
057 -       34  x<>y
058 -    43 32  RTN
059 - 42,21, 0  LBL 0
060 -       24  COS
061 -        1  1
062 -       34  x<>y
063 -       30  -
064 -        2  2
065 -       10  /
066 -    43 32  RTN
067 - 42,21, 1  LBL 1
068 -        2  2
069 -       20  x
070 -        1  1
071 -       34  x<>y
072 -       30  -
073 -    43 24  ACOS
074 -    43 32  RTN
075 - 42,21, 2  LBL 2
076 -       23  SIN
077 -    43 24  ACOS
078 -    43 32  RTN
079 - 42,21,12  LBL B
080 -    45  2  RCL 2
081 -    45  0  RCL 0
082 -       40  +
083 -    44  1  STO 1
084 -        3  3
085 -        6  6
086 -        0  0
087 - 43,30, 7  TEST 7
088 -    22  4  GTO 4
089 -       30  -
090 -    44  1  STO 1
091 - 42,21, 4  LBL 4
092 -    45  1  RCL 1
093 -    43 32  RTN

Assume that the coordinates are stored in the following registers:

R5: \(φ_1\)
R6: \(λ_1\)
R8: \(φ_2\)
R9: \(λ_2\)

Let's start with the first few lines:

Code:

002 -    45  5  RCL 5
003 -    32  2  GSB 2
004 -    45  8  RCL 8
005 -    32  2  GSB 2
006 -       30  -

The subroutine 2 calculates the complement of the angle:

Code:

075 - 42,21, 2  LBL 2
076 -       23  SIN
077 -    43 24  ACOS
078 -    43 32  RTN

But instead of calculating \((90^{\circ} - \varphi_1) - (90^{\circ} - \varphi_2)\) we could calculate \(\varphi_2 - \varphi_1\) directly.

Quote:Complementing an angle takes less programming steps with the Sin aCos method than subtracting the angle from 90°.

It is used again when calculating \(\cos \left(\varphi_1\right)\cos \left(\varphi_2\right)\):

Code:

008 -    45  5  RCL 5
009 -    32  2  GSB 2
010 -       23  SIN
011 -    45  8  RCL 8
012 -    32  2  GSB 2
013 -       23  SIN
014 -       20  x

But here we could just as well use COS instead of SIN.

Then I had a look at the implementation of the haversine function:

Code:

059 - 42,21, 0  LBL 0
060 -       24  COS
061 -        1  1
062 -       34  x<>y
063 -       30  -
064 -        2  2
065 -       10  /
066 -    43 32  RTN

Unfortunately the variant where cancellation for small angles occur was chosen.
This invalidates the reason to use the haversine function.

So I came up with this improved version:

Code:

001 - 42,21,13  LBL C       
002 -    45  5  RCL 5       
003 - 45,30, 8  RCL- 8      
004 -    32  8  GSB 8       
005 -    45  5  RCL 5       
006 -       24  COS         
007 -    45  8  RCL 8       
008 -       24  COS         
009 -       20  x           
010 -    45  6  RCL 6       
011 - 45,30, 9  RCL- 9      
012 -    32  8  GSB 8       
013 -       20  x           
014 -       40  +           
015 -    32  9  GSB 9       
016 -    44  7  STO 7       
017 -    43 32  RTN         
018 - 42,21, 8  LBL 8       
019 -        2  2           
020 -       10  /           
021 -       23  SIN         
022 -    43 11  x^2         
023 -    43 32  RTN         
024 - 42,21, 9  LBL 9       
025 -       11  SQRT        
026 -    43 23  ASIN        
027 -        2  2           
028 -       20  x           
029 -    43 32  RTN

However it calculates only the distance, not the bearing.

When I analysed the accuracy of the different approaches I noticed that with 10 digits it doesn't really matter.
The coordinates are a couple of minutes appart which translates to a few nautical miles.

Here's a typical example:

\(φ_1\) = 41.6258
\(φ_2\) = 41.6683
\(λ_1\) = -71.9950
\(λ_2\) = -71.8650

However we'd notice the difference if either the differences were smaller or then we would use a calculator with fewer digits.

Can we do better?

Great Circle Navigation



To match the drawing with our problem consider the following:

\(
\begin{align}
A &= \lambda_2-\lambda_1 \, \text{(the angle at w)} \\
b &= 90^{\circ} - φ_1 \\
c &= 90^{\circ} - φ_2 \\
\end{align}
\)

Here \(a\) is the great-circle distance while \(C\) is the bearing.

Spherical Law of Sines

\(
\begin{align}
\frac{\sin A}{\sin a} = \frac{\sin B}{\sin b} = \frac{\sin C}{\sin c}
\end{align}
\)

Spherical Law of Cosines

\(
\begin{align}
\cos c=\cos a\cos b+\sin a\sin b\cos C
\end{align}
\)

With a little rearrangement we get:

\(
\begin{align}
\cos C&= \frac{\cos c-\cos a\cos b}{\sin a\sin b} \\
\end{align}
\)

This leads to:

\(
\begin{align}
\sin C &= \sin c \frac{\sin A}{\sin a} \\
\\
\cos C &= \frac{\cos c - \cos a \cos b}{\sin a \sin b} \\
\\
\cos a & = \cos b \cos c + \sin b \sin c \cos A \\
\end{align}
\)

Or then by multiplying the first two equations by \( \sin a \):

\(
\begin{align}
\sin a \sin C &= \sin c \sin A \\
\\
\sin a \cos C &= \frac{\cos c - \cos a \cos b}{ \sin b} \\
\\
\cos a & = \cos b \cos c + \sin b \sin c \cos A \\
\end{align}
\)

The goal is to use a rectangular to polar transformation to calculate \(a\) and \(C\) from:

\(
\begin{align}
y &= \sin a \sin C \\
x &= \sin a \cos C \\
z &= \cos a \\
\end{align}
\)

For this we have to substitute \( \cos a \) in the 2nd equation:

\(
\begin{align}
\sin a \cos C
&= \frac{\cos c - \cos a \cos b}{ \sin b} \\
\\
&= \frac{\cos c - (\cos b \cos c + \sin b \sin c \cos A) \cos b}{ \sin b} \\
\\
&= \frac{\cos c - \cos^2 b \cos c - \sin b \cos b \sin c \cos A}{ \sin b} \\
\\
&= \frac{\sin^2 b \cos c - \sin b \cos b \sin c \cos A}{ \sin b} \\
\\
&= \sin b \cos c - \cos b \sin c \cos A \\
\end{align}
\)

We can define the following polar to rectangular transformation:

\(
\begin{align}
u &= \sin c \cos A \\
v &= \sin c \sin A \\
w &= \cos c
\end{align}
\)

Pluging all this into the equations leads to:

\(
\begin{align}
y &= \sin a \sin C &= \sin c \sin A & = v \\
x &= \sin a \cos C &= \sin b \cos c - \cos b \sin c \cos A &= w \sin b - u \cos b \\
z &= \cos a &= \cos b \cos c + \sin b \sin c \cos A &= w \cos b + u \sin b \\
\end{align}
\)

This turns out to be a rotation by the angle \(φ_1 = 90^{\circ} - b\):

\(
\begin{align}
\begin{bmatrix}
x \\
z \\
\end{bmatrix}
&=
\begin{bmatrix}
\sin b & - \cos b \\
\cos b & \sin b \\
\end{bmatrix}
\cdot
\begin{bmatrix}
w \\
u \\
\end{bmatrix}
\end{align}
\)

Programs

Here's the program for the HP-42S:

Code:

00 { 37-Byte Prgm }   #
01▸LBL "NAV"          # λ_2 φ_2 λ_1 φ_1
02 X<>Y               # φ_2 λ_2 λ_1 φ_1
03 R↓                 # λ_2 λ_1 φ_1 φ_2
04 X<>Y               # λ_1 λ_2 φ_1 φ_2
05 -                  # A=λ_2-λ_1 φ_1 φ_2 φ_2
06 X<>Y               # φ_1 A φ_2 φ_2
07 R↑                 # φ_2 φ_1 A φ_2
08 1                  # 1 φ_2 φ_1 A
09 →REC               # sin(c) cos(c) φ_1 A
10 R↑                 # A sin(c) cos(c) φ_1
11 X<>Y               # sin(c) A cos(c) φ_1
12 →REC               # u=sin(c)cos(A) v=sin(c)sin(A) w=cos(c) φ_1
13 X<>Y               # v u w φ_1
14 R↓                 # u w φ_1 v
15 X<>Y               # w u φ_1 v
16 →POL               # r φ φ_1 v
17 R↓                 # φ φ_1 v r
18 +                  # φ+φ_1 v r r
19 R↑                 # r φ+φ_1 v r
20 →REC               # x z y r
21 X<>Y               # z x y r
22 R↓                 # x y r z
23 →POL               # r C r z
24 R↑                 # z r C
25 →POL               # 1 a C
26 R↓                 # a C
27 60                 # 60 a C
28 ×                  # 60a C
29 END                #

This is the same program for the HP-15C:

Code:

001 - 42,21,11  LBL A      
002 -       34  x<>y       
003 -       33  Rv         
004 -       34  x<>y       
005 -       30  -          
006 -       34  x<>y       
007 -    43 33  R^         
008 -        1  1          
009 -    42  1  ->R        
010 -    43 33  R^         
011 -       34  x<>y       
012 -    42  1  ->R        
013 -       34  x<>y       
014 -       33  Rv         
015 -       34  x<>y       
016 -    43  1  ->P        
017 -       33  Rv         
018 -       40  +          
019 -    43 33  R^         
020 -    42  1  ->R        
021 -       34  x<>y       
022 -       33  Rv         
023 -    43  1  ->P        
024 -    43 33  R^         
025 -    43  1  ->P        
026 -       33  Rv         
027 -        6  6          
028 -        0  0          
029 -       20  x          
030 -    43 32  RTN

But similar programs work for other calculators as well as longs as they provide →REC, →POL, R↑ and R↓.

Example

41.6258
ENTER
-71.9950
ENTER
41.6683
ENTER
-71.8650
GSB A

6.361974738
x<>y
66.32743461

Pros

• It doesn't use registers.
• No ambiguity when the initial course is close to due east or west.
• No cancellation with small angles.
• The program is smaller and faster.
  

Cons

• The bearing can be negative. If that bothers you add 360.
  

Cancellation

The cancellation is avoided since with →REC both cosine and sine are calculated.
This allows to keep track of small angles.

You can test this with a small angle, say 1.23456789e-5.
Compare the result when calculating COS ACOS with 1 →REC →POL R↓.

Thank you

As always, I'm very grateful for your findings of these programs in the wild.

References

• Haversine Navigation
• Haversine Navigation Excel Workbook
• Haversine formula
• The spherical law of sines
• Spherical law of cosines
• Great-circle distance
• Calculate distance, bearing and more between Latitude/Longitude points