Geodesic

# Terrestrial Geodesic Distance for the HP-41

Overview

0°) Sphere
1°) Andoyer's Formulae

a) Geodesic Distance only

a-1) Program#1 ( First order accuracy )
a-2) Program#2 ( Second order accuracy )
a-3) Nearly Antipodal Points

b) Geodesic Distance and Azimuths

b-1) Program#1
b-2) Program#2
b-3) Program#3
b-4) Program#4
b-5) Program#5

c) Forward Geodesic Problem

2°) Vincenty's Formulae

a) Geodesic Distance and Azimuths
b) Forward Geodesic Problem

3°) Numerical Integration
4°) Series Expansions

a) Program#1
b) Program#2 ( better precision )

5°) Euclidean Distance

a) Geocentric Coordinates
b) 3D-Distance

------------------------------------------------------------------------------------------------------------

Latest Update:

-Paragraph 0°)

------------------------------------------------------------------------------------------------------------

-The programs listed in paragraphs 1 to 4 solve geodesic problems  between 2 points.
-These points are defined by their geographic coordinates:  ( L1 , b1 )  ( L2 , b2 )   where   Li = Longitudes ,  bi = latitudes.
-The geodesic distance D between these 2 points P1 & P2 ( on the Earth's surface , at mean sea level ) is expressed in kilometers.
-The azimuths  are expressed in  ° ' "
-An ellipsoidal model of the Earth is used with

a = equatorial radius = 6378.137 km     ( DE403 )
f  = Earth's flattening = 1/298.257       ( IERS 1992 )

-WGS84 uses   a = 6378.137 km  and  f = 1/298.257223563
-More recent determinations suggest:

a = 6378.13661 km & f = 1/298.256421   Change the corresponding lines if you prefer these values.

-If you only want to compute the geodesic distance, the longitudes can be measured positively westward or positively eastward from the meridian of Greenwich.

-Otherwise, the longitudes are measured positively Eastwards and the azimuths are reckoned clockwise from North.

-In paragraph 5, the observers' heights are taken into account to compute the Euclidean 3D-distance.

0°)  Sphere

-With  a = 6371 km

Formula:

Sin2 µ/2  = sin2 G  cos2 H  +  cos2 F  sin2 H      F = ( b + b' ) / 2
Cos2 µ/2 = cos2 G  cos2 H  +  sin2 F  sin2 H      G = ( b - b' ) / 2             H = ( L - L' ) / 2

D = a.µ

Data Registers: /
Flags: /
Subroutines: /

 01 LBL "DGT"  02 DEG  03 X<> T            04 HMS-  05 HR  06 2  07 /  08 1  09 P-R 10 R^  11 HR  12 STO 00            13 R^  14 HR  15 ST+ 00  16 -  17 2  18 ST/ 00 19 /  20 X<>Y  21 P-R  22 X^2  23 X<>Y  24 X^2  25 RCL 00            26 R^  27 P-R 28 X^2  29 ST+ Z  30 X<> Z            31 SQRT  32 X<>Y  33 X^2  34 R^  35 + 36 SQRT  37 R-P  38 X<>Y  39 D-R  40 ST+ X            41 6371  42 *  43 END

( 60 bytes / SIZE 001 )

 STACK INPUTS OUTPUTS T L ( ° ' " ) / Z b ( ° ' " ) / Y L' ( ° ' " ) / X b' ( ° ' " ) D ( km )

Example:   ( L , b ) = ( -77°03'56" , 38°55'17 )    ( L' , b' ) = ( 2°20'14" , 48°50'11" )

-77.0356    ENTER^
38.5517   ENTER^
2.2014    ENTER^
48.5011   XEQ "DGT"  >>>>   D = 6165.597256 km

1°) Andoyer's Formulae

a) Geodesic Distance only

a-1) Program#1 ( First order accuracy )

-Let                    L = ( L2 - L1 )/2  ;   F = ( b1 + b2 )/2  ;  G = ( b2 - b1 )/2
-We calculate     S = sin2 G  cos2 L  +  cos2 F  sin2 L  ;   µ = Arcsin S1/2  ;  R = ( S.(1-S) )1/2

-Then,                D = a.{ 2µ + f.[ (3R-µ)/(1-S)  sin2 F  cos2 G  - (3R+µ)/S  sin2 G  cos2 F ] }

Data Registers: /
Flags: /
Subroutines: /

-Synthetic register M may be replaced by any data register like R00 ...

 01  LBL "TGD" 02  DEG 03  HR 04  X<> T  05  HMS- 06  HR 07  2 08  / 09  SIN 10  X^2 11  RDN 12  HR 13  + 14  2 15  / 16  ST- Y 17  COS 18  X^2 19  ST* T  20  RDN 21  SIN 22  X^2 23  STO M       24  ST* Y 25  - 26  - 27  ENTER^ 28  SQRT 29  ASIN 30  D-R 31  RCL M  32  R^ 33  ST* M  34  + 35  RCL M       36  - 37  1 38  ST- Y 39  R^ 40  ST/ M 41  - 42  ST/ Y 43  LASTX        44  * 45  SQRT 46  3 47  * 48  R^ 49  ST+ T  50  - 51  ST* Y 52  R^ 53  + 54  0 55  X<> M        56  * 57  + 58  298.257 59  / 60  - 61  6378.137 62  * 63  END

( 98 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) / Y L2 ( ° ' " ) / X b2 ( ° ' " ) D ( km )

Example:   Calculate the geodesic distance between the U.S. Naval Observatory at Washington (D.C.)    L1 = 77°03'56"0 W ;  b1 = 38°55'17"2 N
and the Observatoire de Paris  L2 = 2°20'13"8 E ;  b2 = 48°50'11"2 N

77.0356    ENTER^
38.55172  ENTER^
-2.20138   ENTER^
48.50112  XEQ "TGD"  >>>>   D = 6181.6156 km   ( execution time = 5 seconds )

-In this example, the error is about  6 meters.
-If the points are not ( nearly ) antipodal, the relative error is less than  10 -5 and the maximum error is of the order of 60 meters in absolute value:
with ( 0° , -40° ) ( 120° , 40° )   error = 66 meters
-However, errors can reach several kilometers for nearly antipodal points.
-If the 2 points are on the equator:  ( L1 , 0 ) ( L2 , 0 ) , Andoyer's formula gives a perfect accuracy provided  |  L2 - L1 |  <  179°23'47"377

a-2) Program#2 ( Second order accuracy )

-This program uses a formula which is second-order in the flattening f:

-With   L = ( L2 - L1 )/2  ;   F = ( b1 + b2 )/2  ;  G = ( b2 - b1 )/2
S = sin2 G  cos2 L  +  cos2 F  sin2 L  ;   µ = Arcsin S1/2  ;  R = ( S.(1-S) )1/2

we define  A = ( sin2 G  cos2 F )/S + ( sin2 F  cos2 G )/( 1-S )
and      B = ( sin2 G  cos2 F )/S - ( sin2 F  cos2 G )/( 1-S )

-The geodesic distance is then   D = 2a.µ ( 1 + eps )

where   eps =    (f/2) ( -A - 3B R/µ )
+ (f2/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 }

Data Registers:   R00 = µ/R    R01 = A    R02 = B    R03 = 2S-1    R04 = 2/f = 596.514
Flags: /
Subroutines: /

 01 LBL "TGD1"  02 DEG  03 HR  04 X<> T  05 HMS-  06 HR  07 2  08 /  09 COS  10 X^2  11 RDN  12 HR  13 +  14 2  15 /  16 ST- Y  17 1  18 P-R  19 X^2  20 STO 02            21 RDN 22 X^2  23 STO 01  24 SIGN  25 P-R  26 X^2  27 ST* 01  28 RDN  29 X^2  30 ST* 02  31 RCL Z  32 -  33 *  34 +  35 ST/ 02  36 STO 03  37 SQRT  38 ASIN  39 D-R  40 RCL 03            41 1  42 RCL 03 43 -  44 ST/ 01  45 ST- 03  46 *  47 SQRT  48 /  49 STO 00  50 CLX  51 RCL 02  52 RCL 01  53 ST- 02  54 +  55 STO 01  56 RCL 03  57 RCL 00            58 ST* 03  59 /  60 7.5  61 *  62 STO 04  63 LASTX 64 -  65 2  66 /  67 RCL 03  68 +  69 *  70 4  71 +  72 RCL 03  73 -  74 6  75 RCL 00  76 ST/ Y  77 +  78 RCL 02            79 *  80 STO Z  81 -  82 RCL 01  83 *  84 +  85 RCL 02 86 X^2  87 RCL 04  88 *  89 -  90 596.514  91 STO 04  92 /  93 RCL 01  94 -  95 RCL 02  96 3  97 *  98 RCL 00            99 / 100 - 101 RCL 04 102 / 103 * 104 + 105 12756.274 106 * 107 END

( 144 bytes / SIZE 005 )

 STACK INPUTS OUTPUTS T L ( ° ' " ) / Z b ( ° ' " ) / Y L' ( ° ' " ) / X b' ( ° ' " ) D ( km )

Example:

77.0356    ENTER^
38.55172  ENTER^
-2.20138   ENTER^
48.50112  XEQ "TGD1"  >>>>   D = 6181.621809 km                                        ---Execution time = 7s---

Notes:

-The error is only 15 mm!

-The precision is good ( less than 1.50 meter if | L-L' | < 160° ) except for nearly antipodal points.
-Otherwise, the errors may be larger: for example, ( 0° , 0° ) ( 179° , 1° )  give 19860.3438 km whereas the correct result is 19860.5092 km

-Cf "Terrestrial Geodesic Distance(2)" for a 3rd order Andoyer's formula.

-This version doesn't work if  L - L'  is very small:  DATAERROR line 35
-Here is another version that works well in this case ( but 2 extra bytes )

 01 LBL "TGD1"  02 DEG  03 HR  04 X<> T  05 HMS-  06 HR  07 2  08 /  09 SIN  10 X^2  11 RDN  12 HR  13 +  14 2  15 /  16 ST- Y  17 1  18 P-R  19 X^2  20 STO 02            21 X<>Y  22 X^2 23 STO 01  24 RDN  25 X<>Y  26 1  27 P-R  28 X^2  29 ST* 01  30 RDN  31 X^2  32 ST* 02  33 STO T  34 -  35 *  36 +  37 ST/ 02  38 STO 03  39 SQRT  40 ASIN  41 D-R  42 RCL 03            43 1  44 RCL 03 45 -  46 ST/ 01  47 ST- 03  48 *  49 SQRT  50 /  51 STO 00  52 CLX  53 RCL 02  54 RCL 01  55 ST- 02  56 +  57 STO 01  58 RCL 03  59 RCL 00  60 ST* 03  61 /  62 7.5  63 *  64 STO 04            65 LASTX  66 - 67 2  68 /  69 RCL 03  70 +  71 *  72 4  73 +  74 RCL 03  75 -  76 6  77 RCL 00  78 ST/ Y  79 +  80 RCL 02  81 *  82 STO Z  83 -  84 RCL 01            85 *  86 +  87 RCL 02  88 X^2 89 RCL 04  90 *  91 -  92 596.514  93 STO 04  94 /  95 RCL 01  96 -  97 RCL 02  98 3  99 * 100 RCL 00 101 / 102 - 103 RCL 04           104 / 105 * 106 + 107 12756.274 108 * 109 END

( 146 bytes / SIZE 005 )

 STACK INPUTS OUTPUTS T L ( ° ' " ) / Z b ( ° ' " ) / Y L' ( ° ' " ) / X b' ( ° ' " ) D ( km )

Example:

77.0356    ENTER^
38.55172  ENTER^
-2.20138   ENTER^
48.50112  XEQ "TGD1"  >>>>   D = 6181.621809 km                                        ---Execution time = 7s---

a-3) Nearly Antipodal Points

-"NAPP" uses an intermediate point and adjust it so that the sum of the geodesic distances is minimum.
-It calls "TGD1" as a subroutine. "TGD1" may be replaced by "TGD" if a lesser accuracy is acceptable.

Data Registers:   R00 thru R04 are used by "TGD1"  ,   R05 to R13: temp   R10 = sum of the 2 distances.
Flags:   F01 & F02
Subroutine:  "TGD1"

 01  LBL "NAPP"  02  STO 06   03  X<> T  04  HMS-  05  HR  06  ABS  07  PI  08  R-D  09  X>Y?  10  CLX  11  ST+ X  12  -  13  ABS  14  HMS  15  STO 07  16  X<> L  17  2  18  /  19  HMS 20  STO 08  21  X<>Y  22  STO 05          23  20  24  STO 11  25  CLX  26  STO 09  27  XEQ 02  28  STO 10  29  LBL 00  30  CF 01  31  CF 02  32  LBL 01  33  VIEW 10  34  FS? 02  35  GTO 00  36  RCL 09  37  RCL 11  38  + 39  XEQ 02  40  RCL 10  41  X<=Y?  42  GTO 00         43  X<>Y  44  STO 10   45  RCL 13  46  STO 09  47  SF 01  48  GTO 01  49  LBL 00  50  FS? 01  51  GTO 00  52  RCL 09  53  RCL 11  54  -  55  XEQ 02  56  RCL 10  57  X<=Y? 58  GTO 00  59  X<>Y  60  STO 10  61  RCL 13  62  STO 09   63  SF 02  64  GTO 01         65  LBL 00  66  PI  67  ST/ 11  68  .01  69  RCL 11  70  X>Y?  71  GTO 00  72  RCL 10  73  RTN  74  LBL 02  75  STO 13  76  CLST 77  RCL 05  78  RCL 08  79  RCL 13   80  HMS  81  XEQ "TGD1"  82  STO 12  83  RCL 08  84  RCL 13  85  HMS  86  RCL 07  87  RCL 06  88  XEQ "TGD1"  89  RCL 12  90  +  91  RTN  92  END

( 138 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) / Y L2 ( ° ' " ) / X b2 ( ° ' " ) D ( km )

Example:      L1 = b1 = 0  ;   L2 = 179°51'   b2 = 0

0      ENTER^  ENTER^
179.51  ENTER^
0      XEQ "NAPP"  >>>>   D = 20001.85456 km                             ---Execution time = 4mn30s---

Notes:

-It's obviously a slow method ... but it works !
-Line 68 , the value 0.01 seems small enough for nearly antipodal points.
-Otherwise, a smaller constant could be preferable.

b) Geodesic Distance and Azimuths

-With the same notations as above, the forward azimuth Az1 in the direction from P1 towards P2 is obtained by the following formulae:    Az1 = Az1sph + dAz

where  dAz = (f/2).( cos2 F - sin2 G ).[ ( 1+µ/R ).( sin G cos F )/S + ( 1-µ/R ).( sin F cos G )/(1-S) ].( sin 2L )

and    Tan Az1sph = cos b2 sin 2L / ( cos b1 sin b2 - sin b1 cos b2 cos 2L )     >>>  the R-P function returns the angle in the proper quadrant

Az1sph is the forward azimuth in a spherical model of the Earth,
dAz is a correction that takes the Earth's flattening f into account.

-The back azimuth ( from P2 to P1 ) is obtained by replacing L by -L and exchanging b1 & b2 in the formulas above.

b-1) Program#1

-This routine computes D with the 1st-order formula used by "TGD"

Data Registers:   R00 = 2L , µ , µ/R         R03 = ( sin F cos G )/(1-S)               R06 = (f/2).( cos2 F - sin2 G ).( sin 2L )
R01 = Az1sph               R04 =  ( sin G cos F )/S , temp          R07 = 1 - S
R02 = Az2sph               R05 = S
Flags: /
Subroutines: /

 01  LBL "TGD+"   02  DEG   03  HR   04  STO 03    05  X<> T   06  HMS-   07  HR   08  STO 00   09  X<>Y   10  HR   11  STO 01   12  COS   13  STO 04    14  CHS   15  P-R   16  R^   17  SIN   18  ST* 04   19  *   20  RCL 01    21  SIN   22  STO 05   23  R^   24  ST+ 01   25  COS   26  STO 02 27  *   28  +   29  R-P   30  RCL 00    31  SIN   32  STO 06          33  LASTX   34  R^   35  X<> 02   36  P-R   37  RCL 05   38  *   39  RCL 04    40  X<>Y   41  -   42  R-P   43  X<>Y   44  X<> 01   45  2   46  ST/ 00   47  /   48  ST- 03   49  1   50  P-R   51  STO 04    52  X^2 53  X<>Y   54  X<> 03   55  1   56  P-R   57  ST* 03    58  RDN   59  ST* 04   60  X^2   61  STO Z   62  -   63  ST* 06   64  RCL 00          65  SIN   66  X^2   67  *   68  +   69  STO 05   70  SQRT   71  ASIN   72  D-R   73  STO 00   74  RCL 05   75  1   76  RCL 05    77  -   78  STO 07 79  *   80  SQRT   81  ST/ 00   82  3   83  *   84  ST- Z   85  +   86  RCL 04    87  X^2   88  *   89  RCL 05          90  ST/ 04   91  /   92  X<>Y   93  RCL 03   94  X^2   95  *   96  RCL 07   97  ST/ 03   98  /   99  + 100  596.514 101  ST/ 06 102  / 103  - 104  RCL 00 105  RCL 04 106  ST* Y 107  + 108  STO 04 109  RCL 00  110  RCL 03 111  ST* Y 112  - 113  ST- 04 114  + 115  RCL 06        116  R-D 117  ST* 04 118  * 119  RCL 02 120  + 121  HMS 122  RCL 01 123  RCL 04  124  + 125  HMS 126  R^ 127  12756.274 128  * 129  END

( 173 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) Az2 ( ° ' " ) Z b1 ( ° ' " ) Az2 ( ° ' " ) Y L2 ( ° ' " ) Az1 ( ° ' " ) X b2 ( ° ' " ) D ( km )

Example:

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD+"  >>>>    D   =   6181.6156 km                                                                   ---Execution time = 13s---
RDN    Az1 =  51°47'36"76 = forward azimuth      Az1-error = -0"05
RDN    Az2 = -68°09'59"26 = back azimuth          Az2-error = -0"29

Notes:

-Like "TGD" , "TGD+" gives accurate results if the 2 points are not nearly antipodal.
-Az is obtained with an accuracy of the order of a few arcseconds ( the error is smaller than 1" in this example ).
-If the distance doesn't exceed 17000 km, the maximum errors I've remarked are 66 meters for D and 55" for the azimuths.

-Otherwise, the results may still be accurate:       with ( 0° , -60° ) ( 170° , 70° )   D-error = 13m ,    Az1-error = 2" ,       Az2-error = 3"
but in other examples, it may be disappointing:  with ( 0° , 0° )    ( 179° , 1° )     D-error = 3303m , Az1-error = 49'06" , Az2-error = 49'09"

b-2) Program#2

-This program computes D with the 2nd-order formula used by "TGD1"

Data Registers:   R00 = 2L , µ , µ/R         R03 = ( sin F cos G )/(1-S)               R06 = (f/2).( cos2 F - sin2 G ).( sin 2L )          R08 = A
R01 = Az1sph               R04 =  ( sin G cos F )/S , temp          R07 = 1 - S , temp , 2/f = 596.514                 R09 = B
R02 = Az2sph               R05 = S , 2S-1
Flags: /
Subroutines: /

 01  LBL "TGD1+"   02  DEG   03  HR   04  STO 03    05  X<> T   06  HMS-   07  HR   08  STO 00    09  X<>Y   10  HR   11  STO 01   12  COS   13  STO 04   14  CHS   15  P-R   16  R^   17  SIN   18  ST* 04   19  *   20  RCL 01   21  SIN   22  STO 05   23  R^   24  ST+ 01   25  COS   26  STO 02   27  *   28  +   29  R-P   30  RCL 00 31  SIN   32  STO 06    33  LASTX   34  R^   35  X<> 02   36  P-R   37  RCL 05            38  *   39  RCL 04    40  X<>Y   41  -   42  R-P   43  X<>Y   44  X<> 01   45  2   46  ST/ 00   47  /   48  ST- 03   49  1   50  P-R   51  STO 04   52  X^2   53  X<>Y   54  X<> 03   55  1   56  P-R   57  ST* 03   58  RDN   59  ST* 04   60  X^2 61  STO Z   62  -   63  ST* 06   64  RCL 00    65  SIN   66  X^2   67  *   68  +   69  STO 05            70  SQRT   71  ASIN   72  D-R   73  STO 00   74  RCL 05   75  1   76  RCL 05   77  -   78  STO 07   79  *   80  SQRT   81  ST/ 00   82  CLX   83  RCL 04   84  X^2   85  RCL 05   86  ST/ 04   87  /   88  STO 09   89  RCL 03   90  X^2 91  RCL 07    92  ST- 05   93  ST/ 03   94  /   95  ST- 09   96  +   97  STO 08    98  3.75   99  STO 07          100  RCL 00 101  ST/ Y 102  + 103  RCL 05 104  * 105  RCL 07 106  - 107  * 108  6 109  RCL 00 110  ST/ Y 111  + 112  STO 07 113  RCL 09 114  * 115  - 116  4 117  + 118  RCL 00 119  RCL 05 120  * 121  - 122  RCL 08  123  * 124  RCL 05 125  RCL 09  126  * 127  7.5 128  * 129  RCL 00          130  / 131  RCL 07 132  - 133  RCL 09 134  * 135  - 136  596.514 137  STO 07 138  ST/ 06 139  / 140  RCL 08 141  - 142  RCL 09 143  3 144  * 145  RCL 00 146  / 147  - 148  RCL 07 149  / 150  * 151  + 152  RCL 00  153  RCL 04 154  ST* Y 155  + 156  STO 04 157  RCL 00          158  RCL 03 159  ST* Y 160  - 161  ST- 04 162  + 163  RCL 06 164  R-D 165  ST* 04 166  * 167  RCL 02 168  + 169  HMS 170  RCL 01 171  RCL 04 172  + 173  HMS 174  R^ 175  12756.274 176  * 177  END

( 230 bytes / SIZE 010 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) Az2 ( ° ' " ) Z b1 ( ° ' " ) Az2 ( ° ' " ) Y L2 ( ° ' " ) Az1 ( ° ' " ) X b2 ( ° ' " ) D ( km )

Example:

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD1+"  >>>>    D   =   6181.621809 km                              ---Execution time = 15s---
RDN    Az1 =  51°47'36"76 = forward azimuth
RDN    Az2 = -68°09'59"26 = back azimuth

b-3) Program#3

-The formulae that compute the azimuths on a sphere may be re-written to save bytes and execution time, but it uses more data registers:

Tan ( Az1sph  - Az2sph )/2 = [ ( cos L )( cos G ) ] / [ (sin L )( sin F ) ]
Tan ( -Az1sph - Az2sph )/2 = [ ( cos L )( sin G ) ] / [ (sin L )( cos F ) ]

R-P  returns the angles in the proper quadrants.

-The precision is the same as above.

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

 01  LBL "TGD"   02  DEG   03  HR   04  X<> T   05  HMS-   06  HR   07  STO 09    08  2   09  /   10  1   11  P-R   12  STO 10   13  STO 12   14  CLX   15  +   16  STO 11   17  STO 13   18  X^2   19  RDN   20  HR   21  +   22  2   23  /   24  ST- Y   25  1   26  P-R   27  STO 08   28  ST* 13 29  X^2   30  STO 02   31  X<>Y   32  STO 07   33  ST* 11   34  X^2   35  STO 01         36  RDN   37  X<>Y   38  1   39  P-R   40  ST* 07   41  ST* 10   42  X^2   43  ST* 01   44  RDN   45  ST* 08   46  ST* 12   47  X^2   48  ST* 02   49  STO T   50  -   51  STO 06   52  *   53  +   54  ST/ 02   55  STO 03   56  SQRT 57  ASIN   58  D-R   59  RCL 03   60  ST/ 08   61  1   62  RCL 03         63  -   64  ST/ 01   65  ST/ 07   66  ST- 03   67  *   68  SQRT   69  /   70  STO 00   71  CLX   72  RCL 02   73  RCL 01   74  ST- 02   75  +   76  STO 01   77  RCL 03   78  STO 04   79  RCL 00   80  ST* 04   81  /   82  7.5   83  *   84  STO 05 85  LASTX   86  -   87  2   88  /   89  RCL 04         90  +   91  *   92  4   93  +   94  RCL 04   95  -   96  6   97  RCL 00   98  ST/ Y   99  + 100  RCL 02 101  * 102  STO Z 103  - 104  RCL 01 105  * 106  + 107  RCL 02 108  X^2 109  RCL 05 110  * 111  - 112  596.514 113  STO 04 114  ST/ 06 115  / 116  RCL 01  117  - 118  RCL 02       119  3 120  * 121  RCL 00 122  / 123  - 124  RCL 04 125  / 126  * 127  + 128  STO 03 129  RCL 00 130  RCL 07 131  ST* Y 132  - 133  STO 04 134  RCL 00 135  RCL 08 136  ST* Y 137  + 138  RCL 06 139  RCL 09 140  SIN 141  * 142  R-D 143  ST* 04 144  * 145  RCL 12  146  RCL 13        147  R-P 148  RDN 149  - 150  STO 01 151  RCL 10 152  RCL 11 153  R-P 154  CLX 155  RCL 04 156  - 157  ST+ 01 158  - 159  HMS 160  RCL 01 161  HMS 162  RCL 03 163  12756.274 164  * 165  END

( 214 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) Az2 ( ° ' " ) Z b1 ( ° ' " ) Az2 ( ° ' " ) Y L2 ( ° ' " ) Az1 ( ° ' " ) X b2 ( ° ' " ) D ( km )

Example:

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD"  >>>>    D   =   6181.621809 km                                ---Execution time = 11s---
RDN    Az1 =  51°47'36"76 = forward azimuth
RDN    Az2 = -68°09'59"26 = back azimuth

b-4) Program#4

-Second order formulae - given in reference [4] - may also be used to calculate the azimuths
-They employed the parametric latitudes and the program is longer and slower !

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

 01  LBL "TGD"    02  DEG   03  HR   04  STO 04    05  X<> T   06  HMS-   07  HR   08  STO 12   09  2   10  /   11  STO 00    12  RDN   13  HR   14  1   15  P-R   16  6378.137   17  STO 01   18  SIGN   19  298.257   20  1/X   21  STO 02   22  -   23  STO 05   24  /   25  R-P   26  X<>Y   27  STO 03   28  X<> 04   29  1   30  P-R   31  RCL 05   32  /   33  R-P   34  X<>Y   35  ST+ 03 36  RCL 04          37  -   38  2   39  ST/ 03   40  /   41  STO 04    42  SIN   43  X^2   44  ENTER^   45  SIGN   46  STO 06    47  LASTX   48  -   49  STO 05   50  RCL 03   51  SIN   52  X^2   53  ST* 05   54  ST- 06   55  -   56  RCL 00   57  SIN   58  X^2   59  *   60  X<>Y   61  ST* 06   62  +   63  ST/ 06   64  STO 07   65  1   66  X<>Y   67  -   68  ST/ 05   69  RCL 05   70  RCL 06 71  ST+ 05   72  -   73  2   74  ST* 05   75  *   76  STO 06          77  RCL 07   78  SQRT   79  ASIN   80  ST+ X   81  STO 08    82  D-R   83  LASTX   84  SIN   85  ST* 01   86  /   87  STO 09   88  ST+ X   89  X^2   90  STO 10   91  RCL 05   92  2   93  -   94  *   95  RCL 08   96  COS   97  STO 11   98  RCL 06   99  ST+ X 100  * 101  - 102  RCL 06 103  * 104  1 105  RCL 10 106  - 107  RCL 11 108  * 109  RCL 09 110  + 111  RCL 05        112  * 113  RCL 10 114  RCL 11  115  ST+ X 116  * 117  + 118  RCL 05 119  * 120  + 121  RCL 02 122  * 123  16 124  / 125  RCL 06 126  + 127  RCL 05 128  RCL 09 129  * 130  - 131  RCL 02 132  * 133  4 134  / 135  RCL 09 136  + 137  ST* 01 138  4 139  RCL 05 140  - 141  RCL 11 142  * 143  RCL 06 144  - 145  RCL 10  146  RCL 11  147  * 148  RCL 09        149  10 150  * 151  - 152  RCL 05 153  * 154  RCL 10 155  2 156  + 157  RCL 06 158  * 159  - 160  16 161  / 162  RCL 09 163  + 164  RCL 02 165  * 166  RCL 09 167  + 168  RCL 02 169  * 170  * 171  RCL 12 172  1 173  P-R 174  8 175  * 176  RDN 177  * 178  R^ 179  R-P 180  CLX 181  RCL 00         182  + 183  1 184  P-R 185  X<>Y 186  RCL 03  187  X<>Y 188  P-R 189  RCL 04 190  CHS 191  R^ 192  P-R 193  X<> Z 194  R-P 195  R^ 196  R^ 197  X<>Y 198  R-P 199  RDN 200  ENTER^ 201  R^ 202  ST+ Z 203  X<>Y 204  - 205  HMS 206  X<>Y 207  HMS 208  RCL 01 209  END

( 252 bytes / SIZE 013 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) Az2 ( ° ' " ) Y L2 ( ° ' " ) Az1 ( ° ' " ) X b2 ( ° ' " ) D ( km )

Example:

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD"  >>>>    D   =   6181.621790 km                                ---Execution time = 17s---
RDN    Az1 =  51°47'36"81 = forward azimuth
RDN    Az2 = -68°09'58"96 = back azimuth

Warning:

-The formulae that compute the azimuths do not work well if the difference in longitudes is very close to +/-90° !!!
-The variant hereunder doesn't have this defect.

b-5) Program#5

-The formulae ( 2nd order ) may be found in references [6] & [7]

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

 01 LBL "TGD"  02 DEG  03 HR  04 STO 06  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 07  17 RDN  18 HR  19 298.257  20 1/X  21 STO 04  22 SIGN  23 LASTX  24 -  25 STO 03  26 P-R  27 LASTX  28 /  29 R-P  30 X<>Y  31 STO 05  32 1  33 P-R  34 STO 01  35 X<>Y  36 STO 02  37 RCL 06  38 RCL 03  39 P-R  40 LASTX 41 /  42 R-P  43 X<>Y  44 STO 06        45 1  46 P-R  47 STO 08  48 X<>Y  49 STO 09  50 RCL 02  51 *  52 STO 03  53 RCL 01  54 ST* 09  55 RCL 08  56 ST* 02  57 *  58 STO 10  59 RCL 07  60 2  61 /  62 SIN  63 X^2  64 *  65 RCL 06  66 RCL 05  67 -  68 2  69 /  70 SIN  71 X^2  72 +  73 SQRT  74 ASIN  75 ST+ X  76 D-R  77 STO 13  78 LASTX  79 1  80 P-R 81 STO 11  82 X<>Y  83 STO 12  84 RCL 07        85 SIN  86 RCL 10  87 *  88 X<>Y  89 /  90 STO 06  91 X^2  92 SIGN  93 LASTX  94 -  95 STO 05  96 RCL 13  97 X^2  98 RCL 11  99 * 100 RCL 12 101 / 102 STO 14 103 RCL 11 104 RCL 12 105 * 106 STO 15 107 RCL 13 108 5 109 * 110 - 111 4 112 / 113 + 114 RCL 05 115 * 116 RCL 13 117 X^2 118 RCL 12 119 / 120 STO 10 121 LASTX 122 2 123 / 124 + 125 RCL 03       126 * 127 - 128 RCL 04 129 X^2 130 ST* Y 131 LASTX 132 + 133 STO 00 134 RCL 13 135 * 136 + 137 RCL 06 138 * 139 R-D 140 ST+ 07 141 RCL 10 142 RCL 11 143 RCL 15 144 * 145 + 146 RCL 03 147 * 148 RCL 05 149 * 150 RCL 14 151 RCL 11 152 X^2 153 RCL 15 154 * 155 4 156 / 157 + 158 RCL 15 159 RCL 13 160 + 161 8 162 / 163 - 164 RCL 05       165 X^2 166 * 167 - 168 RCL 15 169 RCL 03 170 X^2 171 * 172 - 173 RCL 14 174 RCL 05 175 * 176 + 177 RCL 04 178 X^2 179 * 180 RCL 13 181 RCL 15 182 + 183 RCL 00 184 * 185 RCL 05 186 * 187 - 188 2 189 / 190 RCL 10 191 RCL 04 192 X^2 193 * 194 2 195 / 196 RCL 12 197 RCL 00 198 * 199 - 200 RCL 03 201 * 202 - 203 1 204 RCL 04       205 - 206 * 207 RCL 13 208 + 209 STO 00 210 RCL 07 211 SIN 212 STO 11 213 RCL 01 214 * 215 CHS 216 RCL 02 217 RCL 09 218 RCL 07 219 COS 220 STO 12 221 * 222 - 223 R-P 224 X<>Y 225 HMS 226 RCL 11 227 RCL 08 228 * 229 RCL 02 230 RCL 12 231 * 232 RCL 09 233 X<>Y 234 - 235 R-P 236 RDN 237 HMS 238 RCL 00 239 6378.137 240 * 241 END

( 271 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) Az2 ( ° ' " ) Y L2 ( ° ' " ) Az1 ( ° ' " ) X b2 ( ° ' " ) D ( km )

Example:

-77.0356    ENTER^
38.55172  ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD"  >>>>    D   =   6181.621809 km                                ---Execution time = 17s---
RDN    Az1 =  51°47'36"81 = forward azimuth
RDN    Az2 = -68°09'58"96 = back azimuth

c) Forward Geodesic Problem

-Given the coordinates ( L1 , b1 ) of the first point  P1 , the forward azimuth Az1 and the length D of the geodesic,
"FWD" computes the longitude, the latitude ( L2 , b2 ) of the second point  P2 and the backward azimuth Az2

-The formulas are given in the reference [4]

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

 01  LBL "FWD"   02  DEG   03  R-D   04  STO 01    05  X<>Y   06  HR   07  STO 02    08  R^   09  HR   10  STO 04    11  R^   12  HR   13  1   14  P-R   15  LASTX   16  298.257   17  1/X   18  STO 05    19  -   20  STO 00   21  /   22  R-P   23  SIGN   24  P-R   25  STO 13   26  X<>Y   27  STO 03   28  RCL 02   29  RCL 13   30  P-R   31  STO 07   32  X<>Y   33  STO 06   34  X^2 35  SIGN   36  LASTX   37  -   38  RCL 05    39  *   40  4   41  /   42  ST0 08         43  1   44  -   45  RCL 06    46  RCL 05   47  RCL 06   48  *   49  STO 09   50  *   51  STO 10    52  +   53  *   54  ST/ 01   55  RCL 10   56  2   57  +   58  RCL 08   59  *   60  X<>Y   61  /   62  STO 11   63  RCL 03   64  RCL 13   65  RCL 02   66  COS   67  *   68  R-P 69  CLX   70  ACOS   71  X<>Y   72  -   73  STO 12   74  RCL 01    75  6378.137   76  /   77  STO 01         78  -   79  ST+ X   80  ENTER^   81  COS   82  RCL 11    83  *   84  1   85  X<>Y   86  -   87  X<>Y   88  RCL 01    89  +   90  COS   91  ENTER^   92  X^2   93  .5   94  -   95  RCL 08   96  X^2   97  *   98  RCL 01   99  ST+ X 100  SIN 101  * 102  X<> Z 103  * 104  RCL 11 105  * 106  RCL 01  107  SIN 108  * 109  - 110  R-D 111  RCL 01       112  + 113  STO 10  114  RCL 12  115  ST+ X 116  X<>Y 117  - 118  STO 11 119  RCL 06  120  CHS 121  RCL 03 122  LASTX 123  SIN 124  STO 12 125  * 126  RCL 07 127  RCL 10 128  COS 129  STO 01 130  * 131  - 132  R-P 133  X<>Y 134  STO 05 135  RCL 01 136  RCL 03 137  * 138  RCL 07 139  RCL 12  140  * 141  + 142  RCL 07 143  RCL 01 144  * 145  RCL 03       146  RCL 12  147  * 148  - 149  X^2 150  RCL 06 151  X^2 152  + 153  SQRT 154  RCL 00  155  * 156  R-P 157  X<>Y 158  X<> 02 159  RCL 12 160  P-R 161  RCL 03 162  * 163  RCL 01 164  RCL 13 165  * 166  X<>Y 167  - 168  R-P 169  CLX 170  RCL 08 171  RCL 09 172  ST* Y 173  - 174  RCL 10  175  * 176  + 177  RCL 08  178  RCL 09       179  * 180  RCL 12 181  * 182  RCL 11 183  COS 184  * 185  + 186  ST+ 04 187  RCL 05  188  HMS 189  RCL 02 190  HMS 191  RCL 04 192  360 193  MOD 194  PI 195  R-D 196  X>Y? 197  CLX 198  ST+ X 199  - 200  HMS 201  END

( 234 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) Az2 ( ° ' " ) Z b1 ( ° ' " ) Az2 ( ° ' " ) Y Az1 ( ° ' " ) b2 ( ° ' " ) X D ( km ) L2 ( ° ' " )

where   Az1  is the forward azimuth in the direction from the first point P1 towards the second point P2

Example:     L1 = 10°30'     b1 = 49°41'     Az1 = 12°24'     D = 16000 km

10.30   ENTER^
49.41   ENTER^
12.24   ENTER^
16000   XEQ "FWD"  >>>>     L2 = -177°03'08"01                ---Execution time = 16s---
RDN      b2 =   -14°06'40"72
RDN    Az2 =    -8°15'03"68

Notes:

-Andoyer's formulas have one advantage: they are non-iterative!
-In most cases, their accuracy is satisfactory, all the more that the Earth's irregularities cannot be taken into account.
-If, however, a better precision is required, this can be achieved via Vincenty's formulas below.

2°) Vincenty's Formulae

a) Geodesic Distance and Azimuths

-Let   L =  L2 - L1  = L0   ;   U = Arctan((1-f).tan b1)  ;   V = Arctan((1-f).tan b2)

-We compute recursively:      Ln+1 = L + (1-C).f.(sin µ).{ c + C.sin c [ cos 2d + C.cos c ( -1 + 2.cos2 2d ) ] }

where      sin c = [ ( cos V sin Ln )2 + ( cos U sin V - sin U cos V cos Ln )2 ]1/2
cos c = sin U sin V + cos U cos V cos Ln
µ = Arcsin [ ( cos U cos V sin Ln )/sin c ]  ;  µ = azimuth of the geodesic at the equator
C = (f/16).(cos2 µ ).[ 4 + f.(4-3.cos2 µ ) ]
cos 2d = cos c - 2.sin U sin V / cos2 µ

until   Ln+1 = Ln = lambda

-Then,     D = a.A.(1-f).(c-D)

with     A = 1 + (u/16384).{ 4096 + u.[-768 + u.(320-175.u)] }      ;       B = (u/1024).{ 256 + u.[ -128 + u.(74-47.u)] }
D = (B.sin c).{ cos 2d + (B/4). [ ( cos c ).( -1 + 2.cos2 2d ) - (B/6).(cos 2d).(-3+4.sin2 c ).(-3+4.cos2 2d )] }
u =  f(2-f)/(1-f)2 . cos2 µ

-Actually, the HP-41 works with 10 digits only. So several terms may be omitted without reducing the accuracy,  namely the terms in  u4 and  B3
-The azimuths  Az1 &  Az2 are finally calculated by the formulas:

Tan Az1 = [ cos V sin (lambda) ] / [ cos U sin V - sin U cos V cos (lambda) ]
Tan Az2 = [ -cos U sin (lambda) ] / [ sin U cos V - cos U sin V cos (lambda) ]

Data Registers:  R00 = f  ;  R01 = U  ;  R02 = V  ;  R04 = lambda  ;  R06 = µ = Azequator  ;   R03 , R05 and R07 to R09: temp
Flags:  /
Subroutines: /

 01  LBL "TGD2"   02  DEG   03  HR   04  X<> T   05  HMS-   06  HR   07  360   08  MOD   09  PI   10  R-D   11  X>Y?   12  CLX   13  ST+ X   14  -   15  STO 03    16  STO 04    17  RDN   18  HR   19  298.257   20  1/X   21  STO 00    22  SIGN   23  LASTX   24  -   25  STO 09   26  P-R   27  LASTX   28  /   29  R-P   30  RDN   31  STO 01   32  CLX   33  RCL 09   34  P-R   35  LASTX   36  /   37  R-P   38  X<>Y   39  STO 02 40  LBL 01   41  VIEW 04    42  RCL 02          43  COS   44  STO 05   45  STO 06   46  RCL 04   47  SIN   48  *   49  STO 08   50  RCL 01   51  COS   52  ST* 05   53  ST* 08   54  RCL 02    55  SIN   56  STO 07   57  *   58  RCL 01   59  SIN   60  ST* 07   61  RCL 06    62  *   63  RCL 04   64  COS   65  ST* 05   66  *   67  -   68  R-P   69  RCL 05   70  RCL 07   71  +   72  R-P   73  X<> 08   74  X<>Y   75  STO 05   76  SIN   77  X#0?   78  / 79  ASIN   80  STO 06    81  COS   82  X^2   83  RCL 05          84  COS   85  RCL 07   86  ENTER^   87  +   88  R^   89  X=0?   90  SIGN   91  /   92  -   93  STO 07    94  CLX   95  3   96  *   97  4   98  X<>Y   99  - 100  RCL 00  101  * 102  4 103  + 104  * 105  RCL 00 106  * 107  16 108  / 109  RCL 05 110  SIN 111  RCL 07 112  RCL Z 113  * 114  * 115  R-D 116  RCL 05 117  + 118  RCL 06 119  SIN 120  * 121  RCL 00        122  * 123  1 124  R^ 125  - 126  * 127  RCL 03 128  + 129  ENTER^ 130  X<> 04 131  - 132  ABS 133   E-7 134  XY 221  HMS 222  RCL 03 223  RCL 05 224  RCL 08 225  - 226  R-P 227  RDN 228  HMS 229  RCL 09 230  6378.137 231  * 232  CLD 233  END

( 289 bytes / SIZE 010 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) Az2 ( ° ' " ) Z b1 ( ° ' " ) Az2 ( ° ' " ) Y L2 ( ° ' " ) Az1 ( ° ' " ) X b2 ( ° ' " ) D ( km )

where   Az1  =  forward azimuth  ( from P1 to P2 )
and     Az2  =  back azimuth  ( from P2 to P1 )

Example:

-77.0356    ENTER^
38.55172   ENTER^
2.20138   ENTER^
48.50112  XEQ "TGD2"  >>>>   the successive Ln-values are displayed  and  ( after 47 seconds )

D   =  6181.621787 km
RDN    Az1 =  51°47'36"8132
RDN    Az2 = -68°09'58"9656

-We also have  R06 = µ = azimuth at the equator = 37.74571879°

Notes:

-The accuracy is excellent, provided the 2 points are not nearly antipodal.
-Otherwise, the algorithm doesn't converge! It happens when | Ln | exceeds 180° ( in register R04 )
-If you don't want to compute the azimuths, replace lines 193 to 229 by  RCL 09  *

b) Forward Geodesic Problem

-Given the coordinates ( L1 , b1 ) of the first point  P1 , the forward azimuth Az1 and the length D of the geodesic,
"FWD" computes the longitude and the latitude ( L2 , b2 ) of the second point  P2

-The formulas are given in the reference [3]  ( the terms in  u4 and  B3 have been neglected )

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

 01  LBL "FWD"   02  DEG   03  R-D   04  STO 01    05  RDN   06  HR   07  STO 02    08  RDN   09  HR   10  X<>Y   11  HR   12  STO 03   13  CLX   14  SIGN   15  P-R   16  LASTX   17  298.257   18  1/X   19  STO 00    20  -   21  STO 09   22  ST/ 01   23  /   24  R-P   25  X<>Y   26  STO 04   27  1   28  P-R   29  STO 13   30  RCL 02   31  COS   32  *   33  R-P   34  X<>Y   35  STO 05   36  6378.137 37  ST/ 01   38  RCL 13   39  RCL 02    40  SIN   41  *   42  STO 06         43  ASIN   44  2   45  RCL 00   46  ST- Y   47  *   48  X<>Y   49  COS   50  STO 07   51  RCL 09    52  /   53  X^2   54  *   55  12   56  RCL Y   57  5   58  *   59  X<>Y   60  -   61  *   62  64   63  +   64  *   65  256   66  ST+ Y   67  /   68  ST/ 01   69  CLX   70  37   71  *   72  64 73  -   74  *   75  512   76  /   77  4   78  1/X   79  +   80  *   81  STO 08         82  RCL 01   83  STO 10   84  LBL 01   85  VIEW 10   86  RCL 10   87  RCL 05   88  ST+ X   89  +   90  COS   91  STO 11    92  X^2   93  ST+ X   94  RCL 10   95  COS   96  ST* Y   97  -   98  STO 12   99  RCL 08 100  * 101  4 102  / 103  RCL 11 104  + 105  RCL 10 106  SIN 107  * 108  RCL 08 109  * 110  R-D 111  RCL 01  112  + 113  ENTER^ 114  X<> 10 115  - 116  ABS 117   E-7 118  XY 154  HMS 155  RCL 02  156  RCL 05       157  P-R 158  RCL 04 159  * 160  RCL 01 161  X<>Y 162  - 163  R-P 164  CLX 165  RCL 07  166  X^2 167  STO 07 168  3 169  * 170  4 171  - 172  RCL 00 173  ST* 06 174  ST* 07 175  * 176  4 177  - 178  RCL 07 179  * 180  16 181  / 182  STO 08 183  ST* 05 184  RCL 12  185  * 186  RCL 11 187  - 188  RCL 05  189  * 190  R-D 191  RCL 10       192  + 193  RCL 06 194  * 195  ST- Y 196  RCL 08  197  * 198  - 199  RCL 03 200  + 201  360  202  MOD 203  PI 204  R-D 205  X>Y? 206  CLX 207  ST+ X 208  - 209  HMS 210  CLD 211  END

( 263 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) b2 ( ° ' " ) Z b1 ( ° ' " ) b2 ( ° ' " ) Y Az1 ( ° ' " ) b2 ( ° ' " ) X D ( km ) L2 ( ° ' " )

where   Az1  is the forward azimuth in the direction from the first point P1 towards the second point P2

Example:     L1 = 10°30'     b1 = 49°41'     Az1 = 12°24'     D = 16000 km       Calculate  L2 & b2

10.30   ENTER^
49.41   ENTER^
12.24   ENTER^
16000   XEQ "FWD"  >>>  the successive sigma-approximations are displayed and eventually, it yields:

L2 = -177°03'07"987     X<>Y    b2 =  -14°06'40"7154       ( execution time = 25 seconds )

-This method works for antipodal points too and - more generally - for any length D !

3°) Numerical Integration

-The rigorous formulas are:

D/a = s = §r1r2  r [ ( 1-e2.r2 )/( 1 - r2 )/ ( r2 - k2 ) ] 1/2 dr                                                             ( § = Integral )

where  e = f.(2-f) = the eccentricity of the ellipsoid,
(a.r)  is the radius of the parallel of latitude b ,   r = ( cos b ).( 1 - e2.sin2 b ) -1/2

and  k  is the solution of    | L1 - L2 | =  §r1r2  (k/r) [ ( 1-e2.r2 )/( 1 - r2 )/ ( r2 - k2 ) ] 1/2 dr

-Actually,  a given geodesic cuts all the parallels at an angle phi such that   r.cos(phi) = k = Constant
-In order to remove the singularities of these integrals, I've made the following change of argument:

sin 2µ = (2/( 1 - k2 )).[ ( 1 - r2 ).( r2 - k2 ) ] 1/2    or   sin µ = +/- [ ( 1 - r2 )/( 1 - k2) ] 1/2   ;    cos µ = +/- [ ( r2 - k2 )/( 1 - k2 ) ] 1/2

-And the integrals become:

s = §µ1µ2  ( 1 - e2.( cos2 µ + k2.sin2 µ ) )1/2 dµ            ( E )

| L1 - L2 | =  §µ1µ2  ( k/( cos2 µ + k2.sin2 µ ) ).( 1 - e2.( cos2 µ + k2.sin2 µ ) )1/2

-However, this last integral is difficult if the geodesic passes near the Pole, so I've split it into 2 parts which eventually yields:

| L1 - L2 | =  [ Arctan ( k.tan µ ) ]µ1µ2 -   §µ1µ2  k.e2/[ 1 + ( 1 - e2.( cos2 µ + k2.sin2 µ ) )1/2 ] dµ     ( E' )

-"TGD3"  solves equation (E') by the secant method, thus giving k , µ1 , µ2
-And Gaussian integration produces  s   through equation (E).

Data Registers:  R00 thru R08 are used by "GL3"    R00 = "S" ;  R01 = µ1 ;  R02 = µ2
R03 = 2  for solving (E')      ( line 64 )              ( number of subintervals over which
R03 = 4  for evaluating  s    ( line 191 )               the Gauss-Legendre formula is applied )

These values produce ( almost ) 10 significant figures for the Earth.
They should be increased if you want to use this program for another planet with a greater flattening.

R09 = -e2  ;  R10 = r1 ( or -r1 if µ2 > 90° )  ;  R11 = sign(b1).(1-r1)1/2  ;  R12 = r2 ( calling r2 < r1 ) ;  R13 = (1-r2)1/2
R14 & R15 = the successive k-values  ;  R16 = f(k) ;  R17 = | L1 - L2 |
Flag:  F10
Subroutine:  "GL3"  Gauss-Legendre 3-point formula ( cf "Numerical Integration for the HP-41" )

-Line 32 =  e2  with  f = 1/298.257
-Line 71:    if  | L1 - L2 | < 179°23'47"377 ,  s = | L1 - L2 | radians
-Lines 195 to 198 may look strange, but they are useful to correct the inaccuracy in the µ-values, especially if  k is close to 1.

 01  LBL "TGD3"   02  DEG   03  HR   04  X<> T   05  HMS-   06  HR   07  ABS   08  PI   09  R-D   10  X>Y?   11  CLX   12  ST+ X   13  -   14  ABS   15  STO 17    16  RDN   17  HR   18  ST* Z   19  ABS   20  X<>Y   21  ABS   22  X>Y?   23  X<>Y   24  RCL Z   25  SIGN   26  *   27  COS   28  LASTX   29  SIN   30  STO 11   31  X^2   32  .006694385   33  CHS   34  STO 09   35  * 36  1   37  +   38  SQRT   39  ST/ 11   40  /   41  STO 10           42  X<>Y   43  COS   44  LASTX   45  SIN   46  STO 13    47  X^2   48  RCL 09    49  *   50  1   51  +   52  SQRT   53  ST/ 13   54  /   55  STO 12   56  STO 14   57  STO 15   58  SIGN   59  RCL 09   60  +   61  SQRT   62  ST* 11   63  ST* 13   64  2   65  STO 03   66  "S"   67  ASTO 00   68  RCL 13   69  X#0?   70  GTO 00 71  179.3964936   72  RCL 17   73  X<=Y?   74  GTO 04   75  LBL 00   76  SF 10   77  XEQ 01   78  STO 16   79  SIGN   80  ST* 10   81  .9   82  ST* 15   83  GTO 03   84  LBL 01   85  XEQ 02   86  RCL 15    87  *   88  RCL 09   89  *   90  RCL 02   91  1   92  P-R   93  RCL 15   94  ST* Z   95  RDN   96  R-P   97  RDN   98  +   99  RCL 01 100  1 101  P-R 102  RCL 15 103  ST* Z 104  RDN 105  R-P 106  RDN 107  - 108  RCL 17 109  - 110  RTN 111  LBL 02  112  RCL 13         113  RCL 10 114  SIGN 115  RCL 12 116  X^2 117  RCL 15  118  X^2 119  - 120  SQRT 121  * 122  R-P 123  X<>Y 124  STO 02 125  RCL 11 126  RCL 10 127  X^2 128  RCL 15 129  X^2 130  - 131  SQRT 132  R-P 133  X<>Y 134  STO 01 135  XEQ "GL3" 136  RTN 137  LBL "S" 138  1 139  P-R 140  X^2 141  X<>Y 142  RCL 15 143  * 144  X^2 145  + 146  RCL 09         147  * 148  1 149  + 150  SQRT 151  FC? 10 152  RTN 153  1 154  + 155  1/X 156  RTN 157  LBL 03 158  CLX 159  SIGN 160   E-10 161  - 162  RCL 15  163  ABS 164  X>Y? 165  X<>Y 166  STO 15 167  RCL 12 168  X 16 175  - 176  X#0? 177  / 178  RCL 14         179  RCL 15  180  STO 14 181  - 182  * 183  ST+ 15 184  ABS 185   E-9 186  X

( 295 bytes / SIZE 018 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) / Y L2 ( ° ' " ) / X b2 ( ° ' " ) D ( km )

Example1:

77.0356    ENTER^
38.55172  ENTER^
-2.20138   ENTER^
48.50112  XEQ "TGD3"  >>>>   the successive k-values are displayed  and eventually,    D = 6181.621797 km   ( in 2m45s ) , the last digit should be a 4.

-We also have in register R15    k = 0.612158197

Example2:      L1 = b1 = 0  ;   L2 = 179°51'   b2 = 0

0      ENTER^  ENTER^
179.51  ENTER^
0      XEQ "TGD3"  >>>>   D = 20001.85464 km ( in 2m10s ) , the last digit should be a 3.
and  R15 = k = 0.248743230

-Unlike  TGD & TGD2 , TGD3 produces exact results for nearly antipodal points too ... But it is much slower!

-If the 2 points are on the equator and  | L1 - L2 |  <  179°23'47"377  we have  k = 1  whereas  k = 0  if  | L1 - L2 | = 180°
-The maximum geodesic distance occurs with ( 0 ; 0 ) ( 180 ; 0 )  and  Dmax = 6378.137 * 3.136328278 = 20003.93143 km

4°)  Series Expansions

a)  Program#1

-"TGD4" employs the same integrals as "TGD3", but they are evaluated by a series expansion: the program runs faster and no subroutine is used.
-The first neglected term is proportinal to e6 :

Formulae:

| L1 - L2 | =  [ Arctan ( k.tan µ ) ]µ1µ2 - k.(e2/2).( µ2 - µ1 ) - k.(e4/8).[ ( µ2 - µ1 ) - ( 1 - k2).( µ2 - µ1 )/2 + (1/4).( 1 - k2).( Sin 2µ2 - Sin 2µ1 ) ]

s / a  = ( µ2 - µ1 ) - (e2/2).[ ( µ2 - µ1 ) - ( 1 - k2).( µ2 - µ1 )/2 + ( 1 - k2).( Sin 2µ2 - Sin 2µ1 )/4 ]
- (e4/8).[ ( µ2 - µ1 ) - ( 1 - k2).( µ2 - µ1 ) + ( 1 - k2).( Sin 2µ2 - Sin 2µ1 )/2 + (3/8).( 1 - k2)2.( µ2 - µ1 )
- ( 1 - k2)2.( Sin 2µ2 - Sin 2µ1 )/4 + ( 1 - k2)2.( Sin 4µ2 - Sin 4µ1 )/32 ]

Data Registers:  R00 = | L1 - L2 | ;  R01 = µ1 ;  R02 = µ2 ;  R03 & R04 = the successive k-values ; R05 = f(k)
R06 = r1 ( or -r1 if µ2 > 90° )  ;  R07 = sign(b1).(1-r1)1/2  ;  R08 = r2 ( calling r2 < r1 ) ;  R09 = (1-r2)1/2
R10 = -e2 ;  R11 = Sin 2µ2 - Sin 2µ1 ;  R12 = µ2 - µ1 ;  R13 = 1 - k2  ;  R14 = temp ;  R15 = 10 -9
Flags:  /
Subroutines:  /

-Line 32 may be replaced by  .38356  D-R  thus saving 3 bytes.

 01  LBL "TGD4"   02  DEG   03  HR   04  X<> T   05  HMS-   06  HR   07  ABS   08  PI   09  R-D   10  X>Y?   11  CLX   12  ST+ X   13  -   14  ABS   15  STO 00    16  RDN   17  HR   18  ST* Z   19  ABS   20  X<>Y   21  ABS   22  X>Y?   23  X<>Y   24  RCL Z   25  SIGN   26  *   27  1   28  P-R   29  X<>Y   30  STO 07    31  X^2   32  .006694385   33  CHS   34  STO 10   35  *   36  1   37  +   38  SQRT   39  ST/ 07   40  / 41  STO 06   42  SIGN   43  P-R   44  X<>Y   45  STO 09   46  X^2   47  RCL 10          48  *   49  1   50  +   51  SQRT   52  ST/ 09   53  /   54  STO 03    55  STO 04    56  STO 08   57  SIGN   58  RCL 10   59  +   60  SQRT   61  ST* 07   62  ST* 09   63   E-9   64  STO 15   65  XEQ 01   66  STO 05   67  SIGN   68  ST* 06   69  .9   70  ST* 04   71  GTO 03   72  LBL 01   73  XEQ 02   74  RCL 02    75  ST+ X   76  SIN   77  RCL 01   78  ST+ X   79  SIN   80  - 81  R-D   82  STO 11   83  2   84  /   85  -   86  1   87  RCL 04          88  X^2   89  -   90  STO 13   91  *   92  RCL 12   93  ST+ X   94  -   95  RCL 10    96  *   97  STO 14   98  8   99  / 100  RCL 12 101  + 102  RCL 04 103  * 104  RCL 10 105  * 106  2 107  / 108  RCL 02  109  1 110  P-R 111  RCL 04  112  ST* Z 113  RDN 114  R-P 115  RDN 116  + 117  RCL 01 118  1 119  P-R 120  RCL 04 121  ST* Z 122  RDN 123  R-P 124  RDN 125  - 126  RCL 00       127  - 128  RTN 129  LBL 02 130  RCL 09 131  RCL 06 132  SIGN 133  RCL 08  134  X^2 135  RCL 04  136  X^2 137  - 138  SQRT 139  * 140  R-P 141  X<>Y 142  STO 02 143  RCL 07 144  RCL 06 145  X^2 146  RCL 04  147  X^2 148  - 149  SQRT 150  R-P 151  RDN 152  STO 01  153  - 154  STO 12 155  RTN 156  LBL 03 157  SIGN 158  RCL 15 159  - 160  RCL 04 161  ABS 162  X>Y? 163  X<>Y 164  STO 04 165  RCL 08  166  X 05 173  - 174  X#0? 175  / 176  RCL 03  177  RCL 04 178  STO 03 179  - 180  * 181  ST+ 04 182  ABS 183  RCL 15  184  X

( 293 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) / Y L2 ( ° ' " ) / X b2 ( ° ' " ) D ( km )

Example1:

77.0356    ENTER^
38.55172  ENTER^
-2.20138   ENTER^
48.50112  XEQ "TGD4"  >>>>   the successive k-values are displayed  and eventually,    D = 6181.621792 km   ( in 64s ) , the last digit should be a 4.

-We also have in register R04    k = 0.612158197

Example2:      L1 = b1 = 0  ;   L2 = 179°51'   b2 = 0

0      ENTER^  ENTER^
179.51  ENTER^
0      XEQ "TGD4"  >>>>   D = 20001.85474 km ( in 40s ) , error = 11cm.
and  R04 = k = 0.248743896

-The maximum error is about 13 centimeters.

b)  Program#2  ( better precision )

-The first neglected term is proportinal to e8

Formulae:

| L1 - L2 | =  [ Arctan ( k.tan µ ) - k e2 { µ [ 1/2 + e2 ( 1 + k2 ) / 16 + e4 ( 3 + 2 k2 + 3 k4 ) / 128 ]
+ ( Sin 2.µ ) [ e2 ( 1 - k2 ) / 32 + e4 ( 1 - k4 ) / 64 ]
+ ( Sin 4.µ ) [ e4 ( 1 - 2 k2 + k4 ) / 512 ] }  ]µ1µ2

s / a     =  {  µ  [ 1 - e2 ( 1 + k2 ) / 4 - e4 ( 3 + 2 k2 + 3 k4 ) / 64 - e6 ( 5 + 3 k2 + 3 k4 + 5 k6 ) / 256 ]
+ ( Sin 2.µ ) [ -e2 ( 1 - k2 ) / 8 - e4 ( 1 - k4 ) / 32 - e6 ( 15 + 3 k2 - 3 k4 - 15 k6 ) / 1024 ]
+ ( Sin 4.µ ) [ - e4 ( 1- 2 k2 + k4 ) / 256 - 3 e6 ( 1 - k2 - k4 + k6 ) / 1024 ]
+ ( Sin 6.µ ) [ -e6 ( 1 - 3 k2 + 3 k4 - k6 ) / 3072 }µ1µ2

Data Registers:  R00 = | L1 - L2 | ;  R01 = µ1 ;  R02 = µ2 ;  R03 & R04 = the successive k-values ; R05 = f(k)
R06 = r1 ( or -r1 if µ2 > 90° )  ;  R07 = sign(b1).(1-r1)1/2  ;  R08 = r2 ( calling r2 < r1 ) ;  R09 = (1-r2)1/2
R10 = e2 ;  R11 = Sin 2µ2 - Sin 2µ1 ;  R12 = µ2 - µ1 ;  R13 = 1 - k2  ;  R14 = k2 ;  R15 = 10 -9
Flags:  /
Subroutines:  /

-Line 72  is a three-byte GTO 02

 01 LBL "TGD5"  02 DEG  03 HR  04 X<> T  05 HMS-  06 HR  07 ABS  08 PI  09 R-D  10 X>Y?  11 CLX  12 ST+ X  13 -  14 ABS  15 STO 00  16 RDN  17 HR  18 ST* Z  19 ABS  20 X<>Y  21 ABS  22 X>Y?  23 X<>Y  24 RCL Z  25 SIGN  26 *  27 1  28 P-R  29 X<>Y  30 STO 07  31 X^2  32 .006694385  33 STO 10  34 CHS  35 *  36 1  37 +  38 SQRT  39 ST/ 07  40 /  41 STO 06            42 SIGN  43 P-R  44 X<>Y  45 STO 09  46 X^2  47 RCL 10  48 *  49 CHS  50 1 51 +  52 SQRT  53 ST/ 09  54 /  55 STO 03  56 STO 04  57 STO 08  58 SIGN  59 RCL 10  60 -  61 SQRT  62 ST* 07  63 ST* 09  64  E-9  65 STO 15  66 XEQ 01  67 STO 05  68 SIGN  69 ST* 06  70 .9  71 ST* 04  72 GTO 02  73 LBL 01  74 RCL 09  75 RCL 06  76 SIGN  77 RCL 08  78 X^2  79 RCL 04  80 X^2  81 -  82 SQRT  83 *  84 R-P  85 X<>Y  86 STO 02  87 RCL 07  88 RCL 06  89 X^2  90 RCL 04            91 X^2  92 -  93 SQRT  94 R-P  95 RDN  96 STO 01  97 -  98 STO 12  99 LASTX 100 4 101 * 102 SIN 103 RCL 02 104 4 105 * 106 SIN 107 - 108 R-D 109 STO 16 110 1 111 RCL 04 112 X^2 113 STO 14 114 - 115 STO 13 116 RCL 10 117 * 118 X^2 119 * 120 512 121 / 122 1 123 RCL 14 124 X^2 125 - 126 RCL 10 127 * 128 RCL 13 129 ST+ X 130 + 131 RCL 10 132 * 133 64 134 / 135 RCL 02 136 ST+ X 137 SIN 138 RCL 01 139 ST+ X 140 SIN 141 - 142 R-D 143 STO 11           144 * 145 - 146 RCL 14 147 3 148 * 149 2 150 + 151 RCL 14 152 * 153 3 154 + 155 RCL 10 156 * 157 8 158 / 159 RCL 14 160 + 161 RCL 10 162 ST* Y 163 + 164 8 165 / 166 RCL 12 167 ST* Y 168 + 169 2 170 / 171 - 172 RCL 04 173 * 174 RCL 10 175 * 176 RCL 02 177 1 178 P-R 179 RCL 04 180 ST* Z 181 RDN 182 R-P 183 RDN 184 + 185 RCL 01 186 1 187 P-R 188 RCL 04 189 ST* Z 190 RDN 191 R-P 192 RDN 193 - 194 RCL 00           195 - 196 RTN 197 LBL 02 198 SIGN 199 RCL 15 200 - 201 RCL 04 202 ABS 203 X>Y? 204 X<>Y 205 STO 04 206 RCL 08 207 X 05 214 - 215 X#0? 216 / 217 RCL 03 218 RCL 04 219 STO 03 220 - 221 * 222 ST+ 04 223 ABS 224 RCL 15 225 X

( 429 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T L1 ( ° ' " ) / Z b1 ( ° ' " ) / Y L2 ( ° ' " ) / X b2 ( ° ' " ) D ( km )

Example1:

77.0356    ENTER^
38.55172  ENTER^
-2.20138   ENTER^
48.50112  XEQ "TGD5"  >>>>   the successive k-values are displayed  and eventually,    D = 6181.621797 km                ---Execution time = 87s---

-We also have in register R04    k = 0.612158197

Example2:      L1 = b1 = 0  ;   L2 = 179°51'   b2 = 0

0      ENTER^  ENTER^
179.51  ENTER^
0      XEQ "TGD5"  >>>>   D = 20001.85464 km                              ---Execution time = 57s---

and  R04 = k = 0.248743230

Notes:

-Roundoff errors may decrease the precision.
-But with free42, it yields:  6181.6217939  &  20001.8546318

5°) Euclidean Distance

a) Geocentric Coordinates

-If the latitude b and the height h of an observer O are known, the following routine calculates
the geocentric latitude b' and the distance rho = OC  where C is the center of the Earth.

Formulae:

(rho)  sin b' = a(1-f) sin u + h sin b                  a = 6378.137 km
(rho) cos b' =  a cos u + h cos b                      f = 1/298.257

where   tan u = (1-f) tan b

Data Registers: /
Flags: /
Subroutines: /

 01  LBL "GEOC" 02  DEG 03  HR 04  RCL Y  05  STO M  06  CLX 07  SIGN 08  P-R 09  ST* M  10  X<>Y 11  ST* Z  12  298.257         13  1/X 14  SIGN 15  ST- L 16  X<> L           17  CHS 18  ST* Y 19  X<> Z  20  R-P 21  CLX 22  6378.137 23  P-R 24  ST+ M          25  RDN 26  * 27  + 28  0 29  X<> M          30  R-P 31  X<>Y 32  HMS 33  END

( 65 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y h (km) rho (km) X b ( ° ' " ) b' ( ° ' " )

Example:   Find the geocentric coordinates of the Mount Palomar Observatory     b = 33°21'22"4     h = 1706 m = 1.706 km

1.706      ENTER^
33.21224  XEQ "GEOC"  >>>>    b'  = 33°10'47"12
X<>Y   rho = 6373.4156 km

Notes:

-If you prefer h & rho in meters, replace line 22 by 6378137
-In several applications, (rho) sin b' & (rho) cos b' are more useful than rho or b':  Simply delete lines 30-31-32

(rho) sin b'  will be inY-register,
(rho) cos b' will be in X-register.

b) 3D-Distance

-"EDST" computes the geocentric rectangular coordinates of 2 points and their Euclidean 3D-distance
knowing their longitudes, latitudes and heights.

Data Registers:               R00 is unused                                   ( Registers R01 thru R06 are to be initialized before executing "EDST" )

•  R01 = L1 ( ° ' " )   •  R04 = L2 ( ° ' " )        R07 = x1     R10 = x2
•  R02 = b1 ( ° ' " )   •  R05 = b2 ( ° ' " )        R08 = y1     R11 = y2         xi  yi  zi   in km
•  R03 = h1 ( km )    •  R06 = h2  ( km )        R09 = z1     R12 = z2
Flags: /
Subroutine:   "GEOC"

 01  LBL "EDST" 02  RCL 03 03  RCL 02  04  XEQ "GEOC" 05  HR 06  X<>Y 07  P-R 08  RCL 01  09  HR 10  X<>Y 11  P-R 12  STO 07  13  RDN 14  STO 08           15  X<>Y 16  STO 09  17  RCL 06  18  RCL 05 19  XEQ "GEOC" 20  HR 21  X<>Y 22  P-R 23  RCL 04  24  HR 25  X<>Y 26  P-R 27  STO 10 28  RCL 07 29  - 30  X^2 31  X<>Y 32  STO 11  33  RCL 08          34  - 35  X^2 36  + 37  X<>Y 38  STO 12  39  RCL 09          40  - 41  X^2 42  + 43  SQRT 44  END

( 63 bytes / SIZE 013 )

 STACK INPUT OUTPUT X / D (km)

Example:     Find the euclidean distance between the Mount Palomar Observatory:   L1 = -116°51'50"4  ,  b1 = 33°21'22"4  ,  h1 = 1706 m = 1.706 km
and the "Observatoire du Pic du Midi"   L2 = 0°08'32"4  ,  b2 = 42°56'12"0  ,  h2 = 2861 m = 2.861 km

-116.51504    STO 01        0.08324      STO 04
33.21224    STO 02       42.56120     STO 05
1.706        STO 03         2.861         STO 06

XEQ "EDST"  >>>>   D = 8585.5760 km   ( execution time = 11 seconds )

-And we have the geocentric coordinates:

R07 = x1 = -2410.4237                  R10 = x2 = 4678.8290        ( all in kilometers )
R08 = y1 = -4758.6127                  R11 = y2 =     11.6231
R09 = z1 =   3487.9636                  R12 = z2 = 4324.3023

-In comparison, the geodesic distance is  9410.6525 km  and the "loxodromic distance" is  10284.7558 km

Note:

-Delete lines 20-21-22 and 05-06-07 if you have deleted lines 30-31-32 in the "GEOC" listing.
-This reduces the execution time to 8.6 seconds.

References:

[1]  Jean Meeus - "Astronomical Algorithms" - Willmann-Bell -  ISBN 0-943396-35-2
[2]  Jacqueline Lelong-Ferrand - "Geometrie Differentielle" - Masson - ( in French )
[3]  Vincenty's formula
[4]  Paul D. Thomas - "Spheroidal Geodesics, Reference Systems & Local Geometry" - US Naval Oceanographic Office.
[5]  Paul D. Thomas - "Mathematical Models for Navigation Systems" - TR-182 - US Naval Oceanographic Office.
[6]  Emanuel M. Sodano - "General Non-Iterative Solution of the Inverse and Direct Geodetic Problems"
[7]  Richard H. Rapp - "Geometric Geodesy, Volume 11"