The Museum of HP Calculators

# Loxodromics for the HP-41

This program is Copyright © 2005 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

## Overview

-A loxodromic curve is a line on the surface of the Earth which cuts all meridians at the same angle µ  ( the azimuth )
-The 2 following programs compute D = the length of these curves and the azimuth µ

1°) Spherical Model of the Earth
2°) Ellipsoidal model of the Earth

-The azimuths are measured positively clockwise from North
-And the longitudes are measured positively Eastwards from the meridian of Greenwich

1°) Spherical Model of the Earth

Formulae:        L2 - L1 =  ( Tan µ ). Ln [ ( Tan ( 45° + b2/2 ) ) / (Tan  ( 45° + b1/2 ) ) ]         where  Li = longitudes ;  bi = latitudes
D =  R.( b2 - b1 )/Cos µ                                                                                    R = mean radius of the Earth = 6371 km

-If  cos µ = 0  we have  D = R.( L2 - L1 ).cos b      where  b = b1 = b2

Data Registers: /
Flags: /
Subroutines: /

-Synthetic registers M & N may be replaced by any standard registers.

01  LBL "LOX"
02  DEG
03  HR
04  STO M
05  STO N
06  X<> T
07  HMS-
08  HR
09  360
10  MOD
11  PI
12  R-D
13  X>Y?
14  CLX
15  ST+ X
16  -
17  X<>Y
18  HR
19  ST- M
20  2
21  ST/ T
22  /
23  45
24  ST+ T
25  +
26  TAN
27  R^
28  TAN
29  X<>Y
30  /
31  LN
32  R-D
33  R-P
34  X<>Y
35  HMS
36  RCL M
37  LASTX
38  COS
39  X#0?
40  GTO 00
41  +
42  X<> N
43  COS
44  R^
45  *
46  1
47  LBL 00
48  /
49  ABS
50  D-R
51  6371
52  *
53  CLA
54  END

( 78 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) / Y L2 ( ° ' " ) µ ( ° ' " ) X b2 ( ° ' " ) D ( km )

Example1:

U.S. Naval Observatory at Washington (D.C.)    L1 = 77°03'56"0 W =  -77°03'56"0   ;   b1 = 38°55'17"2 N =  +38°55'17"2
The Observatoire de Paris:                                  L2 =   2°20'13"8 E  =   +2°20'13"8    ;   b2 = 48°50'11"2 N =  +48°50'11"2

-77.03560  ENTER^
38.55172  ENTER^
2.20138  ENTER^
48.50112  XEQ "LOX"  >>>>   D = 6436.5499 km     X<>Y    µ = 80°08'14"     ( in 4 seconds )

Example2:

0    ENTER^
30   ENTER^
120  ENTER^
30      R/S      >>>>  D = 11555.7158 km    X<>Y    µ = 90°

2°) Ellipsoidal Model of the Earth

-Higher accuracy is obtained if we take the Earth's flattening into account.

-Equatorial radius of the Earth = a = 6378.137 km
-Eccentricity of the ellipsoid = e = ( 0.006694385 )1/2  which corresponds to a flattening  f = 1/298.257    ( e2 = 2.f - f 2 )

Formulae:

L2 - L1 =  ( Tan µ ).{ Ln [ ( Tan ( 45° + b2/2 ) ) / (Tan  ( 45° + b1/2 ) ) ] - (e/2).Ln [ ( 1 + e.sin b2 ).( 1 - e.sin b1 )/( 1 - e.sin b2 )/( 1 + e.sin b1 ) ] }
D =  [ a.( 1 - e2 )/ cos µ ] §b1b2  ( 1 - e2 sin2 t ) -3/2 dt     ( § = Integral )

-One could use any integration formula but "LOX2" evaluates this integral by a series expansion up to the terms in  e6
( the first neglected term is proportional to  e8 )

D ~ [ a.( 1 - e2 )/ cos µ ] [ ( 1 + 3e2/4 + 45e4/64 + 175e6/256 ).( b2 - b1 ) - ( 3e2/8 + 15e4/32 + 525e6/1024 ).( sin 2b2 - sin 2b1 )
+ ( 15e4/256 + 105e6/1024 ).( sin 4b2 - sin 4b1 ) - ( 35e6/3072 ).( sin 6b2 - sin 6b1 ) ]

-However, we cannot use this formula if   cos µ = 0 , in this case ( i-e  if  b1 = b2 ) we have:

D = a.( cos b ).( L2 - L1 ) / ( 1 - e2 sin2 b ) -1/2    the loxodromic line is the parallel of latitude b = b1 = b2

Data Registers:    R00 = e = 0.0066943851/2  ;  R01 = b1  ;  R02 = b2  ;  R03 =  L2 - L1  ;  R04 =  µ   ( in degrees )
Flags: /
Subroutines: /

01  LBL "LOX2"
02  DEG
03  HR
04  STO 02
05  X<> T
06  HMS-
07  HR
08  360
09  MOD
10  PI
11  R-D
12  X>Y?
13  CLX
14  ST+ X
15  -
16  STO 03
17  X<>Y
18  HR
19  STO 01
20  2
21  ST/ T
22  /
23  45
24  ST+ T
25  +
26  TAN
27  R^
28  TAN
29  /
30  LN
31  STO 04
32  1
33  RCL 01
34  SIN
35  .006694385
36  SQRT
37  STO 00
38  *
39  ST+ Y
40  1
41  -
42  /
43  1
44  RCL 00
45  RCL 02
46  SIN
47  *
48  ST- Y
49  1
50  +
51  /
52  *
53  CHS
54  SQRT
55  LN
56  RCL 00
57  *
58  RCL 04
59  -
60  RCL 03
61  X<>Y
62  R-D
63  R-P
64  X<>Y
65  STO 04
66  RCL 01
67  6
68  *
69  SIN
70  RCL 02
71  6
72  *
73  SIN
74  -
75  3
76  *
77  RCL 02
78  4
79  *
80  SIN
81  RCL 01
82  4
83  *
84  SIN
85  -
86  2657
87  *
88  +
89  RCL 02
90  ST+ X
91  SIN
92  RCL 01
93  ST+ X
94  SIN
95  -
96  2531555
97  *
98  -
99  RCL 02
100  RCL 01
101  -
102  D-R
103  1005052504
104  *
105  +
106   E9
107  /
108  1
109  RCL 00
110  X^2
111  -
112  *
113  RCL 04
114  COS
115  X#0?
116  GTO 00
117  RCL 03
118  D-R
119  RCL 01
120  COS
121  *
122  1
123  RCL 01
124  SIN
125  RCL 00
126  *
127  X^2
128  -
129  SQRT
130  LBL 00
131  /
132  RCL 04
133  HMS
134  X<>Y
135  ABS
136  6378.137
137  *
138  END

( 194 bytes / SIZE 005 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) / Y L2 ( ° ' " ) µ ( ° ' " ) X b2 ( ° ' " ) D ( km )

Example1:

-77.03560  ENTER^
38.55172  ENTER^
2.20138  ENTER^
48.50112  XEQ "LOX2"  >>>>   D = 6453.389979 km     X<>Y    µ = 80°10'15"31     ( in 12 seconds )
the exact value is  D = 6453.389986 km

-Compare this value with the length of the corresponding geodesic:  6181.621794 km
-Rhumb-sailing is easier but slower than orthodromy.

Example2:

0    ENTER^
30   ENTER^
120  ENTER^
30      R/S       >>>>  D = 11578.35364 km    X<>Y    µ = 90°

-Due to roundoff errors, there is a loss of significant digits if the 2 latitudes are very close but not equal.
-For instance:  ( 0 ; 30° ) ( 120 ; 30°00'01" )  gives  D = 11578.956 km  instead of  11578.340
( µ = 89°59'59"45  is correctly computed, however )
-Otherwise, the accuracy is of the order of a few centimeters.
-If you don't want to calculate D , simply replace lines 66 to 137 by a single HMS.
-In order to avoid a DATA ERROR  & an OUT OF RANGE if the latitudes are +/- 90°
key in  89°59'59"9999 instead of 90° , the difference is negligible.