hp41programs

Geodesic2

Terrestrial Geodesic Distance(2) for the HP-41


Overview
 

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

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

 3°)  With Vincenty's Formula
 4°)  Minimizing the Sum of 2 Distances
 

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

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

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

0°)  Method
 

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

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

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

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

    D = (A/A') d
 

1°)  With Andoyer's Formula of order 2
 

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

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

   OP = ( Sin W ) / ( 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
 

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

-Finally, the geodesic distance d on the WGS84 ellipsoid is calculated by 2nd order Andoyer's formula
 and the geodesic distance D on the triaxial ellipsoid is approximated by   D = (A/A') d
 

Andoyer's Formula of order 2:
 

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

              S = 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
 

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

-Line 38 is a three-byte GTO 14
 
 

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

 
      ( 471 bytes / SIZE 016 )
 
 

      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 "TGD"   >>>>   D = 12138.65748 km                                                      ---Execution time = 47s---
                                             X<>Y  d  = 12138.68425 km

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

Notes:

-Longitudes are geodetic and positive East.
-If the points are not nearly antipodal, the precision is of the order of one meter.

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

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

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

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

2°)  With Sodano's Formula of order 3
 

     a) Sodano's Formula ( Ellipsoid of Revolution )
 

-"SOD" uses a 3rd order formula in the flattening to calculate the geodesic distance on an ellipsoid of revolution ( WGS84 ) - cf references [2] & [3]

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

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

          +  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    µ = 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 l  where  l = 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 )
 
 

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

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

 
     ( 374 bytes / SIZE 014 )
 
 

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

Notes:

-If  | L' - L | < 120° , the RMS error is about 1.4 mm and the maximum error about 9 mm
-If  | L' - L | < 170° , -----------------------  12 mm  ----------------------------- 22 cm
-If  | L' - L | < 175° , -----------------------   5  cm   -----------------------------  2 m

-When Vincenty's iterative method requires more than 3 iterations, even a 3rd order formula may produce an error of 2 meters:
-For example, with  Long1 = 0° , Lat1 = 2° &  Long2 = 175° , Lat2 = -4°

   "SOD"  gives  19434.163468 km  instead of  19434.161314 km  returned by Vincenty's formula !

-These estimations were obtained with free42.
-With an HP41 - which works with 10 digits - the errors may be slightly larger...
 

     b) Triaxial Ellipsoid
 

-This variant employs Sodano's formula for the last calculation and Andoyer's formula of order 2 to estimate  (A/A')
 

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

-Line 38 is a three-byte GTO 14
 
 

 01 LBL "TGD"
 02 DEG
 03 HR
 04 STO 15        
 05 RDN
 06 HR
 07 14.92911
 08 STO 00
 09 +
 10 STO 14
 11 RDN
 12 HR
 13 STO 13
 14 X<>Y
 15 HR
 16 RCL 00
 17 +
 18 STO 12
 19 6378.172
 20 STO 01
 21 .07
 22 -
 23 STO 02
 24 21.349686
 25 -
 26 STO 03
 27 XEQ 01
 28 STO 17
 29 RCL 01
 30 RCL 02
 31 +
 32 2
 33 /
 34 STO 01
 35 STO 02
 36 XEQ 01
 37 ST/ 17
 38 GTO 14
 39 LBL 01
 40 RCL 12
 41 RCL 13
 42 XEQ 02
 43 RCL 08
 44 STO 05
 45 RCL 09
 46 STO 06
 47 RCL 10
 48 STO 07
 49 RCL 14
 50 RCL 15
 51 XEQ 02
 52 RCL 05
 53 RCL 08
 54 *
 55 RCL 06
 56 RCL 09
 57 *
 58 RCL 07
 59 RCL 10
 60 *
 61 +
 62 +
 63 STO 04
 64 RCL 01
 65 ST/ 05
 66 ST/ 08
 67 RCL 02
 68 ST/ 06
 69 ST/ 09
 70 RCL 03
 71 ST/ 07
 72 ST/ 10
 73 RCL 05
 74 X^2
 75 RCL 06
 76 X^2
 77 RCL 07
 78 X^2
 79 +
 80 +
 81 STO 11
 82 RCL 08
 83 X^2
 84 RCL 09
 85 X^2
 86 RCL 10
 87 X^2
 88 +
 89 +
 90 X<> 10        
 91 RCL 07 
 92 *
 93 RCL 06
 94 RCL 09
 95 *
 96 RCL 05
 97 RCL 08
 98 *
 99 +
100 +
101 STO 09
102 ST+ 09
103 RCL 04
104 RCL 11
105 *
106 -
107 ST+ X
108 RCL 04
109 RCL 09
110 *
111 RCL 10
112 -
113 RCL 11
114 RCL 04
115 X^2
116 *
117 -
118 RCL 11
119 RCL 04
120 ACOS
121 STO 04
122 SIN
123 STO 08
124 ST/ Z
125 *
126 +
127 X#0?
128 /
129 ATAN
130 2
131 /
132 90
133 STO 06
134 X<>Y
135 X<0?
136 +
137 STO 05
138 XEQ 03
139 X<> 06
140 RCL 05
141 +
142 XEQ 03
143 STO 07
144 RCL 05
145 CHS
146 ST+ 04
147 1
148 P-R
149 RCL 07
150 RCL 06
151 /
152 X^2
153 STO 09
154 *
155 R-P
156 X<> 04
157 1
158 P-R
159 RCL 09
160 *
161 R-P
162 RDN
163 ST- Z
164 +
165 COS
166 STO 05
167 X^2
168 ST+ X
169 X<>Y
170 ABS
171 ENTER
172 D-R
173 STO 10
174 SIGN
175 P-R
176 X<>Y
177 RCL 10        
178 /
179 STO 11
180 *
181 ST* Y
182 -
183 15
184 *
185 1
186 +
187 4
188 /
189 RCL 07
190 CHS
191 RCL 06
192 ST+ Y
193 ST+ X
194 /
195 STO 09
196 *
197 RCL 05
198 3
199 *
200 RCL 11
201 *
202 -
203 RCL 09
204 ST* Y
205 -
206 RCL 10
207 ST* Y
208 +
209 RCL 06
210 *
211 RTN
212 LBL 02
213 1
214 P-R
215 X<>Y
216 RCL 03
217 X^2
218 *
219 STO 10
220 X^2
221 X<> Z
222 X<>Y
223 P-R
224 RCL 01
225 X^2
226 *
227 STO 08
228 X^2
229 X<>Y
230 RCL 02
231 X^2
232 *
233 STO 09
234 X^2
235 +
236 +
237 SQRT
238 ST/ 08
239 ST/ 09
240 ST/ 10
241 RTN
242 LBL 03
243 ENTER
244 SIN
245 STO 07
246 X^2
247 RCL 10
248 *
249 X<>Y
250 RCL 04
251 -
252 SIN
253 ST* 07
254 X^2
255 RCL 11
256 *
257 +
258 RCL 07
259 RCL 09
260 *
261 -
262 SQRT
263 RCL 08        
264 X<>Y
265 /
266 RTN
267 LBL 14
268 RCL 14
269 RCL 12
270 -
271 STO 07
272 RCL 13
273 RCL 03
274 CHS
275 RCL 01
276 STO 16
277 ST+ Y
278 /
279 STO 04
280 SIGN
281 LASTX
282 -
283 STO 03
284 P-R
285 LASTX
286 /
287 R-P
288 X<>Y
289 STO 05
290 1
291 P-R
292 STO 01
293 X<>Y
294 STO 02
295 RCL 15
296 RCL 03
297 P-R
298 LASTX
299 /
300 R-P
301 X<>Y
302 STO 06
303 1
304 P-R
305 STO 08
306 X<>Y
307 STO 09
308 RCL 02
309 *
310 STO 03
311 RCL 01
312 ST* 09
313 RCL 08
314 ST* 02
315 *
316 STO 10
317 RCL 07
318 2
319 STO 14
320 /
321 SIN
322 X^2
323 *
324 RCL 06
325 RCL 05
326 -
327 RCL 14
328 /
329 SIN
330 X^2
331 +
332 SQRT
333 ASIN
334 ST+ X
335 D-R
336 STO 13
337 LASTX
338 1
339 P-R
340 STO 11
341 X<>Y
342 STO 12
343 RCL 07
344 SIN
345 RCL 10
346 *
347 X<>Y
348 /
349 X^2
350 STO 10        
351 SIGN
352 LASTX
353 -
354 STO 05
355 3
356 *
357 4
358 -
359 RCL 05
360 *
361 STO 02
362 RCL 03
363 ST* Y
364 +
365 RCL 11
366 RCL 12
367 /
368 STO 08
369 *
370 RCL 12
371 ST+ X
372 /
373 RCL 10
374 X^2
375 RCL 05
376 *
377 6
378 /
379 -
380 3
381 RCL 05
382 ST+ X
383 -
384 RCL 05
385 X^2
386 STO 00
387 *
388 RCL 05
389 -
390 RCL 08
391 X^2
392 *
393 RCL 14
394 /
395 +
396 RCL 10
397 RCL 03
398 X^2
399 STO 15
400 *
401 RCL 12
402 X^2
403 STO 09
404 ST+ X
405 /
406 +
407 RCL 13
408 *
409 RCL 02
410 RCL 05
411 -
412 RCL 14
413 +
414 4
415 /
416 STO 02
417 RCL 05
418 *
419 RCL 08
420 *
421 +
422 RCL 02
423 RCL 03
424 *
425 RCL 12
426 /
427 -
428 RCL 13
429 X^2
430 STO 07
431 *
432 RCL 00
433 RCL 14
434 /
435 RCL 15        
436 8
437 *
438 STO 06
439 +
440 RCL 05
441 *
442 RCL 00
443 -
444 STO 02
445 RCL 06
446 ST- 02
447 +
448 RCL 11
449 RCL 12
450 *
451 STO 01
452 *
453 RCL 02
454 RCL 13
455 *
456 +
457 16
458 STO 02
459 /
460 -
461 RCL 10
462 RCL 13
463 *
464 4
465 *
466 RCL 10
467 RCL 01
468 ST* Y
469 +
470 -
471 RCL 00
472 *
473 RCL 11
474 X^2
475 *
476 RCL 02
477 /
478 +
479 RCL 05
480 RCL 01
481 ST* Y
482 +
483 RCL 10
484 RCL 13
485 *
486 1.5
487 STO 08
488 *
489 -
490 RCL 14
491 /
492 RCL 03
493 RCL 09
494 *
495 RCL 12
496 *
497 +
498 RCL 03
499 *
500 RCL 05
501 *
502 RCL 11
503 *
504 +
505 RCL 09
506 X^2
507 RCL 00
508 ST* Y
509 -
510 RCL 15
511 +
512 RCL 03
513 *
514 RCL 12
515 *
516 RCL 01
517 RCL 05
518 *
519 X^2
520 LASTX
521 *
522 6
523 /
524 +
525 RCL 14        
526 /
527 +
528 RCL 03
529 RCL 12
530 *
531 X^2
532 LASTX
533 *
534 RCL 08
535 /
536 -
537 RCL 04
538 *
539 RCL 05
540 RCL 11
541 *
542 RCL 03
543 -
544 RCL 07
545 *
546 RCL 10
547 *
548 RCL 12
549 ST+ X
550 /
551 +
552 RCL 13
553 RCL 01
554 RCL 11
555 X^2
556 *
557 ST+ X
558 -
559 RCL 00
560 ST* Y
561 RCL 06
562 -
563 RCL 01
564 *
565 +
566 RCL 02
567 /
568 +
569 RCL 01
570 RCL 03
571 *
572 RCL 05
573 *
574 RCL 11
575 *
576 RCL 14
577 /
578 +
579 RCL 04
580 *
581 RCL 03
582 RCL 12
583 *
584 +
585 RCL 01
586 RCL 13
587 +
588 RCL 05
589 *
590 RCL 14
591 /
592 -
593 RCL 04
594 *
595 RCL 13
596 +
597 RCL 16
598 *
599 STO Y
600 RCL 17
601 *
602 END

 
      ( 699 bytes / SIZE 018 )
 
 

      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 "TGD"   >>>>   D = 12138.65755 km                                                      ---Execution time = 59s---
                                             X<>Y  d  = 12138.68432 km
 

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

4°)  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:  "DGT"  ( cf paragraph 2-b above, using Sodano's 3rd order formula )
 
 
 

 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 "DGT"
 36 STO 31
 37 RCL 23
 38 RCL 24
 39 RCL 27
 40 RCL 35
 41 XEQ "DGT"
 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 RCL 13
 98 PI
 99 R-D
100 /
101 RCL 28
102 RCL 30
103 HMS-
104 HR
105 ABS
106 X>Y?
107 GTO 10
108 RCL 30
109 RCL 31
110 CLD
111 END

 
      ( 181 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                    ---Execution time = 43m02s---
                                        X<>Y  b" ~ 75°34'

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

   40  CHS  ENTER^  ENTER^
      120       ENTER^
       50        XEQ "D+D"  >>>>   D = 18095.49438  km                    ---Execution time = 40m56s---
                                        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 !
-"DGT" is called 44 ( respectively 42 ) times in the examples above.
-With free42, you get the results 20001.89954 & 18095.49449 almost instantaneously...

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

References:

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

-Many thanks to Professor Richard H. Rapp who sent me reference [3] to find the typo in Sodano's formula !