hp41programs

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  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 dµ

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


       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<Y?
208 STO 04
209 VIEW 04
210 XEQ 01
211 ENTER
212 ENTER
213 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<Y?
226 GTO 02
227 RCL 01
228 6
229 *
230 SIN
231 RCL 02
232 6
233 *
234 SIN
235 -
236 R-D
237 RCL 10
238 RCL 13
239 *
240 ST* Y
241 X^2
242 *
243 3072
244 /
245 RCL 14          
246 ENTER
247 ST* Y
248 ST- Y
249 ST* Y
250 -		 251 RCL 10
252 ST* Y
253 +
254 .75
255 *
256 RCL 13
257 X^2
258 +
259 RCL 10
260 16
261 /
262 X^2
263 *
264 RCL 16
265 *
266 +
267 RCL 14
268 15
269 *
270 3
271 +
272 RCL 14
273 *
274 3
275 -
276 RCL 14
277 *
278 15
279 -
280 RCL 10
281 *
282 32
283 /
284 RCL 14
285 X^2
286 +
287 RCL 10
288 ST* Y
289 -
290 4
291 /
292 RCL 13
293 -
294 RCL 10          
295 *
296 8
297 /
298 RCL 11
299 *
300 +		 301 RCL 14
302 5
303 *
304 3
305 +
306 RCL 14
307 *
308 3
309 +
310 RCL 14
311 *
312 5
313 +
314 RCL 10
315 *
316 4
317 /
318 RCL 14
319 3
320 *
321 2
322 +
323 RCL 14
324 *
325 3
326 +
327 +
328 RCL 10
329 *
330 16
331 /
332 RCL 14
333 +
334 RCL 10
335 ST* Y
336 +
337 4
338 /
339 RCL 12
340 ST* Y
341 -
342 -
343 RCL 04
344 RCL 05          
345 *
346 -
347 D-R
348 6378.137
349 *
350 CLD
351 END
 

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