(15C) Haversine Navigation

[SlideRule]

An extract from Haversine Navigation 8 pgs.

"Abstract
This is to provide a practical means for cruise planning using … the hand held calculator …
…
Background
…
The spherical trigonometry equations taught left an ambiguity when the initial course was close to due east or west. 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. Today there are dozens of websites that will give you the course and distance between two points on the earth given the Latitude and Longitude of each. Many of these sites report the distance in kilometers which require a definition of the radius of the earth. The haversine method herein described uses a unit radius; thus, a distance as an angle readily resolves to nautical miles; therefore, I am going to limit this to a simple definition of the haversine and practical topics for cruise planning … on the handheld
calculator …
…
2. Calculator
The equations shown earlier are easily implemented into any handheld calculator with memory storage. I have chosen to show an application using the HP 15C calculator that uses RPN … You only have to program your device once. The RPN programming steps are shown …"

BEST!
SlideRule

[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
  

[C.Ret]

Hi there,

Once again, my thanks for this very well documented topic. Note that this is always the case. It is therefore always with great pleasure that I scrutinize all this with the keenest interest. Opportunities to study and learn so effectively are unfortunately increasingly rare.

If I understand all this correctly, the final multiplication by 60 gives the displacement in nautical miles.

So I will also check the following assertion:

(11-05-2022 11:13 PM)Thomas Klemm Wrote:  But similar programs work for other calculators as well as longs as they provide →REC, →POL, R↑ and R↓.

So I'm going to try to translate these clever codes for a SHARP EL-5150 calculator which also has the instructions →REC, →POL, ↑ and ↓ on its keyboard.

   

This is a good thing done, not without effort, hesitation and perseverance. But note that the instructions ↑ and ↓ are not used in the algebraic expressions of the SHARP EL-5150 but are used to select the correct slot in the AER formula pool.

Thanks.

[Thomas Klemm]

(11-06-2022 05:14 PM)C.Ret Wrote:  If I understand all this correctly, the final multiplication by 60 gives the displacement in nautical miles.

This is correct.

From Nautical mile:

Quote:Historically, it was defined as the meridian arc length corresponding to one minute (\(\frac{1}{60}\) of a degree) of latitude.

While not exactly 1852m it is probably close enough.



Even though I've never used this calculator the program is quite readable.
Both →REC and →POL appear to use infix notation which is rather unusual.

Here's what I came up with for the HP-48:

Code:

\<< \-> \Gh1 \Gl1 \Gh2 \Gl2
    \<< 1 \Gl2 \Gl1 - 0 ACOS \Gh2 -
        SPHERE \->V3 RECT V\->
        \-> u v w
        \<< w  u \->V2
            SPHERE V\-> \Gh1 + \->V2
            RECT V\-> v SWAP \->V3
            SPHERE V\-> RECT
        \>>
    \>> ROT DROP 60 *
\>>

Since the latitude is measured from the north-pole we have to calculate the complement \(c = 90^{\circ} - \theta_2\).
And then the transformation between spherical and rectangular coordinates seems a bit awkward.
Thus I'm not super happy with it.