hp41programs

Geodesic2

Terrestrial Geodesic Distance(2) for the HP-41


Overview
 

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

   a)  Ellipsoid of Revolution
   b)  Triaxial Ellipsoid ( with parametric coordinates )
   c)  
Triaxial Ellipsoid ( without parametric coordinates )

 3°)  Andoyer's Formula of order 2 + a few terms of order 3

   a)  Ellipsoid of Revolution

   b)  Triaxial Ellipsoid ( without parametric coordinates )
   c)  
Triaxial Ellipsoid ( with parametric coordinates )

 4°)  With Vincenty's Formula
 5°)  Minimizing the Sum of 2 Distances
 6°)  Other Formulae for the Great Elliptic Distance

   a)  Formula of order 1 + Andoyer's formula of order 2
   b)  Formula of order 1 + Andoyer's formula of order 3
   c)  F
ormula of order 2 + Andoyer's formula of order 3

 7°)  Geodesic + Great Elliptic Distance ( Ellipsoid of Revolution )



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

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

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

 

0°)  Method
 

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

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

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

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

    D ~  Δ + ( d - Δ' ) 

-In other words, the part of the formula giving the great elliptic arc distance on the ellipsoid of revolution in the WGS84 geodesic distance
  is replaced by the great elliptic arc distance on the triaxial ellipsoid.

 

1°)  With Andoyer's Formula of order 2
 

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

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

   OP = ( Sin W ) / ( A2 + B2 + C2 )1/2    for a point  P  on the ellipsoid  and   with  µ = ( OM , OP )

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

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

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

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

    P  =  -2 [ ( x Sin W ) ( x Cos W - x' ) / a2 + ( y Sin W ) ( y Cos W - y' ) / b2 + ( z Sin W ) ( z Cos W - z' ) / c2 ]
    Q =  [ x2 Sin2 W - ( x Cos W - x' )2 ] / a2 + [ y2 Sin2 W - ( y Cos W - y' )2 ] / b2 + [ z2 Sin2 W - ( z Cos W - z' )2 ] / c2

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

    A2 + B2 + C2 = T Sin2 ( W - µ ) + S Sin2 µ + 2 U Sin µ Sin ( W - µ )
   where  S = x'2 / a2 + y'2 / b2 + z'2 / c2  ,   T =  x2 / a2 + y2 / b2 + z2 / c2  ,   U =  xx' / a2 + yy' / b2 + zz' / c2
    and   P = 2 U - 2 T Cos W  ,  Q = T Sin W + ( 2 U Cos W - S - T Cos2 W ) / Sin W 

>>> But we can save a few bytes if we measure the angle µ from the bisector of the angle ( OM , ON )

 The formulae become:       A2 + B2 + C2 = T Sin2 ( W/2 - µ ) + S Sin2 ( W/2 + µ ) + 2 U Sin ( W/2 - µ ) Sin ( W/2 + µ )

  and    Tan 2µ = [ Sin W ( T - S ) ] / [ Cos W ( T + S ) - 2.U ]
 

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

-We can transform the usual formulas employing parametric latitudes into formulas employing geocentric latitudes ( cf page 91 of reference [4] ). It yields:

  D/a = W + (f/2) (-W+(Cos2F)(SinW))  +(f2 /4) [W/4+(SinW)(CosW)/4-(Cos2 (2F))(SinW)(CosW)/2]

     a = major axis of the great ellipse , f = eccentricity
   2F = b1+b2 where b = geocentric "latitudes"

 ( formulas of order 2 give a good accuracy to compute great elliptic distance for the Earth )

-Finally, the geodesic distance d on the WGS84 ellipsoid is calculated by 2nd order Andoyer's formula

 and the geodesic distance D on the triaxial ellipsoid is approximated by   D ~  Δ + ( d - Δ' )
 

Andoyer's Formula of order 2:
 

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

              S = sin2 G  cos2 H  +  cos2 F  sin2 H  ;   µ = 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 on the biaxial ellipsoid 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 }

   and   f = flattening = 1/298.25722

And we also have:

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

  ( here,   µ = 2 Arcsin S1/2 )

-But with a formula of order 2,  d - Δ' is computed by the same formula if  µ , A , B  are computed with geodetic latitudes !


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


 
 

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

 
      ( 358 bytes / SIZE 016 )
 
 

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

   Where  D is the geodesic distance on the triaxial ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                                                ---Execution time = 23s---
                                            LASTX  d  = 12138.65875 km = great elliptic distance = R10

-The exact result ( computed with Jacobi's method ) is D = 12138.657551942....
 

Notes:

-Longitudes are geodetic and positive East.

-When line 99 is executed ( STO 07 ),  the semi axes of the great ellipse =  R15/R06 & R15/R07.
-With the example above,  6378.151662 km & 6368.318461 km

-If the points are not nearly antipodal, the precision is of the order of one meter.

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

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

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

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

2°)  With Andoyer's Formula of order 3
 

     a) Ellipsoid of Revolution
 

-Sodano's formula is a 3rd order formula in the flattening to calculate the geodesic distance on an ellipsoid of revolution ( WGS84 ) - cf references [2] & [3]

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

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

          +  f2 [ p ((m-1)/2) µ2 Csc µ + m ((1-m)/2) µ2 Cot µ + ((m2-8.p2)/16) Sin µ Cos µ + (m2/16) µ - (m2/8) Sin µ Cos3 µ + (mp/2) Sin µ Cos2 µ ]

          +  f3 [ µ2 Csc µ (p/4) (-2+5m-3m2) + µ3 Csc µ Cot µ (p/2) (1-4m+3m2) - µ2 Cot µ (m/4) (-2+5m-3m2) + µ3 (m/6) (-1+2m-m2)
                + µ3 Cot2 µ (m/2) (-1+3m-2m2) + Sinµ Cos µ (-8 p2 + m2 - 8 p2m - m3/2 )/16 + µ  (8 p2 + m2 - 8 p2m - m3/2 )/16
                + µ3 Csc2 µ (p2/2) (1-m) + Sin µ Cos3 µ (m2/16) (-2+m) + µ Cos2 µ (m2/4) (1-m) + Sin µ Cos2 µ (pm/2) (1+m) + p2m Sin3 µ Cos µ
                + µ Cos µ (3pm/4) (-1+m) + (pm2/2) Sin5 µ + Sin µ (p/2) (p2-m2) + (m3/12) Sin3 µ Cos3 µ - (2p3/3) Sin3 µ ]  }

  where    µ = 2 Arc Sin [ sin2 ((B'-B)/2)  cos2 ((L'-L)/2)  +  cos2 ((B+B')/2)  sin2 ((L'-L)/2) ]1/2     L , L' = longitudes  &  B , B' = parametric latitudes

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

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

 
-But Sodano's formula may be re-written as an Andoyer's formula. 

 ( In reference [1] - page 87 - Paul D. Thomas gives the key to perform this transformation )

  With      L , L' = longitudes  &  p , p' = parametric latitudes  ,  f = earth flattening

               Tan p = (1-f) Tan b  where  b = geodetic latitude

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

              S = sin2 G  cos2 H  +  cos2 F  sin2 H ,  µ = 2 Arc Sin S1/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 )

-After "some" calculations, it yields:

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

                + (f2/4) [ A µ2 Cot µ + B µ2 / Sin µ + A2 ( µ + Sin µ Cos µ - 4 µ2 Cot µ ) / 4 - A.B µ2 / Sin µ - ( B2/2 ) Sin µ Cos µ ]

                + (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 + (1/2) Sin µ Cos µ + µ /2 )
                          + A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 + (1/4) Sin µ Cos µ - µ / 4 )
                          + B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ + µ Cos µ )
                        + A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ - µ Cos µ + (1/2) Sin µ Cos2 µ )
                        + B2 ( - Sin µ Cos µ + µ3 / Sin2 µ  + µ ) + A B2 ( (1/2) Sin µ Cos µ - µ3 / Sin2 µ - µ )
                        + B3 ( (2/3) Sin3 µ - (1/2) Sin µ ) ] }
 

Data Registers:   R00 = µ , R01 = A , R02 = B , R04 thru R09: temp
Flags: /
Subroutines: /
 
 

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

 
     ( 268 bytes / SIZE 010 )
 
 

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

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

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

Notes:

-It is certainly possible to express the formula without calculating the parametric latitudes, but I've not found how.
-Without parametric latitudes, Andoyer's formula of order 2 gives a shorter & faster program,
  but it's not at all sure with Andoyer's formula of order 3...

-We can also use an M-code routine to calculate the parametric latitudes:

( "PRM" doesn't check for alpha data )

0BD  "M"
012   "R"
010   "P"
2A0   SETDEC
0F8   C=X
22D   C=
060   1/C
2BE   C=-C
02E   C
0FA  ->
0AE  AB
001   C=
060   AB+1
128   L=C
0B8  C=Y
070   C
209   =
048  Tan(C)
138  C=L
13D  C=
060   AB*C
10E   A=C
0B8   C=Y
11E   A=C MS
0AE  A<>C
070   C
2A9  =
040  Atan(C)
0E8  X=C
078  C=Z
070  C
209  =
048  Tan(C)
138   C=L
13D  C=
060   AB*C
10E   A=C
078   C=Z
11E   A=C MS
0AE  A<>C
070   C
2A9  =
040  Atan(C)
068  Z=C
0A8 Y=C
3E0  RTN


           STACK           INPUTS         OUTPUTS
               Z                b              bparam
               Y                b'             bparam
               X               1/f             b'param
 

 

     b) Triaxial Ellipsoid  ( With Parametric Coordinates )
 

-The following "TGD" employs Andoyer's formula of order 2 for the great elliptic distance and Andoyer's formula of order 3 to estimate  d - Δ'   


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

                      + (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 + Sin µ Cos µ - µ  )
                          + A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 - Sin µ Cos µ + µ )
                          + B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ + µ Cos µ - 2 Sin µ )
                        + A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ - µ Cos µ + 2 Sin µ )
                        + B2 ( - 2 Sin µ Cos µ + µ3 / Sin2 µ  + µ ) + A B2 ( 2 Sin µ Cos µ - µ3 / Sin2 µ - µ )

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


 
 

 01 LBL "TGD"
 02 DEG
 03 HR
 04 STO 14
 05 RDN
 06 HR
 07 14.92911
 08 STO 11
 09 +
 10 STO 13
 11 RDN
 12 HR
 13 STO 12
 14 X<>Y
 15 HR
 16 ST+ 11
 17 6378.172
 18 STO 01
 19 .07
 20 -
 21 STO 02
 22 21.349686
 23 -
 24 STO 03
 25 RCL 11
 26 RCL 12
 27 XEQ 01
 28 STO 05
 29 RCL 09
 30 STO 06
 31 RCL 10
 32 STO 07
 33 RCL 13
 34 ST- 11
 35 RCL 14
 36 XEQ 01
 37 XEQ 00
 38 STO 04
 39 RCL 01
 40 ST/ 05
 41 ST/ 08
 42 RCL 02
 43 ST+ 01
 44 ST/ 06
 45 ST/ 09
 46 RCL 03          
 47 ST/ 07
 48 ST/ 10
 49 RCL 08
 50 X^2
 51 RCL 09
 52 X^2
 53 RCL 10
 54 X^2
 55 +
 56 +
 57 STO 13
 58 RCL 05
 59 X^2
 60 RCL 06
 61 X^2
 62 RCL 07
 63 X^2
 64 +
 65 +
 66 STO 00
 67 -
 68 RCL 04
 69 ACOS
 70 STO 02
 71 SIN
 72 STO 15
 73 *
 74 RCL 08
 75 XEQ 00
 76 ST+ X
 77 STO 09
 78 RCL 00
 79 RCL 13
 80 +
 81 RCL 04
 82 *
 83 -
 84 X#0?
 85 /
 86 ATAN
 87 STO 08
 88 2
 89 ST/ 01
 90 ST/ 11
 91 /
 92 STO 05
 93 XEQ 02
 94 STO 06
 95 RCL 05
 96 90
 97 +
 98 XEQ 02
 99 STO 07
100 RCL 02          
101 D-R
102 RCL 08
103 COS
104 STO 05
105 X^2
106 ST+ X
107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X
121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 GTO 03
132 LBL 00
133 RCL 05
134 *
135 RCL 09
136 RCL 06
137 *
138 +
139 RCL 10
140 RCL 07
141 *
142 +
143 RTN
144 LBL 01
145 1
146 P-R
147 X<>Y
148 RCL 03
149 X^2
150 *
151 STO 10          
152 RDN
153 P-R
154 RCL 01
155 X^2
156 *
157 STO Z
158 X^2
159 X<>Y
160 RCL 02
161 X^2
162 *
163 STO 09
164 X^2
165 RCL 10
166 X^2
167 +
168 +
169 SQRT
170 ST/ 09
171 ST/ 10
172 /
173 STO 08
174 RTN
175 LBL 02
176 RCL 02
177 2
178 /
179 +
180 ENTER
181 SIN
182 ENTER
183 X^2
184 RCL 13
185 *
186 RCL 02
187 R^
188 -
189 SIN
190 ST* Z
191 X^2
192 RCL 00
193 *
194 +
195 X<>Y
196 RCL 09
197 *
198 +
199 SQRT
200 RTN
201 LBL 03
202 RCL 14
203 1
204 STO 08          
205 STO 09
206 P-R
207 RCL 01
208 RCL 03
209 -
210 RCL 01
211 /
212 STO 04
213 SIGN
214 LASTX
215 -
216 STO 05
217 /
218 R-P
219 X<> 12
220 RCL 09
221 P-R
222 RCL 05
223 /
224 R-P
225 RDN
226 +
227 2
228 ST/ 04
229 /
230 ST- Y
231 RCL 09
232 P-R
233 X^2
234 STO T
235 RDN
236 X^2
237 STO 14
238 SIGN
239 P-R
240 X^2
241 ST* 14
242 CLX
243 +
244 X^2
245 ST* T
246 -
247 RCL 11
248 COS
249 X^2
250 *
251 -
252 /
253 STO 02
254 RCL 14
255 RCL 09
256 LASTX
257 STO 14          
258 -
259 /
260 ST- 02
261 +
262 ST- 09
263 X<> 14
264 SQRT
265 ASIN
266 ST+ X
267 D-R
268 STO 00
269 LASTX
270 RCL 08
271 P-R
272 STO 03
273 STO 06
274 STO 13
275 RDN
276 STO 05
277 /
278 STO 07
279 ST/ Z
280 X^2
281 ST* 03
282 RCL 08
283 +
284 R^
285 ST- 08
286 ST+ X
287 -
288 RCL 02
289 *
290 2
291 RCL 07
292 ST- 03
293 ST/ Y
294 -
295 RCL 06
296 -
297 RCL 03
298 ST+ 03
299 ST* 13
300 4
301 *
302 -
303 RCL 09
304 ST* Z
305 *
306 RCL 03
307 -
308 RCL 14          
309 *
310 -
311 RCL 03
312 -
313 RCL 02
314 *
315 RCL 13
316 RCL 00
317 ST* 01
318 ST+ X
319 X^2
320 3
321 ST* Z
322 /
323 STO 11
324 +
325 STO 12
326 RCL 08
327 -
328 RCL 14
329 *
330 RCL 12
331 ST+ X
332 -
333 RCL 13
334 +
335 RCL 08
336 +
337 RCL 14
338 *
339 RCL 13
340 ST+ X
341 +
342 RCL 11
343 +
344 RCL 14
345 *
346 -
347 RCL 04
348 *
349 RCL 06
350 RCL 07
351 *
352 RCL 14
353 ST* Y
354 -
355 RCL 07
356 ENTER
357 1/X
358 -
359 RCL 02
360 *
361 +
362 RCL 09          
363 *
364 +
365 RCL 04
366 X^2
367 *
368 RCL 01
369 *
370 RCL 15
371 +
372 END

 
      ( 466 bytes / SIZE 016 )
    
 

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

   Where  D is the geodesic distance on the triaxial ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                      ---Execution time = 28s--- 

                                     LASTX = R10 =  12138.65875 km = great elliptic arc distance

Notes:

-The exact result ( computed with Jacobi's method ) is  12138.657551942....

-When line 99 is executed ( STO 07 ),  the semi axes of the great ellipse =  R15/R06 & R15/R07.
-With the example above,  6378.151662 km & 6368.318461 km

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


     c) Triaxial Ellipsoid  ( Without Parametric Coordinates )


-The formula employed in the previous paragraph may be expressed without transforming geodetic coordinates into parametric coordinates.
-After some calculations, it yields:


   (d - Δ') / a  =  (f2/4) µ [ (A2-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]                 where  µ , A , B are computed with geodetic latitudes

              + (f3/8) [ A ( 4 µ - 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 )
                         + A2 ( 9 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 - Sin µ Cos µ - 13 µ  )
                         + A3 ( -7 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 + Sin µ Cos µ + 9 µ )
                         + B ( -2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ + 4 Sin µ )
                         + A.B ( 9 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ - 3 µ Cos µ - 12 Sin µ )
                        + A2 B ( -7 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ + 3 µ Cos µ + 8 Sin µ )
                        + B2 ( 2 Sin µ Cos µ + µ3 / Sin2 µ  - 3 µ )
                        + A B2 ( -2 Sin µ Cos µ - µ3 / Sin2 µ + 3 µ )

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


 
 

 01 LBL "TGD"
 02 DEG
 03 HR
 04 STO 14
 05 RDN
 06 HR
 07 14.92911
 08 STO 11
 09 +
 10 STO 13
 11 RDN
 12 HR
 13 STO 12
 14 X<>Y
 15 HR
 16 ST+ 11
 17 6378.172
 18 STO 01
 19 .07
 20 -
 21 STO 02
 22 21.349686
 23 -
 24 STO 03
 25 RCL 11
 26 RCL 12
 27 XEQ 01
 28 STO 05
 29 RCL 09
 30 STO 06
 31 RCL 10
 32 STO 07
 33 RCL 13
 34 ST- 11
 35 RCL 14
 36 XEQ 01
 37 XEQ 00
 38 STO 04
 39 RCL 01
 40 ST/ 05
 41 ST/ 08
 42 RCL 02
 43 ST+ 01
 44 ST/ 06
 45 ST/ 09
 46 RCL 03          
 47 ST/ 07
 48 ST/ 10
 49 RCL 08
 50 X^2
 51 RCL 09
 52 X^2
 53 RCL 10
 54 X^2
 55 +
 56 +
 57 STO 13
 58 RCL 05
 59 X^2
 60 RCL 06
 61 X^2
 62 RCL 07
 63 X^2
 64 +
 65 +
 66 STO 00
 67 -
 68 RCL 04
 69 ACOS
 70 STO 02
 71 SIN
 72 STO 15
 73 *
 74 RCL 08
 75 XEQ 00
 76 ST+ X
 77 STO 09
 78 RCL 00
 79 RCL 13
 80 +
 81 RCL 04
 82 *
 83 -
 84 X#0?
 85 /
 86 ATAN
 87 STO 08
 88 2
 89 ST/ 01
 90 ST/ 11
 91 /
 92 STO 05
 93 XEQ 02
 94 STO 06
 95 RCL 05
 96 90
 97 +
 98 XEQ 02
 99 STO 07
100 RCL 02          
101 D-R
102 RCL 08
103 COS
104 STO 05
105 X^2
106 ST+ X
107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X
121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 RCL 14
132 ENTER
133 GTO 03
134 LBL 00
135 RCL 05
136 *
137 RCL 09
138 RCL 06
139 *
140 +
141 RCL 10
142 RCL 07
143 *
144 +
145 RTN
146 LBL 01
147 1
148 P-R
149 X<>Y
150 RCL 03          
151 X^2
152 *
153 STO 10
154 RDN
155 P-R
156 RCL 01
157 X^2
158 *
159 STO Z
160 X^2
161 X<>Y
162 RCL 02
163 X^2
164 *
165 STO 09
166 X^2
167 RCL 10
168 X^2
169 +
170 +
171 SQRT
172 ST/ 09
173 ST/ 10
174 /
175 STO 08
176 RTN
177 LBL 02
178 RCL 02
179 2
180 /
181 +
182 ENTER
183 SIN
184 ENTER
185 X^2
186 RCL 13
187 *
188 RCL 02
189 R^
190 -
191 SIN
192 ST* Z
193 X^2
194 RCL 00
195 *
196 +
197 X<>Y
198 RCL 09
199 *
200 +
201 SQRT
202 RTN
203 LBL 03
204 RCL 12          
205 +
206 2
207 /
208 ST- Y
209 1
210 STO 07
211 P-R
212 X^2
213 STO T
214 RDN
215 X^2
216 STO 08
217 SIGN
218 P-R
219 X^2
220 ST* 08
221 CLX
222 +
223 X^2
224 ST* T
225 -
226 RCL 11
227 COS
228 X^2
229 *
230 -
231 /
232 STO 02
233 RCL 08
234 RCL 07
235 LASTX
236 STO 08
237 -
238 /
239 ST- 02
240 +
241 ST- 07
242 X<> 08
243 SQRT
244 ASIN
245 ST+ X
246 D-R
247 STO 00
248 LASTX
249 1
250 P-R
251 STO 09
252 STO 12
253 STO 13          
254 RDN
255 /
256 STO 06
257 ST/ 12
258 X^2
259 ST* 09
260 3
261 -
262 RCL 12
263 ST+ X
264 +
265 RCL 02
266 *
267 3
268 RCL 06
269 ST* Y
270 RCL 09
271 +
272 ST+ X
273 STO 10
274 STO 11
275 ST+ X
276 +
277 ST- 09
278 RCL 13
279 3
280 *
281 -
282 8
283 RCL 06
284 /
285 -
286 RCL 07
287 ST* Z
288 *
289 RCL 10
290 +
291 4
292 STO 14
293 RCL 06
294 /
295 ST- 11
296 -
297 RCL 08
298 *
299 +
300 RCL 11
301 -
302 RCL 02
303 *
304 RCL 09
305 RCL 13          
306 ST* 10
307 *
308 RCL 12
309 +
310 9
311 +
312 RCL 00
313 ST+ X
314 X^2
315 3
316 /
317 ST- 14
318 -
319 RCL 07
320 *
321 RCL 10
322 -
323 RCL 14
324 ST- 10
325 +
326 RCL 08
327 *
328 RCL 10
329 +
330 RCL 08
331 *
332 -
333 RCL 01
334 RCL 03
335 -
336 RCL 01
337 ST+ X
338 /
339 STO 04
340 *
341 RCL 06
342 RCL 13
343 *
344 RCL 08
345 ST* Y
346 -
347 RCL 06
348 ENTER
349 1/X
350 -
351 RCL 02
352 *
353 +
354 RCL 07
355 *
356 +
357 RCL 04          
358 X^2
359 *
360 RCL 00
361 *
362 RCL 01
363 *
364 RCL 15
365 +
366 END

 
      ( 458 bytes / SIZE 016 )
    
 

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

   Where  D is the geodesic distance on the triaxial ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                      ---Execution time = 25s--- 

                                     LASTX = R15 =  12138.65875 km = great elliptic arc distance

Note:

-When line 99 is executed ( STO 07 ),  the semi axes of the great ellipse =  R15/R06 & R15/R07.
-With the example above,  6378.151662 km & 6368.318461 km


3°)  With Andoyer's Formula of order 2 + A Few Terms of Order 3
 

     a)  Ellipsoid of Revolution


-This version employs Andoyer's formula of order 2 ( where A , B , ... are computed with geodetic latitudes ) + a few terms of order 3


Formulae:

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

              S = sin2 G  cos2 H  +  cos2 F  sin2 H  ;   µ = 2 Arcsin S1/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 )

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

                + (f2/4) [ A ( 4.µ + µ2 Cot µ ) + B ( 6.Sin µ + µ2 / Sin µ ) - A2 ( 15.µ + 15 Sin µ Cos µ + 4 µ2 Cot µ ) / 4 - A.B ( 6.Sin µ + µ2 / Sin µ ) + ( B2/2 ) 15 Sin µ Cos µ ]

                + (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 )
                          + A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 )
                          + B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ )
                          + A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ  )
                          + B2 ( µ3 / Sin2 µ  ) + A B2 ( - µ3 / Sin2 µ )
                      

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

 
 

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

 
      ( 220 bytes / SIZE 008 )
 
 

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

   Where  D is the geodesic distance on the WGS84 ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D =  12138.68425 km                                          ---Execution time = 10s---

 
Notes:

-The flattening value ( x2: line 123 ) may be replaced by 596.514
-It will save 2 bytes without changing a lot the precision.

-The precision is the same as Andoyer's formula of order 2 for relatively long distances and almost the same as Andoyer's formula of order 3 for very long distances.
-For instance,  with L = 0°  b = 5°  L' = 175°  b' = 3°  this program gives  18968.13264 km  whereas Vincenty's formula gives 18968.13070 km
-So, an error of 1.94 m

-Here are a few maximum errors ( in meters ) for several  L' - L  values

-The second row = maximum error of Andoyer's formula of order 2 without any term of order 3
-Third row = maximum error with "TGD" above: Andoyer of order 2 + a few of order 3.
-Fourth row = maximum error with complete Andoyer's formula of order 3

-I found them with free42.


     | L'-L |
  90°
 150°
 155°
 160°
 165°
 170°
 171°
 172°
 173°
 174°
 175°
 176°
     And2
 0.16
 0.52
 0.85
 1.50
 3.03
 7.64
 9.76
 12.5
 16.9
 23.4
 34.4
 55.6
 And2afew3  0.16
 0.20
 0.26
 0.34
 0.51
 0.88
 1.00
 1.15
 1.33
 1.57
 1.94
 3.14
     And3
  ~ 0
0.008
 0.01
 0.03
 0.08
 0.27
 0.39
 0.52
 0.82
 1.36
 2.05
 4.38


-Thus, if  L'-L < 171° ,  the errors remain smaller than 1 meter  ( except for nearly antipodal points ) !

-Surprisingly, if L' - L > 174°,  the complete Andoyer's ( or Sodano's ) formula of order 3 may be less accurate than Andoyer's formula of order 2 + main terms of order 3 !

 
     b)  Triaxial Ellipsoid  ( Without Parametric Coordinates )

 

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

              S = sin2 G  cos2 H  +  cos2 F  sin2 H  ;   µ = 2 Arcsin S1/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 )

   d - Δ' = a (f2/4) µ [ (A2-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]  

                + a (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 )
                             + A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 )
                            + B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ )
                            + A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ  )
                            + B2 ( µ3 / Sin2 µ  ) + A B2 ( - µ3 / Sin2 µ )
                      

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


 
 01 LBL "TGD"
 02 DEG
 03 HR
 04 STO 14
 05 RDN
 06 HR
 07 14.929
 08 STO 11
 09 +
 10 STO 13
 11 RDN
 12 HR
 13 STO 12
 14 X<>Y
 15 HR
 16 ST+ 11
 17 6378.172
 18 STO 01
 19 .07
 20 -
 21 STO 02
 22 21.349686
 23 -
 24 STO 03
 25 RCL 11
 26 RCL 12
 27 XEQ 01
 28 STO 05
 29 RCL 09
 30 STO 06
 31 RCL 10
 32 STO 07
 33 RCL 13
 34 ST- 11
 35 RCL 14
 36 XEQ 01
 37 XEQ 00
 38 STO 04
 39 RCL 01
 40 ST/ 05
 41 ST/ 08
 42 RCL 02          
 43 ST+ 01
 44 ST/ 06
 45 ST/ 09
 46 RCL 03
 47 ST/ 07
 48 ST/ 10
 49 RCL 08
 50 X^2
 51 RCL 09
 52 X^2
 53 RCL 10
 54 X^2
 55 +
 56 +
 57 STO 13
 58 RCL 05
 59 X^2
 60 RCL 06
 61 X^2
 62 RCL 07
 63 X^2
 64 +
 65 +
 66 STO 00
 67 -
 68 RCL 04
 69 ACOS
 70 STO 02
 71 SIN
 72 STO 15
 73 *
 74 RCL 08
 75 XEQ 00
 76 ST+ X
 77 STO 09
 78 RCL 00
 79 RCL 13
 80 +
 81 RCL 04
 82 *
 83 -
 84 X#0?
 85 /
 86 ATAN
 87 STO 08
 88 2
 89 ST/ 01
 90 ST/ 11
 91 /
 92 STO 05          
 93 XEQ 02
 94 STO 06
 95 RCL 05
 96 90
 97 +
 98 XEQ 02
 99 STO 07
100 RCL 02
101 D-R
102 RCL 08
103 COS
104 STO 05
105 X^2
106 ST+ X
107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X
121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 RCL 14
132 ENTER
133 GTO 03
134 LBL 00
135 RCL 05
136 *
137 RCL 09
138 RCL 06
139 *
140 +
141 RCL 10          
142 RCL 07
143 *
144 +
145 RTN
146 LBL 01
147 1
148 P-R
149 X<>Y
150 RCL 03
151 X^2
152 *
153 STO 10
154 RDN
155 P-R
156 RCL 01
157 X^2
158 *
159 STO Z
160 X^2
161 X<>Y
162 RCL 02
163 X^2
164 *
165 STO 09
166 X^2
167 RCL 10
168 X^2
169 +
170 +
171 SQRT
172 ST/ 09
173 ST/ 10
174 /
175 STO 08
176 RTN
177 LBL 02
178 RCL 02
179 2
180 /
181 +
182 ENTER
183 SIN
184 ENTER
185 X^2
186 RCL 13
187 *
188 RCL 02          
189 R^
190 -
191 SIN
192 ST* Z
193 X^2
194 RCL 00
195 *
196 +
197 X<>Y
198 RCL 09
199 *
200 +
201 SQRT
202 RTN
203 LBL 03
204 RCL 12
205 +
206 2
207 /
208 ST- Y
209 1
210 STO 07
211 P-R
212 X^2
213 STO T
214 RDN
215 X^2
216 STO 08
217 SIGN
218 P-R
219 X^2
220 ST* 08
221 CLX
222 +
223 X^2
224 ST* T
225 -
226 RCL 11
227 COS
228 X^2
229 *
230 -
231 /
232 STO 02
233 RCL 08
234 RCL 07          
235 LASTX
236 STO 08
237 -
238 /
239 ST- 02
240 +
241 ST- 07
242 X<> 08
243 SQRT
244 ASIN
245 ST+ X
246 D-R
247 STO 00
248 LASTX
249 1
250 P-R
251 STO 04
252 STO 05
253 STO 09
254 RDN
255 /
256 ST* 05
257 STO 06
258 X^2
259 ST* 09
260 RCL 02
261 *
262 RCL 06
263 ST- 09
264 RCL 09
265 ST* 04
266 ST+ 09
267 4
268 *
269 +
270 RCL 07
271 ST* Z
272 *
273 RCL 09
274 +
275 RCL 08
276 *
277 +
278 RCL 09
279 -
280 RCL 02
281 *
282 RCL 04          
283 RCL 00
284 ST+ X
285 X^2
286 3
287 ST* Z
288 /
289 STO Z
290 +
291 STO 09
292 RCL 08
293 *
294 RCL 09
295 ST+ X
296 -
297 RCL 04
298 ST+ Z
299 ST+ Z
300 +
301 RCL 08
302 ST* 05
303 *
304 +
305 RCL 08
306 ST- 05
307 *
308 -
309 RCL 01
310 RCL 03
311 -
312 RCL 01
313 ST+ X
314 /
315 STO 04
316 *
317 RCL 05
318 RCL 06
319 ENTER
320 1/X
321 -
322 RCL 02
323 *
324 +
325 RCL 07
326 *
327 +
328 RCL 04
329 X^2
330 *
331 RCL 00          
332 *
333 RCL 01
334 *
335 RCL 15
336 +
337 END

 
      ( 427 bytes / SIZE 016 )
 
 

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

   Where  D is the geodesic distance on the triaxial ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                                                ---Execution time = 25s---
                                            LASTX  d  = 12138.65875 km = great elliptic distance = R10

-The exact result ( computed with Jacobi's method ) is D = 12138.657551942....


Notes:  

-When line 99 is executed ( STO 07 ),  the semi axes of the great ellipse =  R15/R06 & R15/R07.
-With the example above,  6378.151662 km & 6368.318461 km

-The precision is not so good as the programs listed in paragraphs 2b) & 2c)
-For example:

-If   L = -40° & L' = 120° ,  the root mean square error = 15 cm  and the maximum | error | = 55 cm
-If   L = -95° &  L' = 65° ,  rms = 31 cm and max | error |  = 87 cm


     c)  Triaxial Ellipsoid  ( With Parametric Coordinates )


-The method employed in paragraph above is not really rigorous !
-A better accuracy will be obtained if we use the same formulas but after transforming the geodetic coordinates into parametric coordinates to calculate  µ , A , B ... etc ...


   d - Δ' = a (f2/4) µ [ (A2-A) ( 1 - µ Cot µ ) - B.(A-1) ( µ / Sin µ - (Sin µ) / µ ) ]  

                + a (f3/8) [ A ( 2 µ2 Cot µ - 2 µ3 Cot2 µ - 8 µ3 / 6 ) + A2 ( -5 µ2 Cot µ + 5 µ3 Cot2 µ + 8 µ3 / 3 )
                             + A3 ( 3 µ2 Cot µ - 3 µ3 Cot2 µ - 8 µ3 / 6 )
                            + B ( 2 µ2 / Sin µ - 2 µ3 Cos µ / Sin2 µ ) + A.B ( -5 µ2 / Sin µ + 6 µ3 Cos µ / Sin2 µ )
                            + A2 B ( 3 µ2 / Sin µ - 4 µ3 Cos µ / Sin2 µ  )
                            + B2 ( µ3 / Sin2 µ  ) + A B2 ( - µ3 / Sin2 µ )
                      

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


 
 01 LBL "TGD"
 02 DEG
 03 HR
 04 STO 14
 05 RDN
 06 HR
 07 14.92911
 08 STO 11
 09 +
 10 STO 13
 11 RDN
 12 HR
 13 STO 12
 14 X<>Y
 15 HR
 16 ST+ 11
 17 6378.172
 18 STO 01
 19 .07
 20 -
 21 STO 02
 22 21.349686
 23 -
 24 STO 03
 25 RCL 11
 26 RCL 12
 27 XEQ 01
 28 STO 05
 29 RCL 09
 30 STO 06
 31 RCL 10
 32 STO 07
 33 RCL 13
 34 ST- 11
 35 RCL 14
 36 XEQ 01
 37 XEQ 00
 38 STO 04
 39 RCL 01
 40 ST/ 05
 41 ST/ 08
 42 RCL 02
 43 ST+ 01
 44 ST/ 06
 45 ST/ 09
 46 RCL 03          
 47 ST/ 07
 48 ST/ 10
 49 RCL 08
 50 X^2
 51 RCL 09
 52 X^2
 53 RCL 10
 54 X^2
 55 +
 56 +
 57 STO 13
 58 RCL 05
 59 X^2
 60 RCL 06
 61 X^2
 62 RCL 07
 63 X^2
 64 +
 65 +
 66 STO 00
 67 -
 68 RCL 04
 69 ACOS
 70 STO 02
 71 SIN
 72 STO 15
 73 *
 74 RCL 08
 75 XEQ 00
 76 ST+ X
 77 STO 09
 78 RCL 00
 79 RCL 13
 80 +
 81 RCL 04
 82 *
 83 -
 84 X#0?
 85 /
 86 ATAN
 87 STO 08
 88 2
 89 ST/ 01
 90 ST/ 11
 91 /
 92 STO 05
 93 XEQ 02
 94 STO 06
 95 RCL 05          
 96 90
 97 +
 98 XEQ 02
 99 STO 07
100 RCL 02
101 D-R
102 RCL 08
103 COS
104 STO 05
105 X^2
106 ST+ X
107 RCL 04
108 RCL 15
109 ST* 05
110 *
111 ST* Y
112 -
113 -
114 4
115 /
116 RCL 06
117 ST/ 15
118 RCL 07
119 ST- Y
120 ST+ X
121 /
122 STO 07
123 *
124 RCL 05
125 -
126 +
127 RCL 07
128 *
129 +
130 ST* 15
131 GTO 03
132 LBL 00
133 RCL 05
134 *
135 RCL 09
136 RCL 06
137 *
138 +
139 RCL 10          
140 RCL 07
141 *
142 +
143 RTN
144 LBL 01
145 1
146 P-R
147 X<>Y
148 RCL 03
149 X^2
150 *
151 STO 10
152 RDN
153 P-R
154 RCL 01
155 X^2
156 *
157 STO Z
158 X^2
159 X<>Y
160 RCL 02
161 X^2
162 *
163 STO 09
164 X^2
165 RCL 10
166 X^2
167 +
168 +
169 SQRT
170 ST/ 09
171 ST/ 10
172 /
173 STO 08
174 RTN
175 LBL 02
176 RCL 02
177 2
178 /
179 +
180 ENTER
181 SIN
182 ENTER
183 X^2
184 RCL 13
185 *
186 RCL 02
187 R^
188 -
189 SIN
190 ST* Z
191 X^2
192 RCL 00          
193 *
194 +
195 X<>Y
196 RCL 09
197 *
198 +
199 SQRT
200 RTN
201 LBL 03
202 RCL 14
203 1
204 STO 07
205 P-R
206 RCL 01
207 RCL 03
208 -
209 RCL 01
210 /
211 STO 04
212 SIGN
213 LASTX
214 -
215 STO 05
216 /
217 R-P
218 X<> 12
219 RCL 07
220 P-R
221 RCL 05
222 /
223 R-P
224 RDN
225 +
226 2
227 ST/ 04
228 /
229 ST- Y
230 RCL 07
231 P-R
232 X^2
233 STO T
234 RDN
235 X^2
236 STO 08
237 SIGN
238 P-R
239 X^2
240 ST* 08
241 CLX
242 +
243 X^2
244 ST* T
245 -
246 RCL 11          
247 COS
248 X^2
249 *
250 -
251 /
252 STO 02
253 RCL 08
254 RCL 07
255 LASTX
256 STO 08
257 -
258 /
259 ST- 02
260 +
261 ST- 07
262 X<> 08
263 SQRT
264 ASIN
265 ST+ X
266 D-R
267 STO 00
268 LASTX
269 1
270 P-R
271 STO 05
272 STO 09
273 STO 13
274 RDN
275 /
276 ST* 05
277 STO 06
278 X^2
279 ST* 09
280 RCL 02
281 *
282 RCL 06
283 ST- 09
284 RCL 09
285 ST+ 09
286 ST* 13
287 4
288 *
289 +
290 RCL 07
291 ST* Z
292 *
293 RCL 09
294 +
295 RCL 08          
296 *
297 +
298 RCL 09
299 -
300 RCL 02
301 *
302 RCL 13
303 RCL 00
304 ST* 01
305 ST+ X
306 X^2
307 3
308 ST* Z
309 /
310 STO Z
311 +
312 STO 09
313 RCL 08
314 *
315 RCL 09
316 ST+ X
317 -
318 RCL 13
319 ST+ Z
320 ST+ Z
321 +
322 RCL 08
323 ST* 05
324 *
325 +
326 RCL 08
327 ST- 05
328 *
329 -
330 RCL 04
331 *
332 RCL 05
333 RCL 06
334 ENTER
335 1/X
336 -
337 RCL 02
338 *
339 +
340 RCL 07
341 *
342 +
343 RCL 04          
344 X^2
345 *
346 RCL 01
347 *
348 RCL 15
349 +
350 END

 
      ( 444 bytes / SIZE 016 )
 
 

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

   Where  D is the geodesic distance on the triaxial ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65755 km                                                                                ---Execution time = 27s---
                                         LASTX = d  = 12138.65875 km = great elliptic distance = R15

Notes:  

-When line 99 is executed ( STO 07 ),  the semi axes of the great ellipse =  R15/R06 & R15/R07.
-With the example above,  6378.151662 km & 6368.318461 km

-The precision is the same as the programs listed in paragraphs 2b) & 2c) !!!
-For example:

-If   L = -40° & L' = 120° ,  the root mean square error = 10 cm  and the maximum | error | = 48 cm
-If   L = -95° &  L' = 65° ,  rms = 19 cm and max | error |  = 55 cm

4°)  With Vincenty's Formula
 

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

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

-Line 47 is a three-byte GTO 14
 
 

 01 LBL "TGDV"
 02 DEG
 03 HR
 04 STO 21        
 05 RDN
 06 HR
 07 14.92911
 08 STO 00
 09 +
 10 STO 20
 11 RDN
 12 HR
 13 STO 19
 14 X<>Y
 15 HR
 16 RCL 00
 17 +
 18 STO 18          
 19 6378.172
 20 STO 01
 21 .07
 22 -
 23 STO 02
 24 6356.752314
 25 STO 03
 26 XEQ 01
 27 STO 17
 28 RCL 01
 29 RCL 02
 30 +
 31 2
 32 /
 33 STO 01
 34 STO 02
 35 XEQ 01
 36 ST- 17
 37 RCL 01
 38 STO 06
 39 RCL 03
 40 STO 07
 41 RCL 21
 42 RCL 19
 43 RCL 20
 44 RCL 18
 45 XEQ 13
 46 ENTER
 47 GTO 14
 48 LBL 01
 49 RCL 18
 50 RCL 19
 51 XEQ 02
 52 RCL 08
 53 STO 05
 54 RCL 09
 55 STO 06
 56 RCL 10
 57 STO 07
 58 RCL 20
 59 RCL 21
 60 XEQ 02
 61 RCL 05
 62 RCL 08
 63 *
 64 RCL 06
 65 RCL 09
 66 *
 67 RCL 07 
 68 RCL 10        
 69 *
 70 +
 71 +
 72 STO 00
 73 ACOS
 74 STO 04
 75 RCL 01
 76 ST/ 05
 77 ST/ 08
 78 RCL 02
 79 ST/ 06
 80 ST/ 09
 81 RCL 03          
 82 ST/ 07
 83 ST/ 10
 84 RCL 05
 85 X^2
 86 RCL 06
 87 X^2
 88 RCL 07
 89 X^2
 90 +
 91 +
 92 STO 11
 93 RCL 08
 94 X^2
 95 RCL 09
 96 X^2
 97 RCL 10
 98 X^2
 99 +
100 +
101 X<> 10
102 RCL 07
103 *
104 RCL 06
105 RCL 09
106 *
107 RCL 05
108 RCL 08
109 *
110 +
111 +
112 STO 09
113 RCL 00
114 *
115 ST+ X
116 RCL 10
117 -
118 RCL 11
119 RCL 00
120 X^2
121 *
122 -
123 RCL 11
124 RCL 04
125 SIN
126 STO 08
127 ST/ Z
128 *
129 +
130 RCL 09
131 RCL 11
132 RCL 00        
133 *
134 -
135 ST+ X
136 X<>Y
137 X#0?
138 /
139 ATAN
140 2
141 /
142 90
143 STO 00
144 X<>Y
145 X<0?
146 +
147 STO 05
148 XEQ 03
149 STO 06          
150 RCL 05
151 RCL 00
152 +
153 XEQ 03
154 STO 07
155 RCL 06
156 X>Y?
157 GTO 00
158 X<> 07
159 STO 06
160 RCL 00
161 ST+ 05
162 LBL 00
163 PI
164 R-D
165 RCL 05
166 X=Y?
167 CLX
168 STO 05
169 CHS
170 1
171 P-R
172 RCL 07
173 RCL 06
174 /
175 X^2
176 STO 09
177 *
178 R-P
179 CLX
180 RCL 04
181 RCL 05
182 -
183 1
184 P-R
185 RCL 09
186 *
187 R-P
188 CLX
189 ENTER
190 LBL 13
191 -
192 360
193 MOD
194 PI
195 R-D
196 STO 16        
197 X>Y?
198 CLX
199 ST+ X
200 -
201 STO 11
202 STO 04
203 CLX
204 RCL 07
205 CHS
206 RCL 06
207 ST+ Y
208 /
209 STO 08
210 SIGN
211 LASTX
212 -
213 STO 15          
214 P-R
215 LASTX
216 /
217 R-P
218 RDN
219 STO 09
220 CLX
221 RCL 15
222 P-R
223 LASTX
224 /
225 R-P
226 X<>Y
227 STO 10
228 LBL 12
229 RCL 16
230 RCL 04
231 X>Y?
232 ACOS
233 RCL 10
234 COS
235 STO 05
236 STO 12
237 RCL 04
238 SIN
239 *
240 STO 14
241 RCL 09
242 COS
243 ST* 05
244 ST* 14
245 RCL 10
246 SIN
247 STO 13
248 *
249 RCL 09
250 SIN
251 ST* 13
252 RCL 12
253 *
254 RCL 04
255 COS
256 ST* 05
257 *
258 -
259 R-P
260 RCL 05        
261 RCL 13 
262 +
263 R-P
264 X<> 14
265 X<>Y
266 STO 05
267 SIN
268 X#0?
269 /
270 ASIN
271 STO 12
272 COS
273 X^2
274 RCL 05
275 COS
276 RCL 13          
277 ENTER
278 +
279 R^
280 X=0?
281 SIGN
282 /
283 -
284 STO 13
285 CLX
286 3
287 *
288 4
289 X<>Y
290 -
291 RCL 08
292 *
293 4
294 +
295 *
296 RCL 08
297 *
298 16
299 /
300 RCL 05
301 SIN
302 RCL 13
303 RCL Z
304 *
305 *
306 R-D
307 RCL 05
308 +
309 RCL 12
310 SIN
311 *
312 RCL 08
313 *
314 1
315 R^
316 -
317 *
318 RCL 11
319 +
320 ENTER
321 X<> 04
322 -
323 ABS
324  E-7
325 X<Y?
326 GTO 12
327 2
328 RCL 08        
329 ST- Y
330 *
331 RCL 12
332 COS
333 RCL 15
334 /
335 X^2
336 *
337 12
338 RCL Y
339 5
340 *
341 -
342 *
343 64
344 -
345 *
346 256
347 ST- Y
348 /
349 STO 14          
350 CLX
351 37
352 *
353 64
354 -
355 *
356 512
357 /
358 4
359 1/X
360 +
361 *
362 RCL 13
363 X^2
364 ST+ X
365 1
366 -
367 RCL 05
368 COS
369 4
370 /
371 *
372 *
373 RCL 13
374 +
375 *
376 RCL 05
377 SIN
378 *
379 RCL 05
380 D-R
381 -
382 RCL 14
383 *
384 RCL 15
385 *
386 RCL 06 
387 *
388 RTN
389 LBL 02
390 1
391 P-R
392 X<>Y
393 RCL 03        
394 X^2
395 *
396 STO 10
397 X^2
398 STO 04
399 RDN
400 P-R
401 RCL 01
402 X^2
403 *
404 STO 08
405 X^2
406 ST+ 04
407 X<>Y
408 RCL 02
409 X^2
410 *
411 STO 09          
412 X^2
413 RCL 04
414 +
415 SQRT
416 ST/ 08
417 ST/ 09
418 ST/ 10
419 RTN
420 LBL 03
421 ENTER
422 SIN
423 STO 07
424 X^2
425 RCL 10
426 *
427 X<>Y
428 RCL 04
429 -
430 SIN
431 ST* 07
432 X^2
433 RCL 11
434 *
435 +
436 RCL 07
437 RCL 09
438 ST+ X
439 *
440 -
441 SQRT
442 RCL 08
443 X<>Y
444 /
445 RTN
446 LBL 14
447 RCL 17
448 +
449 END

 
      ( 564 bytes / SIZE 022 )
 
 

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

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

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

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

Notes:

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

5°)  Minimizing the Sum of 2 Distances
 

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

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

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

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

Data Registers:   R00 thru R35: temp
Flags:  /
Subroutine:  "TGD"  ( cf paragraph 2-b )
 
 
 

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

 
      ( 180 bytes / SIZE 036 )
 
 

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

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

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

        0         ENTER^  ENTER^
    179.51    ENTER^
        0         XEQ "D+D"  >>>>   D = 20001.89956  km
                                        X<>Y  b" ~ -75°37'

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

   40  CHS  ENTER^  ENTER^
      120       ENTER^
       50        XEQ "D+D"  >>>>   D = 18095.49451  km
                                        X<>Y  b" ~ 24°58'
 

Notes:

-The HP41 displays the successive approximations.

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

-This program is very slow with a true HP41, so a good emulator is not superfluous !
-With free42, you get the results almost instantaneously...

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


6°)  Other Formulae for the Great Elliptic Distance

 
     a) Formula of order 1 + Andoyer's Formula of order 2


-This variant employs Andoyer's formula of order 2 to compute the geodesic distance on WGS84 ellipsoid of revolution,
  but the following formula of order 1 to calculate the distances   Δ = MN on the great ellipses:


   With       u = OM / OM = ( x , y , z )  ,  v = ON / ON = ( x' , y' ,z' )   ,   W = angle ( OM , ON )

     and         f = (a-b) / a  &  f ' = (a-c) / c


    Δ  ~  a [ W - {  ( f y2 + f y'2 + f ' z2 + f ' z'2 )  ( W/2 - (1/4) Sin 2.W ) + ( f y.y' + f ' z.z' ) ( Sin W - W Cos W ) } / Sin2 W ]

 
-This formula is an approximation of the integral:       Δ  = §0W [ R2 + ( dR/dµ )2 ] 1/2 dµ 


  with       R / Sin W = 1 / ( A2 + B2 + C2 )1/2   where

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


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

-Line 36 is a three-byte GTO 03
 
 

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

 
      ( 318 bytes / SIZE 015 )
 
 

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

   Where  D is the geodesic distance on the triaxial ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65746 km                                                      ---Execution time = 25s---
                                 LASTX   -       d  = 12138.68425 km

  where  d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )

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

Notes:

-The distance is computed by  D ~ d + Δ1 - Δ2  
-It might also be calculated by  D ~ d . Δ12   ( replace line 35  by  ST/ 14  and line 251 by  *  )
-With example above, the results are the same.

-This program is shorter and slightly faster than the other versions.
-Though the precision is not the best, it's often accurate enough
-However, a formula of order 3 will give a better precision:
 

     b) Formula of order 1 + Andoyer's Formula of order 3


-The same method is used for the great ellipse arc distances, but Andoyer's formula of order 3 is employed to compute the geodesic distance on WGS84 ellipsoid of revolution:



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

-Line 36 is a three-byte GTO 03
 
 

 01 LBL "TGD"
 02 DEG
 03 HR
 04 STO 14
 05 RDN
 06 HR
 07 14.929
 08 STO 11
 09 +
 10 STO 13
 11 RDN
 12 HR
 13 STO 12
 14 X<>Y
 15 HR
 16 ST+ 11
 17 6378.172
 18 STO 05
 19 .07
 20 -
 21 STO 06
 22 21.349686
 23 -
 24 STO 07
 25 XEQ 01
 26 STO 04
 27 RCL 05
 28 RCL 06
 29 +
 30 2
 31 /
 32 STO 05
 33 STO 06
 34 XEQ 01
 35 ST- 04
 36 GTO 03
 37 LBL 01
 38 RCL 11
 39 RCL 12
 40 XEQ 02
 41 STO 08
 42 RCL 02
 43 STO 09          
 44 RCL 03
 45 STO 10
 46 RCL 13
 47 RCL 14
 48 XEQ 02
 49 RCL 08
 50 *
 51 RCL 02
 52 RCL 09
 53 *
 54 STO 01
 55 +
 56 RCL 03
 57 RCL 10
 58 *
 59 STO 08
 60 +
 61 STO 00
 62 RCL 03
 63 X^2
 64 RCL 10
 65 X^2
 66 +
 67 RCL 05
 68 RCL 07
 69 -
 70 ST* 08
 71 STO 10
 72 *
 73 RCL 02
 74 X^2
 75 RCL 09
 76 X^2
 77 +
 78 RCL 05
 79 STO 09
 80 ST/ 10
 81 RCL 06
 82 -
 83 ST* 01
 84 *
 85 +
 86 RCL 00
 87 ENTER
 88 ACOS
 89 D-R
 90 ST* 09
 91 LASTX
 92 SIN
 93 STO 02          
 94 /
 95 ST* 00
 96 -
 97 *
 98 2
 99 /
100 RCL 00
101 RCL 01
102 RCL 08
103 +
104 ST* Y
105 -
106 +
107 RCL 02
108 /
109 RCL 09
110 +
111 RTN
112 LBL 02
113 1
114 P-R
115 X<>Y
116 RCL 07
117 X^2
118 *
119 STO 03
120 RDN
121 P-R
122 RCL 05
123 X^2
124 *
125 STO 01
126 X^2
127 X<>Y
128 RCL 06
129 X^2
130 *
131 STO 02
132 X^2
133 RCL 03
134 X^2
135 +
136 +
137 SQRT
138 ST/ 01
139 ST/ 02
140 ST/ 03
141 RCL 01
142 RTN
143 LBL 03
144 RCL 12          
145 1
146 STO 08
147 P-R
148 LASTX
149 RCL 10
150 -
151 STO 07
152 /
153 R-P
154 X<> 14
155 RCL 08
156 P-R
157 RCL 07
158 /
159 R-P
160 RDN
161 +
162 2
163 ST/ 10
164 /
165 ST- Y
166 RCL 08
167 P-R
168 X^2
169 STO 02
170 STO T
171 RDN
172 X^2
173 STO 01
174 SIGN
175 P-R
176 X^2
177 ST* 01
178 RDN
179 X^2
180 ST* 02
181 -
182 RCL 11
183 RCL 13
184 -
185 2
186 /
187 COS
188 X^2
189 *
190 -
191 ST/ 02
192 RCL 02          
193 RCL 01
194 RCL 08
195 R^
196 STO 01
197 -
198 /
199 ST- 02
200 +
201 X<> 01
202 SQRT
203 ASIN
204 ST+ X
205 D-R
206 STO 00
207 LASTX
208 RCL 08
209 P-R
210 STO 03
211 STO 06
212 RDN
213 STO 11
214 /
215 STO 07
216 ST/ Z
217 X^2
218 ST* 03
219 RCL 08
220 +
221 ENTER
222 R^
223 ST+ 08
224 ST- Z
225 2
226 /
227 STO 09
228 -
229 RCL 01
230 *
231 -
232 RCL 11
233 X^2
234 1.5
235 ST/ Y
236 FRC
237 -
238 RCL 07
239 /
240 RCL 02          
241 *
242 +
243 RCL 02
244 *
245 RCL 06
246 X^2
247 RCL 07
248 ST- 03
249 ST+ X
250 /
251 RCL 03
252 4
253 ST/ 08
254 *
255 RCL 06
256 +
257 RCL 07
258 +
259 STO Z
260 -
261 RCL 01
262 *
263 +
264 RCL 03
265 ST+ X
266 STO 11
267 +
268 RCL 01
269 *
270 +
271 RCL 11
272 -
273 RCL 02
274 *
275 RCL 03
276 RCL 06
277 *
278 STO 03
279 RCL 00
280 ST+ X
281 X^2
282 3
283 ST* Z
284 /
285 STO Z
286 +
287 RCL 08
288 +
289 STO 11          
290 RCL 01
291 *
292 RCL 11
293 ST+ X
294 -
295 RCL 03
296 ST+ Z
297 ST+ Z
298 +
299 RCL 01
300 *
301 +
302 RCL 01
303 *
304 -
305 RCL 10
306 *
307 RCL 08
308 RCL 06
309 RCL 07
310 *
311 STO Z
312 -
313 RCL 01
314 *
315 +
316 RCL 02
317 ST* 09
318 RCL 07
319 ST- 09
320 *
321 -
322 RCL 01
323 *
324 +
325 RCL 02
326 RCL 09
327 *
328 -
329 RCL 10
330 *
331 RCL 02
332 RCL 07
333 /
334 -
335 RCL 01
336 -
337 RCL 10
338 *
339 RCL 00          
340 ST* Y
341 +
342 RCL 05
343 *
344 RCL 04
345 +
346 END

 
      ( 429 bytes / SIZE 015 )
 
 

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

   Where  D is the geodesic distance on the triaxial ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65753 km                                                      ---Execution time = 30s---
                                 LASTX   -       d  = 12138.68432 km

  where  d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )


Notes:

-Surprisingly, the precision is almost as good as the results given by the program listed in §2°)b)
-For example, when the difference in the longitudes is 165°, the error can reach about 83 cm instead of 75 cm.

 

     c) Formula of order 2 + Andoyer's Formula of order 3



-A formula of order 2 is used for the great ellipse arc distances and Andoyer's formula of order 3 is employed to compute the geodesic distance on WGS84 ellipsoid of revolution:


               m = [ y2 ( 2.f + 3.f2 ) + z2 ( 2.f ' + 3.f '2 ) ] / ( 2 Sin2 W )
    With    n  = [ y'2 ( 2.f + 3.f2 ) + z'2 ( 2.f ' + 3.f '2 ) ] / ( 2 Sin2 W )         [ In fact,  2.f + 3.f2 & 2.f ' + 3.f '2  may be replaced by  (a2/b2 - 1)/2  and  (a2/c2 - 1)/2 ]
               p  = [ y.y' ( 2.f + 3.f2 ) + z.z' ( 2.f ' + 3.f '2 ) ] / ( Sin2 W )

    Δ / (a.Sin W)  ~   W / Sin W + ( p2/2 - m - n ) ( W/Sin W - Cos W ) / 2  + (p/2) ( W Cot W - 1 )
                           +  ( m2 + n2 ) [ (13/16) ( W/Sin W - Cos W ) + (1/8) Sin2 W Cos W ) ]
                           + 3 ( p2/2 ) [ (-1/4) W Sin W+ (3/8) ( W/Sin W - Cos W ) ]
                           +  m n  [ (-7/4) W Sin W + (13/8) ( W/Sin W -  Cos W ]
                           +  p ( m + n ) [ (3/4) Sin4 W + ( Cos W ) { (-13/8) ( W/Sin W - Cos W ) + (3/4) Sin2 W  Cos W } ]


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

-Line 36 is a three-byte GTO 03
 
 

 01 LBL "TGD"
 02 DEG
 03 HR
 04 STO 14
 05 RDN
 06 HR
 07 14.92911
 08 STO 11
 09 +
 10 STO 13
 11 RDN
 12 HR
 13 STO 12
 14 X<>Y
 15 HR
 16 ST+ 11
 17 6378.172
 18 STO 05
 19 .07
 20 -
 21 STO 06
 22 21.349686
 23 -
 24 STO 07
 25 XEQ 01
 26 STO 15
 27 RCL 05
 28 RCL 06
 29 +
 30 2
 31 /
 32 STO 05
 33 STO 06
 34 XEQ 01
 35 ST- 15
 36 GTO 03
 37 LBL 01
 38 RCL 11
 39 RCL 12
 40 XEQ 02
 41 STO 08
 42 RCL 02
 43 STO 09
 44 RCL 03
 45 STO 10
 46 RCL 13
 47 RCL 14
 48 XEQ 02
 49 RCL 08
 50 *
 51 RCL 02
 52 RCL 09
 53 *
 54 STO 01
 55 +
 56 RCL 03
 57 RCL 10          
 58 *
 59 STO 08
 60 +
 61 STO 00
 62 RCL 09
 63 X^2
 64 RCL 05
 65 RCL 06
 66 /
 67 X^2
 68 1
 69 -
 70 ST* 01
 71 STO 09
 72 *
 73 RCL 05
 74 RCL 07
 75 /
 76 X^2
 77 1
 78 -
 79 STO 04
 80 ST* 08
 81 RCL 10
 82 X^2
 83 *
 84 +
 85 STO 10
 86 RCL 02
 87 X^2
 88 RCL 09
 89 *
 90 RCL 03
 91 X^2
 92 RCL 04
 93 *
 94 +
 95 STO 02
 96 RCL 01
 97 RCL 08
 98 +
 99 RCL 00
100 ACOS
101 STO 08
102 SIN
103 STO 03
104 X^2
105 ST+ X
106 ST/ 02
107 ST/ 10
108 /
109 STO 09
110 RCL 08
111 D-R
112 STO 08
113 RCL 03
114 ST* 08
115 /
116 STO Y
117 RCL 00
118 -
119 2
120 /
121 STO 04
122 RCL 09          
123 X^2
124 ST+ X
125 STO 01
126 RCL 02
127 -
128 RCL 10
129 -
130 *
131 +
132 X<>Y
133 RCL 00
134 *
135 RCL 09
136 ST* Y
137 -
138 +
139 RCL 04
140 3
141 *
142 RCL 08
143 -
144 RCL 01
145 *
146 .75
147 *
148 +
149 RCL 04
150 13
151 *
152 STO 01
153 RCL 08
154 7
155 *
156 -
157 RCL 02
158 *
159 RCL 10
160 *
161 4
162 /
163 +
164 RCL 01
165 RCL 00
166 RCL 03
167 X^2
168 *
169 ST+ 01
170 3
171 *
172 -
173 RCL 00
174 *
175 RCL 03
176 X^2
177 X^2
178 3
179 *
180 -
181 RCL 09
182 *
183 RCL 02
184 RCL 10          
185 +
186 *
187 2
188 /
189 -
190 RCL 01
191 RCL 02
192 X^2
193 RCL 10
194 X^2
195 +
196 *
197 8
198 /
199 +
200 RCL 03
201 *
202 RCL 05
203 *
204 RTN
205 LBL 02
206 1
207 P-R
208 X<>Y
209 RCL 07
210 X^2
211 *
212 STO 03
213 RDN
214 P-R
215 RCL 05
216 X^2
217 *
218 STO 01
219 X^2
220 X<>Y
221 RCL 06
222 X^2
223 *
224 STO 02
225 X^2
226 RCL 03
227 X^2
228 +
229 +
230 SQRT
231 ST/ 01
232 ST/ 02
233 ST/ 03
234 RCL 01
235 RTN
236 LBL 03
237 RCL 05
238 RCL 07
239 -
240 RCL 05
241 /
242 STO 04
243 RCL 12
244 1
245 STO 08
246 P-R
247 LASTX
248 R^
249 -
250 STO 07          
251 /
252 R-P
253 X<> 14
254 RCL 08
255 P-R
256 RCL 07
257 /
258 R-P
259 RDN
260 +
261 2
262 ST/ 04
263 /
264 ST- Y
265 RCL 08
266 P-R
267 X^2
268 STO 02
269 STO T
270 RDN
271 X^2
272 STO 01
273 SIGN
274 P-R
275 X^2
276 ST* 01
277 RDN
278 X^2
279 ST* 02
280 -
281 RCL 11
282 RCL 13
283 -
284 2
285 /
286 COS
287 X^2
288 *
289 -
290 ST/ 02
291 RCL 02
292 RCL 01
293 RCL 08
294 R^
295 STO 01
296 -
297 /
298 ST- 02
299 +
300 X<> 01
301 SQRT
302 ASIN
303 ST+ X
304 D-R
305 STO 00
306 LASTX
307 RCL 08
308 P-R
309 STO 03
310 STO 06
311 RDN
312 STO 10          
313 /
314 STO 07
315 ST/ Z
316 X^2
317 ST* 03
318 RCL 08
319 +
320 ENTER
321 R^
322 ST+ 08
323 ST- Z
324 2
325 /
326 STO 09
327 -
328 RCL 01
329 *
330 -
331 RCL 10
332 X^2
333 1.5
334 ST/ Y
335 FRC
336 -
337 RCL 07
338 /
339 RCL 02
340 *
341 +
342 RCL 02
343 *
344 RCL 06
345 X^2
346 RCL 07
347 ST- 03
348 ST+ X
349 /
350 RCL 03
351 4
352 ST/ 08
353 *
354 RCL 06
355 +
356 RCL 07
357 +
358 STO Z
359 -
360 RCL 01
361 *
362 +
363 RCL 03
364 ST+ X
365 STO 10
366 +
367 RCL 01
368 *
369 +
370 RCL 10
371 -
372 RCL 02
373 *
374 RCL 03
375 RCL 06          
376 *
377 STO 03
378 RCL 00
379 ST+ X
380 X^2
381 3
382 ST* Z
383 /
384 STO Z
385 +
386 RCL 08
387 +
388 STO 10
389 RCL 01
390 *
391 RCL 10
392 ST+ X
393 -
394 RCL 03
395 ST+ Z
396 ST+ Z
397 +
398 RCL 01
399 *
400 +
401 RCL 01
402 *
403 -
404 RCL 04
405 *
406 RCL 08
407 RCL 06
408 RCL 07
409 *
410 STO Z
411 -
412 RCL 01
413 *
414 +
415 RCL 02
416 ST* 09
417 RCL 07
418 ST- 09
419 *
420 -
421 RCL 01
422 *
423 +
424 RCL 02
425 RCL 09
426 *
427 -
428 RCL 04
429 *
430 RCL 02
431 RCL 07
432 /
433 -
434 RCL 01
435 -
436 RCL 04
437 *
438 RCL 00          
439 ST* Y
440 +
441 RCL 05
442 *
443 RCL 15
444 +
445 END

 
      ( 537 bytes / SIZE 016 )
 
 

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

   Where  D is the geodesic distance on the triaxial ellipsoid

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

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65756 km                                                      ---Execution time = 38s---
                                 LASTX   -        d  = 12138.68432 km

  where  d is the geodesic distance on the ellipsoid of revolution ( WGS 84 )


Note:

-Though the formula of order 2 doesn't give a very good precision for the great elliptic distances,
  the differences between the great elliptic distances is accurately computed !

7°)  Geodesic + Great Elliptic Distance



-With a triaxial ellipsoid, we can use programs listed above (  cf §1) §2b) §3b) §4) )
-But with an ellipsoid of revolution, a shoter & faster program may also be written.

-The routine hereunder employs Andoyer's formula of order 3 for the geodesic distance
  and a formula of order 2 for the great elliptic distance ( cf "Great Elliptic Arc Length for the HP41" )


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



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

 
      ( 310 bytes / SIZE 012 )
 
 

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

   Where  D = geodesic distance & d = great elliptic distance

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

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

Notes:

-With the great elliptic distance, the maximum error = 3 cm
-Of course, the precision would be better with the 3rd order formula ( cf "Great Elliptic Arc Length for the HP41" )


Conclusion:

-For a WGS84 ellipsoid, the more accurate results are obtained by Vincenty's formula ( cf "Terrestrial Geodesic Distance(1) for the HP41" paragraph 2 )
-Unfortunately, the iterations take a long time with an HP41.
-So "TGD" listed in §3-a) is often enough accurate.

-Likewise, with a triaxial ellipsoid model of the Earth, the program listed in §4 above is the most accurate, but very slow.
-So, "TGD" listed in §3-b) is satisfactory in many cases... it's up to you to judge !



Remark:

-The values of a & b  that are used in all these programs are those given by WGS84.
-IERS2003 suggests other numbers:  mean equatorial radius = 6378.13661 km & flattening = 1/298.256421 so, polar radius = 6356.751868 km

-For the triaxial ellipsoid, recent evaluations suggest:  equatorial flattening = 1/91026
-So it yields::

   a = 6378.171645 km
   b = 6378.101575 km
   c = 6356.751868 km



References:

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

-Many thanks to Professor Richard H. Rapp who sent me reference [3] !