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:        L₂ - L₁ =  ( Tan µ ). Ln [ ( Tan ( 45° + b₂/2 ) ) / (Tan  ( 45° + b₁/2 ) ) ]              L_i = longitudes ;  b_i = latitudes
                                 D =  R.( b₂ - b₁ )/Cos µ                                                                              R = mean radius of the Earth = 6371 km

-If  cos µ = 0  we have  D = R.( L₂ - L₁ ).cos b      where  b = b₁ = b₂
 

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.)    L₁ = 77°03'56"0 W =  -77°03'56"0   ;   b₁ = 38°55'17"2 N =  +38°55'17"2
         The Observatoire de Paris:                           L₂ =   2°20'13"8 E  =   +2°20'13"8    ;   b₂ = 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    ( e² = 2.f - f ² )
 

Formulae:

    L₂ - L₁ =  ( Tan µ ).{ Ln [ ( Tan ( 45° + b₂/2 ) ) / (Tan  ( 45° + b₁/2 ) ) ] - (e/2).Ln [ ( 1 + e.sin b₂ ).( 1 - e.sin b₁ )/( 1 - e.sin b₂ )/( 1 + e.sin b₁ ) ] }
            D =  [ a.( 1 - e² )/ cos µ ] §_b1^b2  ( 1 - e² sin² 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  e⁶
  ( the first neglected term is proportional to  e⁸ )

            D = [ a.( 1 - e² )/ cos µ ] [ ( 1 + 3e²/4 + 45e⁴/64 + 175e⁶/256 ).( b₂ - b₁ ) - ( 3e²/8 + 15e⁴/32 + 525e⁶/1024 ).( sin 2b₂ - sin 2b₁ )
                   + ( 15e⁴/256 + 105e⁶/1024 ).( sin 4b₂ - sin 4b₁ ) - ( 35e⁶/3072 ).( sin 6b₂ - sin 6b₁ ) ] + .......................

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

           D = a.( cos b ).( L₂ - L₁ ) / ( 1 - e² sin² b ) ^-1/2    the loxodromic line is the parallel of latitude b = b₁ = b₂
 

Data Registers:    R00 = e = 0.006694385^1/2  ;  R01 = b₁  ;  R02 = b₂  ;  R03 =  L₂ - L₁  ;  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.