hp41programs

Geodesic

Terrestrial Geodesic Distance for the HP-41


Overview
 

   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

     c) Forward Geodesic Problem

   2°) Vincenty's Formulae

     a) Geodesic Distance and Azimuths
     b) Forward Geodesic Problem

   3°) Numerical Integration
   4°) A Series Expansion
   5°) Euclidean Distance

    a) Geocentric Coordinates
    b) 3D-Distance
 

-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.
 

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  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       L1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       L2 ( ° ' " )             /
           X       b2 ( ° ' " )        D ( km )

 
Example:

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

-The error is only 15 mm!
-The accuracy is very good ( less than 1 meter ), provided D is not too large:
-Otherwise, the errors may be greater: for example, ( 0° , 0° ) ( 179° , 1° )  give 19860.3438 km whereas the correct result is 19860.5092 km
 

       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° !!!
 

     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  X<Y?
135  GTO 01 
136  2
137  RCL 00 
138  ST- Y
139  *
140  RCL 06 
141  COS
142  RCL 09
143  /
144  X^2
145  *
146  12
147  RCL Y
148  5
149  *
150  -
151  *
152  64
153  -
154  *
155  256
156  ST- Y
157  /
158  STO 08       
159  CLX
160  37
161  *
162  64
163  -
164  *
165  512
166  /
167  4
168  1/X
169  +
170  *
171  RCL 07 
172  X^2
173  ST+ X
174  1
175  -
176  RCL 05 
177  COS
178  4
179  /
180  *
181  *
182  RCL 07
183  +
184  *
185  RCL 05
186  SIN
187  *
188  RCL 05
189  D-R
190  -
191  RCL 08
192  *
193  ST* 09
194  RCL 04
195  SIN
196  STO 03       
197  RCL 01
198  COS
199  STO 05
200  CHS
201  ST* Y
202  RCL 02
203  SIN
204  ST* 05
205  *
206  RCL 04
207  COS
208  STO 08 
209  *
210  RCL 01 
211  SIN
212  ST* 08
213  RCL 02 
214  COS
215  ST* 03
216  ST* 08
217  *
218  +
219  R-P
220  X<>Y
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  X<Y?
119  GTO 01
120  RCL 04 
121  SIN
122  STO 04
123  RCL 10
124  COS
125  STO 01
126  *
127  RCL 13 
128  ST* 01
129  RCL 10
130  SIN
131  STO 05
132  *
133  RCL 02
134  COS
135  STO 13
136  *
137  +
138  RCL 04
139  RCL 05
140  *
141  RCL 01
142  RCL 13
143  *
144  -
145  X^2 
146  RCL 06       
147  X^2
148  +
149  SQRT
150  RCL 09
151  *
152  R-P
153  X<>Y
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<Y?
169  STO 15
170  VIEW 15
171  XEQ 01
172  ENTER^
173  ENTER^
174  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<Y?
187  GTO 03
188  VIEW 15
189  XEQ 01
190  STO 16 
191  4
192  STO 03
193  CF 10
194  XEQ 02
195  RCL 15
196  RCL 16
197  *
198  -
199  LBL 04
200  D-R
201  6378.137
202  *
203  CLA
204  CLD
205  END

 
    ( 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 2mn45s ) , 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 2mn10s ) , 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°) A Series Expansion
 

-"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 = ( µ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<Y?
167  STO 04
168  VIEW 04
169  XEQ 01
170  ENTER^
171  ENTER^
172  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<Y?
185  GTO 03
186  RCL 11
187  ST+ X
188  RCL 02
189  4
190  *
191  SIN
192  RCL 01 
193  4
194  *
195  SIN
196  -
197  4
198  /
199  R-D
200  -
201  RCL 12
202  3
203  *
204  -
205  RCL 13 
206  *
207  8
208  /
209  RCL 11
210  2
211  /
212  -
213  RCL 12        
214  +
215  RCL 13 
216  *
217  RCL 12
218  -
219  RCL 10 
220  X^2
221  *
222  8
223  /
224  RCL 14
225  4
226  /
227  -
228  RCL 12
229  +
230  RCL 04
231  RCL 05 
232  *
233  -
234  D-R
235  6378.137
236  *
237  CLD
238  END

 
  ( 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.
 

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.