Terrestrial Geodesic Distance(2) for the HP-41

Overview
 

 0°)  Method
 1°)  With Andoyer's formula of order 2
 2°)  With Sodano's Formula - order 3

   a)  Sodano's Formula ( Ellipsoid of Revolution )
   b)  Triaxial Ellipsoid

 3°)  With Vincenty's Formula
 4°)  Minimizing the Sum of 2 Distances
 

-These programs evaluate the geodesic distance between 2 points M(L,B) & N(L',B') on a triaxial ellipsoid model of the Earth.
-The latitudes L and the longitudes B are supposed geodetic.

                                  a = 6378.172 km
-The semi-axes are:    b = 6378.102 km                 and the longitude of the equatorial major-axis = 14°92911 West
                                  c = 6356.752314 km

-So, the mean radius of the equator and the flattening correspond to WGS84.
 

0°)  Method
 

-First, the distances of the 2 points M & N are computed on the great ellipses passing by M & N:

   -For the triaxial model, it gives, say A
   -For the WGS84 ellipsoid, it gives  A'

-Then, the geodesic distance d is calculated on the WGS84 ellipsoid of revolution.

-Since the flattening of the Earth equator is very small, the geodesic distance D on the triaxial ellipsoid is very well approximated by the simple relation:

    D = (A/A') d
 

1°)  With Andoyer's Formula of order 2
 

"TGD" determines the elements of the great ellipse (E) passing by M & N as follows:

   With  OM / OM = ( x , y , z )   & ON / ON = ( x' , y' ,z' )      ( the unit vectors defined by OM & ON )     we have

   OP = ( Sin W ) / ( A² + B² + C² )^1/2    for a point  P  on the ellipsoid  and   with  µ = ( OM , OP )

      A = [ x Sin ( W - µ ) + x' Sin µ ] / a
      B = [ y Sin ( W - µ ) + y' Sin µ ] / b           W = the angle ( OM , ON )
      C = [ z Sin ( W - µ ) + z' Sin µ ] / c
 

                         N  P
                        /   /
                      /  /
                    / / µ
               O//-----------------------------M
 

>>> The value of µ for the axes of the great ellipse (E) is found by equating to zero the derivative of OP with respect to µ

-After some calculations, it yields:    ( Sin 2µ ) / ( Cos 2µ ) = P / Q  with

    P  =  -2 [ ( x Sin W ) ( x Cos W - x' ) / a² + ( y Sin W ) ( y Cos W - y' ) / b² + ( z Sin W ) ( z Cos W - z' ) / c² ]
    Q =  [ x² Sin² W - ( x Cos W - x' )² ] / a² + [ y² Sin² W - ( y Cos W - y' )² ] / b² + [ z² Sin² W - ( z Cos W - z' )² ] / c²

>>> However, it's simpler & faster to compute

    A² + B² + C² = T Sin² ( W - µ ) + S Sin² µ + 2 U Sin µ Sin ( W - µ )

   where  S = x'² / a² + y'² / b² + z'² / c²  ,   T =  x² / a² + y² / b² + z² / c²  ,   U =  xx' / a² + yy' / b² + zz' / c²

    and   P = 2 U - 2 T Cos W  ,  Q = T Sin W + ( 2 U Cos W - S - T Cos² W ) / Sin W
 

-Then, the geodesic distances - on the great ellipses corresponding to WGS84 & the triaxial ellipsoid - between M & N ( say A' & A )
   are computed by Andoyer's formula on order 2.

-Finally, the geodesic distance d on the WGS84 ellipsoid is calculated by 2nd order Andoyer's formula
 and the geodesic distance D on the triaxial ellipsoid is approximated by   D = (A/A') d
 

Andoyer's Formula of order 2:
 

-With   H = ( L' - L )/2  ;   F = ( B + B' )/2  ;  G = ( B' - B )/2

              S = sin² G  cos² H  +  cos² F  sin² H  ;   µ = Arcsin S^1/2  ;  R = ( S.(1-S) )^1/2

 we define  A = ( sin² G  cos² F )/S + ( sin² F  cos² G )/( 1-S )
      and      B = ( sin² G  cos² F )/S - ( sin² F  cos² G )/( 1-S )

-The geodesic distance on the biaxial ellipsoid is then:   d = 2a.µ ( 1 + eps )

  where   eps =    (f/2) ( -A - 3B R/µ )
                       + (f²/4) { { - ( µ/R + 6R/µ ).B + [ -15/4 + ( 2S - 1 ) ( µ/R + (15/4) R/µ ) ] A + 4 - ( 2S - 1 ) µ/R } A
                                      - [ (15/2) ( 2S - 1 ).B  R/µ - ( µ/R + 6R/µ ) ] B }

   and   f = flattening = 1/298.25722
 

Data Registers:   R00  to  R15:  temp
Flags: /
Subroutines: /
 

-Line 38 is a three-byte GTO 14
 
 

 
      ( 471 bytes / SIZE 016 )
 
 

   Where  D is the geodesic distance on the triaxial ellipsoid
     and     d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )

Example:    Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
                    and the Palomar Observatory ( Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N )

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65748 km                                                      ---Execution time = 47s---
                                             X<>Y  d  = 12138.68425 km

-The exact results are D = 12138.65757  &  d = 12138.68433 km  ( rounded to 5 decimals )
 

Notes:

-Longitudes are geodetic and positive East.
-If the points are not nearly antipodal, the precision is of the order of one meter.

-With (L,b,L',b') = (70°,0°,270°,10°) "TGD" gives a distance of  17554.6159  instead of  17554.6140 km

-Here are a few results about the error of the results:

  If | L' - L | < 120°  rms error ~  6 cm   and  max error ~  26 cm
  If | L' - L | < 140°  rms error ~  7 cm   and  max error ~  38 cm
  If | L' - L | < 150°  rms error ~  8 cm   and  max error ~  70 cm
  If | L' - L | < 160°  rms error ~ 11 cm  and  max error ~ 168 cm

-A better accuracy may be obtained by the programs listed in "Geodesics on a triaxial ellipsoid for the HP-41" but the program above is faster !
 

2°)  With Sodano's Formula of order 3
 

     a) Sodano's Formula ( Ellipsoid of Revolution )
 

-"SOD" uses a 3rd order formula in the flattening to calculate the geodesic distance on an ellipsoid of revolution ( WGS84 ) - cf references [2] & [3]

-I've rewritten the formula as follows, after multiplying both sides by (1-f)

  Dist = a { µ + f  [ p Sin µ - (m/2) ( µ + Sin µ Cos µ )

          +  f² [ p ((m-1)/2) µ² Csc µ + m ((1-m)/2) µ² Cot µ + ((m²-8.p²)/16) Sin µ Cos µ + (m²/16) µ - (m²/8) Sin µ Cos³ µ + (mp/2) Sin µ Cos² µ ]

          +  f³ [ µ² Csc µ (p/4) (-2+5m-3m²) + µ³ Csc µ Cot µ (p/2) (1-4m+3m²) - µ² Cot µ (m/4) (-2+5m-3m²) + µ³ (m/6) (-1+2m-m²)
                + µ³ Cot² µ (m/2) (-1+3m-2m²) + Sinµ Cos µ (-8 p² + m² - 8 p²m - m³/2 )/16 + µ  (8 p² + m² - 8 p²m - m³/2 )/16
                + µ³ Csc² µ (p²/2) (1-m) + Sin µ Cos³ µ (m²/16) (-2+m) + µ Cos² µ (m²/4) (1-m) + Sin µ Cos² µ (pm/2) (1+m) + p²m Sin³ µ Cos µ
                + µ Cos µ (3pm/4) (-1+m) + (pm²/2) Sin⁵ µ + Sin µ (p/2) (p²-m²) + (m³/12) Sin³ µ Cos³ µ - (2p³/3) Sin³ µ ]  }

  where    µ = Arc Sin [ sin² ((B'-B)/2)  cos² ((L'-L)/2)  +  cos² ((B+B')/2)  sin² ((L'-L)/2) ]^1/2     L , L' = longitudes  &  B , B' = parametric latitudes

               Tan B = (1-f) Tan l  where  l = geodetic latitude

               p = Sin B  Sin B'       m = 1 - [ Cos B Cos B' Sin (L'-L) ]² / Sin² µ   ,   a = 6378.137 km ,  b = 6356.752314 km  ( WGS84 )
 
 

Data Registers:   R00 thru R13: temp
Flags: /
Subroutines: /
 
 
 

 
     ( 374 bytes / SIZE 014 )
 
 

 
Example:    Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
                    and the Palomar Observatory ( Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N )

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "SOD"   >>>>   D = 12138.68432 km                                                      ---Execution time = 19s---
 
 

Notes:

-If  | L' - L | < 120° , the RMS error is about 1.4 mm and the maximum error about 9 mm
-If  | L' - L | < 170° , -----------------------  12 mm  ----------------------------- 22 cm
-If  | L' - L | < 175° , -----------------------   5  cm   -----------------------------  2 m

-When Vincenty's iterative method requires more than 3 iterations, even a 3rd order formula may produce an error of 2 meters:
-For example, with  Long1 = 0° , Lat1 = 2° &  Long2 = 175° , Lat2 = -4°

   "SOD"  gives  19434.163468 km  instead of  19434.161314 km  returned by Vincenty's formula !

-These estimations were obtained with free42.
-With an HP41 - which works with 10 digits - the errors may be slightly larger...
 

     b) Triaxial Ellipsoid
 

-This variant employs Sodano's formula for the last calculation and Andoyer's formula of order 2 to estimate  (A/A')
 

Data Registers:   R00 thru R17: temp
Flags: /
Subroutines: /

-Line 38 is a three-byte GTO 14
 
 

 
      ( 699 bytes / SIZE 018 )
 
 

   Where  D is the geodesic distance on the triaxial ellipsoid
     and     d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )

Example:    Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
                    and the Palomar Observatory ( Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N )
 

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                      ---Execution time = 59s---
                                             X<>Y  d  = 12138.68432 km
 

3°)  With Vincenty's Formula
 

-Vincenty's formula is probably the best method for an ellipsoid of revolution.
-The following version of "TGD" - renamed "TGDV"-  employs this method 3 times.
 

Data Registers:   R00 thru R21: temp
Flags: /
Subroutines: /

-Line 47 is a three-byte GTO 14
 
 

 
      ( 564 bytes / SIZE 022 )
 
 

   Where  D is the geodesic distance on the triaxial ellipsoid
     and     d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )

Example:    Calculate the geodesic distance between the Sydney Observatory ( Long = 151°12'17"8 E , Lat = 33°51'41"1 S )
                    and the Palomar Observatory ( Long = 116°51'50"4 W ,  Lat = 33°21'22"4 N )
 

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGDV"   >>>>   D = 12138.65757 km                                                      ---Execution time = 119s---
                                                X<>Y  d  = 12138.68433 km
 

Notes:

-This version gives the best accuarcy, provided the points are not ( almost ) antipodal.
-For nearly antipodal points, the routine stops and the HP41 displays DATA ERROR ( line 235 )
 

4°)  Minimizing the Sum of 2 Distances
 

-No one of the programs above works well for nearly antipodal points.
-The following one uses the golden search to minimize the sum of 2 geodesic distances:

  Let M(L,b) & N(L',b')  ,  "D+D"  minimizes the sum  D = Dist(M,P) + Dist(P,N)

  where  P(L",b")  with  L" = the longitude of the point on a sphere in the middle of the geodesic MN on the sphere.

-The program starts with  b" = -90° and +90°  and employs the golden search to find the minimum D-value.
 
 

Data Registers:   R00 thru R35: temp
Flags:  /
Subroutine:  "DGT"  ( cf paragraph 2-b above, using Sodano's 3rd order formula )
 
 
 

 
      ( 181 bytes / SIZE 036 )
 
 

   Where  D is the geodesic distance on the triaxial ellipsoid
     and    b" is the latitude of the "midpoint"

Example1:   With (L,b) = (0°,0°)  &  (L',b') = (179.51°,0°)

        0         ENTER^  ENTER^
    179.51    ENTER^
        0         XEQ "D+D"  >>>>   D = 20001.89956  km                    ---Execution time = 43m02s---
                                        X<>Y  b" ~ 75°34'

Example2:   With (L,b) = (-40°,-40°)  &  (L',b') = (120°,50°)

   40  CHS  ENTER^  ENTER^
      120       ENTER^
       50        XEQ "D+D"  >>>>   D = 18095.49438  km                    ---Execution time = 40m56s---
                                        X<>Y  b" ~ 24°58'
 

Notes:

-The HP41 displays the successive approximations.

-The angle b" is not computed accurately: it's not necessary to get a good precision for the distance.
-You can press  XEQ 10  to make another loop, or - for instance - replace lines 98-99 by a larger number.

-This program is very slow with a true HP41, so a good emulator is not superfluous !
-"DGT" is called 44 ( respectively 42 ) times in the examples above.
-With free42, you get the results 20001.89954 & 18095.49449 almost instantaneously...

-This program may also be used for very long geodesics when the points are not almost antipodal, but when Sodano's formula gives an error of a few meters.
-Since Sodano's formula is very accurate when the distance is about 10000 km, the sum will remain very accurate too !
 
 
 
 

References:

[1]  Paul D. Thomas - "Spheroidal Geodesics, Reference Systems & Local Geometry" - US Naval Oceanographic Office.
[2]  Richard H. Rapp - "Geometric Geodesy, Part II - The Ohio State University
[3]  Emanuel M. Sodano & Thelma A. Robinson - "Direct and Inverse Solutions of Geodesics" - Army Map Service, Technical Report n°7

-Many thanks to Professor Richard H. Rapp who sent me reference [3] to find the typo in Sodano's formula !