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.
 
 

 
  ( 78 bytes / SIZE 000 )
 
 

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: /
 
 

 
   ( 194 bytes / SIZE 005 )
 
 

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.