loxodromics

# Loxodromics for the HP-41

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 ) ) ]              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.