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 Andoyer's formula of order 3

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

 4°)  With Vincenty's Formula
 5°)  Minimizing the Sum of 2 Distances
 6°)  A Slightly Different Method

   a)  With Andoyer's formula of order 2
   b)  With Andoyer's formula of order 3
   c)  Almost With Andoyer's formula of order 3
   d)  Geodetic or Geocentric Longitudes & Latitudes ?
 

-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    µ = 2 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 Andoyer's Formula of order 3
 

     a) Andoyer's Formula ( Ellipsoid of Revolution )
 

-Sodano's formula may be re-written as an Andoyer's formula.
 ( In reference [1] - page 87 - Paul D. Thomas gives the key to perform this transformation )

  With      L , L' = longitudes  &  B , B' = parametric latitudes     a = 6378.137 km  ,  b = 6356.752314 km  ( WGS84 )  ,  f = (a-b) / a

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

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

              S = sin² G  cos² H  +  cos² F  sin² H ,  µ = 2 Arc Sin 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 )

-After "some" calculations, it yields:

   Dist = a { µ - (f/2) ( A.µ + B Sin µ )

                + (f²/4) [ A µ² Cot µ + B µ² / Sin µ + A² ( µ + Sin µ Cos µ - 4 µ² Cot µ ) / 4 - A.B µ² / Sin µ - ( B²/2 ) Sin µ Cos µ ]

                + (f³/8) [ A ( 2 µ² Cot µ - 2 µ³ Cot² µ - 8 µ³ / 6 ) + A² ( -5 µ² Cot µ + 5 µ³ Cot² µ + 8 µ³ / 3 + (1/2) Sin µ Cos µ + µ /2 )
                          + A³ ( 3 µ² Cot µ - 3 µ³ Cot² µ - 8 µ³ / 6 + 2 Sin µ Cos³ µ - (9/4) Sin µ Cos µ - µ / 4 + 2 Sin³ µ Cos³ µ + 2 Sin⁵ µ Cos µ )
                          + B ( 2 µ² / Sin µ - 2 µ³ Cos µ / Sin² µ ) + A.B ( -5 µ² / Sin µ + 6 µ³ Cos µ / Sin² µ + µ Cos µ )
                        + A² B ( 3 µ² / Sin µ - 4 µ³ Cos µ / Sin² µ  - 2 Sin³ µ Cos² µ - µ Cos µ - 2 Sin⁵ µ  + 2 Sin µ - (3/2) Sin µ Cos² µ )
                        + B² ( - Sin µ Cos µ + µ³ / Sin² µ  + µ ) + A B² ( (1/2) Sin µ Cos µ - µ - µ³ / Sin² µ )
                        + B³ ( (2/3) Sin³ µ - (1/2) Sin µ ) ] }
 

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

 
     ( 346 bytes / SIZE 015 )
 
 

 
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 "AND"   >>>>   D = 12138.68432 km                                                      ---Execution time = 16s---
 
 

Notes:

-"AND" is slightly shorter & faster than "SOD"
-The formulas are quite equivalent, so the differences in the results are only due to roundoff-errors.

-It is certainly possible to express the formula without calculating the parametric latitudes, but I've not found how.

-In order to calculate the parametric latitudes, we can also use TAN & ATAN instead of P-R & R-P
-But the HP41 gives TAN (-90°) = +9.999999999 E99 instead of -9.999999999 E99
-So, we must use a trick to take this case into account:

-We can replace lines 01 to 32 by
 
 

 
-It costs 2 extra bytes but it saves almost 1 second with an HP41 !
-However, this method cannot be used with free42 because 90 TAN gives OUT OF RANGE
 

     b) Triaxial Ellipsoid
 

-The following "TGD" employs Andoyer's formula of order 3 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
 
 

 
      ( 662 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 = 56s---
                                             X<>Y  d  = 12138.68432 km

Note:

-If you want to use  TAN & ATAN  instead of  P-R & R-P  to calculate the parametric latitudes, replace lines 268 to 290 by
 
 

 
-Just 3 extra bytes, but almost 1 second saved.
 

4°)  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 )
 

5°)  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:  "TGD"  ( cf paragraph 2-b or 3-b above, using 3rd order formulas )
 
 
 

 
      ( 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.89955  km
                                        X<>Y  b" ~ -75°38'

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

   40  CHS  ENTER^  ENTER^
      120       ENTER^
       50        XEQ "D+D"  >>>>   D = 18095.49447  km
                                        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 !
-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 Andoyer's or Sodano's formulas give an error of a few meters.
-Since formulas of order 3 is very accurate when the distance is about 10000 km, the sum will remain very accurate too !

Variant:

-The following routine may be used to find the minimum distance from several values of b"
 
 
 

 
      ( 116 bytes )
 
 

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

   40  CHS  ENTER^  ENTER^
      120       ENTER^
       50        XEQ "D+D"  >>>>   B0=24.1727     i-e  24°17'27"

                               R/S   >>>>    B=?                   if you choose for instance:  b" = 12°  17°  22°  27°  32°

                        12   R/S   >>>>   D=18114.88342      R/S     >>>>    B=?
                        17   R/S   >>>>   D=18102.74082      R/S     >>>>    B=?
                        22   R/S   >>>>   D=18096.49637      R/S     >>>>    B=?
                        27   R/S   >>>>   D=18095.96331      R/S     >>>>    B=?
                        32   R/S   >>>>   D=18101.15912

-Store these 5 results in R01 thru R05 and if you use the programs listed in "Extrema for the HP41" paragraph 1-d)

    5   ENTER^
   22  XEQ "DEX5"  >>>>   b" =  24.96999 = 24°58'12"
                                X<>Y  D = 18095.49452 km

-With a true HP41, it will be less slow than the first "D+D"
-However, it will not be always easy to estimate the proper - and small - interval for b" for nearly antipodal points...
 

6°)  A Slightly Different Method
 

     a) With Andoyer's Formula of Order 2
 

-In the 3 following programs, intead of calculating   D ~ (A/A') d , we compute D ~  A + ( d - A' )
-In other words, the part of the formula giving the great elliptic arc distance on the ellipsoid of revolution in the WGS84 geodesic distance
  is replaced by the great elliptic arc distance on the triaxial ellipsoid.

-The precision is almost the same, but the routines are faster !
 

Formula:

   d - A' = (f²/4) µ [ (A²-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]      where  µ , A , B are computed with parametric latitudes
 
 

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

 
      ( 437 bytes / SIZE 016 )
 
 

   Where  D is the geodesic distance on the triaxial ellipsoid

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.65756 km                                                      ---Execution time = 28s---
 

Notes:

-The root-mean-square error is about 3 cm if  80°  < | L' - L | < 120°  and the maximum error is about 16 cm
-The root-mean-square error is about 6 cm if 110° < | L' - L | < 150°  and the maximum error is about 78 cm
-When the difference in the longitudes is 165°, the error can reach about 4 meters.

-Registers L & R06 contain the great elliptic arc distance.
 

     b) With Andoyer's Formula of Order 3
 

-In this variant, the WGS84 geodesic distance is computed with Andoyer's formula of order 3 ( the other distances with formulas of order 2 )
 

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

 
      ( 628 bytes / SIZE 018 )
 
 

   Where  D is the geodesic distance on the triaxial ellipsoid

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.65757 km                                                      ---Execution time = 35s---
 

Notes:

-The root-mean-square error is about 2 cm if  80°  < | L' - L | < 120°  and the maximum error is about 13 cm
-The root-mean-square error is about 4 cm if 110° < | L' - L | < 150°  and the maximum error is about 34 cm
-When the difference in the longitudes is 165°, the error can reach about 75 cm.

-Registers L & R17 contain the great elliptic arc distance.
 

     c) Almost With Andoyer's Formula of Order 3
 

-We can save about 100 bytes if we only use the larger terms in the 3rd order formula:
-The precision for very long geodesics are of the same order as the program above.
 
 

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

 
      ( 530 bytes / SIZE 016 )
 
 

   Where  D is the geodesic distance on the triaxial ellipsoid

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.65756 km                                                      ---Execution time = 31s---
 

Notes:

-The root-mean-square error is about 3 cm if  80°  < | L' - L | < 120°  and the maximum error is about 15 cm
-The root-mean-square error is about 6 cm if 110° < | L' - L | < 150°  and the maximum error is about 34 cm
-When the difference in the longitudes is 165°, the error can reach about 78 cm.
-I've found these results with free42binary and with latitudes between -87° & +87°

-Registers L & R06 contain the great elliptic arc distance.
 

           d)  Geodetic or Geocentric Longitudes & Latitudes ?
 

-This program employs the same formulas as above, but you can choose the semi-axes of the Earth, and geodetic or geocentric longitudes and latitudes.
 
 

Data Registers:   (•) R00 = L₀                                        ( Registers R00 thru R04 are to be initialized IF Flag F00 is set )

                               (•) R01 = a                    R04 to R16: temp
                               (•) R02 = b
                               (•) R03 = c

Flags:    F00-F03-F04

                               If  F00 is clear, "TGD" uses a = 6378.172 km  ,  b = 6378.102 km  ,  c = 6356.752314 km ,  L₀ = -14°92911
                               If  F00 is set , initialize  R00 to R03

                               If  F04 is clear & F03 is clear , longitudes and latitudes are geodetic
                               If  F04 is clear & F03 is set , longitudes are geocentric but latitudes are geocentric

                               If  F04 is set, longitudes and latitudes are geocentric.

Subroutines: /
 

-Line 190 is a three-byte GTO 14
 
 

 
         ( 592 bytes / SIZE 017 )
 
 

   Where  D is the geodesic distance on the triaxial ellipsoid

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 )

     • If CF 00  CF 03  CF 04  geodetic longitudes & latitudes

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65756 km
                                           LASTX  d = 12138.65876 km = R06 = Great elliptic arc distance

     • If CF 00  SF 03  CF 04   geocentric longitudes & geodetic latitudes

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.70239 km
                                           LASTX  d = 12138.70359 km = R06 = Great elliptic arc distance

     • If CF 00  SF 04   geocentric longitudes & latitudes

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12157.22639 km
                                           LASTX  d = 12157.22759 km = R06 = Great elliptic arc distance
 

Notes:

-The latitudes given on the web are always geodetic, but I don't really know if the longitudes are geocentric or geodetic !
-It's the same for an ellipsoid of revolution, but not for a triaxial ellipsoid...
 
 

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 !