Geodesics & Loxodromes on a Triaxial Ellipsoid for the HP-41

Overview
 

 0°)  2 M-Code Routines
 1°)  Geodesic Distance
 2°)  Loxodromes
 
 

0°)  2 M-Code Routines
 

-In the following programs, the inegrals are evaluated with Gauss-Legendre 7-point formula
-GX7 & GW7 return the abscissas & weights ( the 1st abscissa is zero, so I didn't code it ):

   1  XEQ "GX7"  ->  0.4058451514
   2  XEQ "GX7"  ->  0.7415311856
   3  XEQ "GX7"  ->  0.9491079123

   0  XEQ "GW7"  ->  0.4179591837
   1  XEQ "GW7"  ->  0.3818300505
   2  XEQ "GW7"  ->  0.2797053915
   3  XEQ "GW7"  ->  0.1294849662
 

0B7  "7"
018  "X"
007  "G"
2A0  SETDEC
0F8  C=X
2FC  RCR 13
27E  C=C-1 MS
10E  A=C ALL
04E  C=0 ALL
35C  PT=12
1BE  A=A-1 MS
063  JNC+12d
110
010
150
210
110  C=4.058451514
150
050
150
050
110
0C3  JNC+24d
1BE  A=A-1 MS
063  JNC+12d
1D0
110
050
150
0D0  C=7.415311856
050
050
210
150
190
05B  JNC+11d
250
110
250
050
010  C=9.491079123
1D0
250
050
090
0D0
266  C=C/10
0E8  X=C
3E0  RTN
0B7  "7"
017  "W"
007  "G"
2A0  SETDEC
0F8  C=X
2FC  RCR 13
10E  A=C ALL
04E  C=0 ALL
35C  PT=12
1BE  A=A-1 MS
063  JNC+12d
110
050
1D0
250
150  C=4.179591837
250
050
210
0D0
1D0
12B  JNC+37d
1BE  A=A-1 MS
063  JNC+12d
0D0
210
050
210
0D0  C=3.818300505
010
010
150
010
150
0C3  JNC+24d
1BE  A=A-1 MS
063  JNC+12d
090
1D0
250
1D0
010  C=2.797053915
150
0D0
250
050
150
05B  JNC+11d
050
090
250
110
210  C=1.294849662
110
250
190
190
090
266  C=C/10
0E8  X=C
3E0  RTN

Note:

-The content of X-register is overwritten.
 

1°)  Geodesic Distance
 

>>>  This program is just a slight modification of the last programs listed in "Geodesics on a Triaxial Ellipsoid for the HP41"

 "DGT3" uses an approximate method to evaluate the geodesic distance between 2 points M & N on the Earth
  The error is about a few centimeters for long geodesics, except for nearly antipodal points ( do not use this program in this case )
 

-Let a point P ( on the ellipsoid ) with ( OM , OP' ) = µ , ( OP' , OP ) = h    ,    P' = orthogonal projection of P on the plane (OMN)

    OP = x u + y v + z w       where  u = OM / || OM ||   ,  v = ON / || ON ||  ,  w = u x v          ( x = cross-product )

-We have, with  R = OP  and  W = ( OM , ON )

    x = R Cos h Sin ( W - µ ) / Sin W
    y = R Cos h Sin µ / Sin W
    z = R Sin h / Sin W

-Writing that P(X,Y,Z)  is on the ellipsoid  X²/a² + Y²/b² + Z²/c² = 1,  we find that

  R / Sin W = 1 / ( A² + B² + C² )^1/2   with

      A = [ x₁ Cos h Sin ( W - µ ) + x₂ Cos h Sin µ + ( y₁z₂ - y₂z₁ ) Sin h ] / a
      B = [ y₁ Cos h Sin ( W - µ ) + y₂ Cos h Sin µ + ( z₁x₂ - z₂x₁ ) Sin h ] / b                       where  u ( x₁ , y₁ , z₁ )  &  v ( x₂ , y₂ , z₂ )
      C = [ z₁ Cos h Sin ( W - µ ) + z₂ Cos h Sin µ + ( x₁y₂ - x₂y₁ ) Sin h ] / c
 

-We search an approximation of the function  h(µ) under the form:

   h(µ) = h₁ Sin ( 180° µ / W ) + h₂ Sin ( 360° µ / W ) + h₃ Sin ( 540° µ / W )                 ( the first terms of a Fourier series )

-And the corresponding integral  s = §₀^W [ R² Cos² h  + R² ( dh/dµ )² + ( dR/dµ )²  ] ^1/2 dµ

      =  ( Sin W )    §₀^W [ ( Cos² h ) / T + (1/T) ( dh/dµ )² + ( dr/dµ )²  ] ^1/2 dµ    is evaluated by a Gaussian quadrature ( Gauss-Legendre 7-point formula )

  where   T = ( A² + B² + C² )

    r = 1/sqrt(T)                  dr/dµ = ¶r/¶µ + ( ¶r/¶h ) ( dh/dµ )

    r = 1 / ( A² + B² + C² )^1/2

  the partial derivatives

  ¶r/¶µ = ( -1/sqrt(T³) ) [ (A/a) { -x₁ Cos h Cos ( W - µ ) + x₂ Cos h Cos µ }
                                   + (B/b) { -y₁ Cos h Cos ( W - µ ) + y₂ Cos h Cos µ }
                                   + (C/c) { -z₁ Cos h Cos ( W - µ ) + z₂ Cos h Cos µ } ]
  and

  ¶r/¶h = ( -1/sqrt(T³) ) [ (A/a) { -x₁ Sin h Sin ( W - µ ) - x₂ Sin h Sin µ + ( y₁z₂ - y₂z₁ ) Cos h }
                                    + (B/b) { -y₁ Sin h Sin ( W - µ ) - y₂ Sin h Sin µ + ( z₁x₂ - z₂x₁ ) Cos h }
                                    + (C/c) { -z₁ Sin h Sin ( W - µ ) - z₂ Sin h Sin µ + ( x₁y₂ - x₂y₁ ) Cos h } ]
 

-With  H = W / 1000  , we calculate the following lengths

    A = f(0,0,0)
    B = f(H,0,0)
    C = f(-H,0,0)
    D = f(0,H,0)                     f = length s corresponding to h₁ = -H , 0 , +H ;  h₂ = -H , 0 , +H ;  h₃ = -H , 0 , +H
    E = f(0,-H,0)
    F = f(0,0,H)
    G = f(0,0,-H)
 

-Then, the extremum is estimated by the formula:

   s ~ A - (B-C)² / [ 8.(B+C-2.A) ] - (D-E)² / [ 8.(D+E-2.A) ] - (F-G)² / [ 8.(F+G-2.A) ]
 
 

Data Registers:           •  R00 = L₀               ( If SF 00 , Registers R00 thru R03 are to be initialized before executing "DGT3" )

                                      •  R01 = a                  R04 to R34: temp
                                      •  R02 = b
                                      •  R03 = c

Flags:  F00-F01-F02-F04

   CF 00  6378.171274  6378.101946  6356.751868 & -14.92911 are stored  in registers R01-R02-R03-R00
   SF 00:  you have to initialize R00-R01-R02-R03

   CF 01  CF 02  3 terms in the Fourier series:
   CF 01  SF 02   2 terms in the Fourier series:
   SF 01               1 term in the Fourier series:

   SF 04   geocentric longitudes and latitudes
   CF 04  geodetic longitudes and latitudes

Subroutines:  GX7 & GW7

-Instead if the M-Code routines GX7 & GW7, you can store the coefficients in registers R40 thru R46 ( for example )

   R40 = 0.4179591837
   R41 = 0.3818300505
   R42 = 0.4058451514
   R43 = 0.2797053915
   R44 = 0.7415311856
   R45 = 0.1294849662
   R46 = 0.9491079123

-Then, replace

   lines 158-159 by  RCL IND 24
   lines 148-149 by  RCL IND 24  DSE 24
   lines 141-142 by  46.04
   lines 136 to 138 by  RCL 40
 
 

 
     ( 575 bytes / SIZE 035 )
 
 

 
Example1:      Geodesic distance between the U.S. Naval Observatory at Washington (D.C.)

   l₁ = 77°03'56"0 W ,  b₁ = 38°55'17"2 N   and Cape Town Observatory ( South Africa )
   l₂ = 18°28'35"7 E  ,  b₂ = 33°56'02"5 S

    CF 01  CF 02  CF 04

  -77.0356    ENTER^
   38.55172   ENTER^
   18.28357   ENTER^
  -33.56025  XEQ "DGT3"  >>>>   D₃ =  12709.47906 km = R33                        ---Execution time = 5m46s---
                                            RDN    D₂ = 12709.47906 km
                                            RDN    D₁ = 12709.48047 km
                                            RDN    D₀ = 12709.48074 km = R21

Example2:    With ( -40° , -40° ) ( +120° , +50° ) it yields:

  •  If  CF 01  CF 02  CF 04   we get:

                     D₃ =  18095.49323 km
         RDN    D₂ = 18095.49323 km
         RDN    D₁ = 18095.49858 km
         RDN    D₀ = 18095.55363 km

  •  If  CF 01  CF 02  SF 04   we get:

                     D₃ =  18099.69052 km
         RDN    D₂ = 18099.69052 km
         RDN    D₁ = 18099.69587 km
         RDN    D₀ = 18099.75037 km

Example3:     On the ellipsoid    x² / 41 + y² / 37 + z² / 35 = 1              the geocentric coordinates of the 2 points are:

    l₁ =  9°17'35"55660               l₂ = 63°20'07"3683
    b₁ = -8°11'55"79344              b₂ = 57°46'42"9935

    0   STO 00

   41  SQRT  STO 01
   37  SQRT  STO 02
   35  SQRT  STO 03

    SF 00   CF 01  CF 02  SF 04

    9.173555660  ENTER^
  -8.115579344  ENTER^
   63.20073683  ENTER^
   57.46429935  XEQ "DGT3"  >>>>   D₃ = 8.594823580
                                                 RDN    D₂ = 8.594824326
                                                 RDN    D₁ = 8.594849588
                                                 RDN    D₀ = 8.594880192

-The exact value: about  8.594822580 whence an error of  E-6 with D₃

Notes:

-On the Earth, the maximum error given by D₀ is of the same order as Andoyer's first order formula for an ellipsoid of revolution i-e about 60 meters.

-For the Earth, the difference between D₂ & D₃ is often negligible
-So, you can set F02 without a great loss of precision.

    SF 02  ->  Exectution time = 4m05s
    SF 01  ->  Execution time =  2m29s
 

2°)  Loxodromes
 

-Let  r = semi-major axis of a parallel and  r_m the radius of curvature of the meridian at the point P(L,B)  L = geodetic longitude B = geodetis latitude.

     r² = b² + ( a² - b² ) / [ 1 + (b²/a²) Tan² L ]
   r_m² = c² / { r [ 1 - ( 1 - b²/a² ) Sin² B ]^3/2 }

-Let  a' & b'  the semi-axes of this meridian

    a' = ( a Cos B ) / [ 1 - ( 1 - c²/a² ) Sin² B ]^1/2
    b' = ( b Cos B ) / [ 1 - ( 1 - c²/b² ) Sin² B ]^1/2

-Let r'_m = the radius of curvature of the parallel

    r'_m² = b'² / { a' [ 1 - ( 1 - b'²/a'² ) Sin² L ]^3/2 }

-Then  the azimuth µ of the loxodrome verifies

        r_m Tan µ dB =  r'_m dL
        ( Cos µ ) ds =  r_m dB

-Since r & r'_m depend on L and L is an unknown function of B, we could solve a differential equation to find Tan Az
-But the eccentricity of the equator is very small, and L only appears in a corrective term,
  so we can use the expression L(B) valid for an ellipsoid of revolution

     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₁ ) ] }
 

Data Registers:           •  R00 = L₀               ( If SF 00 , Registers R00 thru R03 are to be initialized before executing "DGT3" )

                                      •  R01 = a                  R04 to R24: temp
                                      •  R02 = b
                                      •  R03 = c

Flags:  F00-F07-F08-F09-F10

   CF 00  6378.171274  6378.101946  6356.751868 & -14.92911 are stored  in registers R01-R02-R03-R00
   SF 00:  you have to initialize R00-R01-R02-R03

Subroutines:  GX7 & GW7

-Instead if the M-Code routines GX7 & GW7, you can store the coefficients in registers R40 thru R46 ( for example )

   R40 = 0.4179591837
   R41 = 0.3818300505
   R42 = 0.4058451514
   R43 = 0.2797053915
   R44 = 0.7415311856
   R45 = 0.1294849662
   R46 = 0.9491079123

-Then, replace

   lines 420-421 by  RCL IND 15
   lines 410-411 by  RCL IND 15  DSE 15
   lines 401 to 403 by  RCL 40
   line  397 by  46.04

-Lines 42-159-216 are three-byte GTOs
 
 

 
     ( 640 bytes / SIZE 025 )
 
 

 
Example1:

   -77.03560  ENTER^
    38.55172  ENTER^
      2.20138  ENTER^
    48.50112  XEQ "LOX3"  >>>>   D = 6453.471591 km                          ---Execution time = 2m40s---
                                            X<>Y   µ  = 80°10'15"8102

Example2:

   0   ENTER^
   0   ENTER^
  70  ENTER^
   0      R/S       >>>>  D = 7792.380608 km                          ---Execution time = 14s---
                       X<>Y   µ = 90°

Example2:

   0  ENTER^
   0  ENTER^
   0  ENTER^
  80    R/S       >>>>  D = 8885.152534 km                          ---Execution time = 15s---
                       X<>Y  µ = 0°

Notes:

-Due to roundoff errors, there is a loss of significant digits if the 2 latitudes are very close but not equal.
-If the latitudes are equal , µ = +/- 90° and if the longitudes are equal , µ = 0 or 180° so the program is much faster
-The azimuths are measured clockwise from North.