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 Andoyer's formula of order 3

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

 4°)  With Vincenty's Formula
 5°)  Minimizing the Sum of 2 Distances
 6°)  A Slightly Different Method

   a)  With Andoyer's formula of order 2
   b)  With Andoyer's formula of order 3
   c)  Almost With Andoyer's formula of order 3
   d)  Geodetic or Geocentric Longitudes & Latitudes ?
 

-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    µ = 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 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 Andoyer's Formula of order 3
 

     a) Andoyer's Formula ( Ellipsoid of Revolution )
 

-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  &  B , B' = parametric latitudes     a = 6378.137 km  ,  b = 6356.752314 km  ( WGS84 )  ,  f = (a-b) / a

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

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

 
     ( 346 bytes / SIZE 015 )
 
 

      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 = 16s---
 
 

Notes:

-"AND" is slightly shorter & faster than "SOD"
-The formulas are quite equivalent, so the differences in the results are only due to roundoff-errors.

-It is certainly possible to express the formula without calculating the parametric latitudes, but I've not found how.

-In order to calculate the parametric latitudes, we can also use TAN & ATAN instead of P-R & R-P
-But the HP41 gives TAN (-90°) = +9.999999999 E99 instead of -9.999999999 E99
-So, we must use a trick to take this case into account:

-We can replace lines 01 to 32 by
 
 

 01 LBL "AND"
 02 DEG
 03 HR
 04 X<> T
 05 HMS-
 06 HR
 07 2
 08 /
 09 SIN
 10 X^2
 11 STO 00
 12 RDN
 13 HR
 14 X<Y?
 15 X<>Y
 16 TAN
 17 PI
 18 SIGN
 19 298.25722
 20 STO 04
 21 ST+ 04
 22 1/X
 23 -
 24 STO 05
 25 * 
 26 ATAN
 27 X<>Y
 28 SIGN
 29 LASTX
 30 ABS
 31 TAN
 32 RCL 05
 33 *
 34 *
 35 ATAN

 
-It costs 2 extra bytes but it saves almost 1 second with an HP41 !
-However, this method cannot be used with free42 because 90 TAN gives OUT OF RANGE
 

     b) Triaxial Ellipsoid
 

-The following "TGD" employs Andoyer's formula of order 3 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 STO 09
153 *
154 R-P
155 X<> 04
156 1
157 P-R
158 RCL 09
159 *
160 R-P
161 RDN
162 ST- Z
163 +
164 COS
165 STO 05        
166 X^2
167 ST+ X
168 CHS
169 X<>Y
170 ABS
171 ENTER
172 D-R
173 STO 10
174 SIGN
175 P-R
176 X<>Y
177 ST* 05
178 *
179 ST* Y
180 +
181 RCL 10
182 +
183 4
184 /
185 RCL 07
186 CHS
187 RCL 06
188 ST+ Y
189 ST+ X
190 /
191 STO 09
192 *
193 RCL 05
194 -
195 RCL 10
196 -
197 RCL 09
198 *
199 RCL 10
200 +
201 RCL 06
202 *
203 RTN
204 LBL 02
205 1
206 P-R
207 X<>Y
208 RCL 03
209 X^2
210 *
211 STO 10
212 X^2
213 X<> Z
214 X<>Y
215 P-R
216 RCL 01
217 X^2
218 *
219 STO 08
220 X^2
221 X<>Y
222 RCL 02
223 X^2
224 *
225 STO 09
226 X^2
227 +
228 +
229 SQRT
230 ST/ 08
231 ST/ 09
232 ST/ 10
233 RTN
234 LBL 03
235 ENTER
236 SIN
237 STO 07
238 X^2
239 RCL 10
240 *
241 X<>Y
242 RCL 04
243 -
244 SIN
245 ST* 07
246 X^2
247 RCL 11        
248 *
249 +
250 RCL 07
251 RCL 09
252 *
253 -
254 SQRT
255 RCL 08
256 X<>Y
257 /
258 RTN
259 LBL 14
260 RCL 14
261 RCL 12
262 -
263 2
264 /
265 SIN
266 X^2
267 STO 00
268 RCL 13
269 1
270 P-R
271 RCL 03
272 CHS
273 RCL 01
274 STO 16
275 ST+ Y
276 /
277 STO 04
278 SIGN
279 LASTX
280 -
281 STO 05
282 /
283 R-P
284 X<> 15
285 1
286 P-R
287 RCL 05
288 /
289 R-P
290 RDN
291 +
292 2
293 ST/ 04
294 /
295 ST- Y
296 1
297 P-R
298 X^2
299 STO 02
300 X<>Y
301 X^2
302 STO 01
303 RDN
304 X<>Y
305 1
306 P-R
307 X^2
308 ST* 01
309 RDN
310 X^2
311 ST* 02
312 STO Z
313 -
314 RCL 00
315 *
316 +
317 ST/ 02
318 STO 05
319 SQRT
320 ASIN
321 ST+ X
322 STO 03
323 D-R
324 STO 00
325 1
326 RCL 05        
327 -
328 ST/ 01
329 RCL 02
330 RCL 01
331 ST- 02
332 +
333 STO 01
334 RCL 03
335 1
336 P-R
337 STO 06
338 X<>Y
339 STO 05
340 /
341 STO 07
342 ENTER
343 X^2
344 RCL 00
345 *
346 -
347 RCL 00
348 X^2
349 STO 09
350 *
351 STO 03
352 3
353 *
354 RCL 00
355 RCL 09
356 *
357 STO 10
358 .75
359 STO 15
360 /
361 STO 13
362 -
363 RCL 05
364 RCL 06
365 *
366 STO 11
367 LASTX
368 X^2
369 *
370 RCL 11
371 ENTER
372 X^2
373 STO 08
374 *
375 +
376 RCL 05
377 ST* 08
378 X^2
379 X^2
380 STO 12
381 RCL 11
382 *
383 +
384 ST+ X
385 +
386 RCL 11
387 3
388 *
389 RCL 15
390 *
391 -
392 RCL 00
393 4
394 /
395 -
396 RCL 01
397 *
398 RCL 03
399 5
400 *
401 -
402 RCL 13
403 ST+ X
404 +
405 RCL 00
406 RCL 11        
407 +
408 2
409 /
410 +
411 RCL 05
412 RCL 05
413 RCL 12
414 *
415 -
416 RCL 08
417 -
418 ST+ X
419 RCL 10
420 RCL 05
421 X^2
422 /
423 STO 12
424 RCL 06
425 *
426 STO 13
427 4
428 *
429 -
430 RCL 09
431 RCL 05
432 /
433 STO 14
434 3
435 *
436 +
437 RCL 00
438 RCL 06
439 *
440 STO 08
441 -
442 RCL 06
443 RCL 11
444 *
445 RCL 15
446 ST+ X
447 STO 06
448 *
449 -
450 RCL 02
451 *
452 +
453 RCL 01
454 *
455 RCL 03
456 ST+ X
457 +
458 RCL 10
459 RCL 15
460 /
461 -
462 RCL 13
463 6
464 *
465 RCL 14
466 5
467 *
468 -
469 RCL 08
470 +
471 RCL 02
472 *
473 +
474 RCL 01
475 *
476 RCL 11
477 2
478 /
479 RCL 00
480 -
481 RCL 12
482 -
483 RCL 01
484 *
485 RCL 05
486 X^2
487 RCL 06        
488 /
489 .5
490 -
491 RCL 05 
492 *
493 RCL 02
494 *
495 +
496 RCL 11
497 -
498 RCL 12
499 +
500 RCL 00
501 +
502 RCL 02
503 *
504 RCL 14
505 RCL 13
506 -
507 ST+ X
508 +
509 RCL 02
510 *
511 +
512 RCL 04
513 *
514 RCL 00
515 RCL 11
516 +
517 4
518 /
519 RCL 07
520 RCL 09
521 *
522 STO 13
523 -
524 RCL 01
525 *
526 RCL 02
527 ST* 05
528 ST* 11
529 RCL 14
530 *
531 -
532 RCL 13
533 +
534 RCL 01
535 *
536 +
537 RCL 14
538 RCL 11
539 2
540 /
541 -
542 RCL 02
543 *
544 +
545 RCL 04
546 *
547 RCL 00
548 RCL 01
549 *
550 -
551 RCL 05
552 -
553 RCL 04
554 *
555 RCL 00
556 +
557 RCL 16
558 *
559 STO Y
560 RCL 17
561 *
562 END

 
      ( 662 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 = 56s---
                                             X<>Y  d  = 12138.68432 km

Note:

-If you want to use  TAN & ATAN  instead of  P-R & R-P  to calculate the parametric latitudes, replace lines 268 to 290 by
 
 

268 RCL 13
269 RCL 15
270 X<Y?
271 X<>Y
272 TAN
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 05
284 *
285 ATAN
286 X<>Y 
287 SIGN
288 LASTX
289 ABS
290 TAN
291 RCL 05
292 *
293 *
294 ATAN

 
-Just 3 extra bytes, but almost 1 second saved.
 

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 or 3-b above, using 3rd order formulas )
 
 
 

 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 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.89955  km
                                        X<>Y  b" ~ -75°38'

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

   40  CHS  ENTER^  ENTER^
      120       ENTER^
       50        XEQ "D+D"  >>>>   D = 18095.49447  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 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 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 !

Variant:

-The following routine may be used to find the minimum distance from several values of 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 RCL 26
 23 HR
 24 SIN
 25 RCL 24
 26 HR
 27 SIN
 28 +
 29 X<>Y
 30 R-P
 31 RDN
 32 HMS
 33 STO 28        
 34 X<>Y
 35 RCL 23
 36 HR
 37 +
 38 HMS
 39 STO 27
 40 X<>Y
 41 " B0="
 42 ARCL X
 43 PROMPT
 44 LBL 10
 45 " B=?"
 46 PROMPT
 47 STO 29        
 48 RCL 27
 49 X<>Y
 50 RCL 25
 51 RCL 26
 52 XEQ "TGD"
 53 STO 30
 54 RCL 23
 55 RCL 24
 56 RCL 27
 57 RCL 29
 58 XEQ "TGD"
 59 RCL 30        
 60 +
 61 " D="
 62 ARCL X
 63 PROMPT
 64 GTO 10
 65 END

 
      ( 116 bytes )
 
 

      STACK      INPUTS
           T      L ( ° ' " )
           Z      b ( ° ' " )
           Y      L' ( ° ' " )
           X      b' ( ° ' " )

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

   40  CHS  ENTER^  ENTER^
      120       ENTER^
       50        XEQ "D+D"  >>>>   B0=24.1727     i-e  24°17'27"

                               R/S   >>>>    B=?                   if you choose for instance:  b" = 12°  17°  22°  27°  32°

                        12   R/S   >>>>   D=18114.88342      R/S     >>>>    B=?
                        17   R/S   >>>>   D=18102.74082      R/S     >>>>    B=?
                        22   R/S   >>>>   D=18096.49637      R/S     >>>>    B=?
                        27   R/S   >>>>   D=18095.96331      R/S     >>>>    B=?
                        32   R/S   >>>>   D=18101.15912

-Store these 5 results in R01 thru R05 and if you use the programs listed in "Extrema for the HP41" paragraph 1-d)

    5   ENTER^
   22  XEQ "DEX5"  >>>>   b" =  24.96999 = 24°58'12"
                                X<>Y  D = 18095.49452 km

-With a true HP41, it will be less slow than the first "D+D"
-However, it will not be always easy to estimate the proper - and small - interval for b" for nearly antipodal points...
 

6°)  A Slightly Different Method
 

     a) With Andoyer's Formula of Order 2
 

-In the 3 following programs, intead of calculating   D ~ (A/A') d , we compute D ~  A + ( d - A' )
-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.

-The precision is almost the same, but the routines are faster !
 

Formula:

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

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

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

 
      ( 437 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.65756 km                                                      ---Execution time = 28s---
 

Notes:

-The root-mean-square error is about 3 cm if  80°  < | L' - L | < 120°  and the maximum error is about 16 cm
-The root-mean-square error is about 6 cm if 110° < | L' - L | < 150°  and the maximum error is about 78 cm
-When the difference in the longitudes is 165°, the error can reach about 4 meters.

-Registers L & R06 contain the great elliptic arc distance.
 

     b) With Andoyer's Formula of Order 3
 

-In this variant, the WGS84 geodesic distance is computed with Andoyer's formula of order 3 ( the other distances with formulas of order 2 )
 

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

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

 
      ( 628 bytes / SIZE 018 )
 
 

      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.65757 km                                                      ---Execution time = 35s---
 

Notes:

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

-Registers L & R17 contain the great elliptic arc distance.
 

     c) Almost With Andoyer's Formula of Order 3
 

-We can save about 100 bytes if we only use the larger terms in the 3rd order formula:
-The precision for very long geodesics are of the same order as the program above.
 
 

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

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

 
      ( 530 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.65756 km                                                      ---Execution time = 31s---
 

Notes:

-The root-mean-square error is about 3 cm if  80°  < | L' - L | < 120°  and the maximum error is about 15 cm
-The root-mean-square error is about 6 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 78 cm.
-I've found these results with free42binary and with latitudes between -87° & +87°

-Registers L & R06 contain the great elliptic arc distance.
 

           d)  Geodetic or Geocentric Longitudes & Latitudes ?
 

-This program employs the same formulas as above, but you can choose the semi-axes of the Earth, and geodetic or geocentric longitudes and latitudes.
 
 

Data Registers:   (•) R00 = L0                                        ( Registers R00 thru R04 are to be initialized IF Flag F00 is set )

                               (•) R01 = a                    R04 to R16: temp
                               (•) R02 = b
                               (•) R03 = c

Flags:    F00-F03-F04

                               If  F00 is clear, "TGD" uses a = 6378.172 km  ,  b = 6378.102 km  ,  c = 6356.752314 km ,  L0 = -14°92911
                               If  F00 is set , initialize  R00 to R03

                               If  F04 is clear & F03 is clear , longitudes and latitudes are geodetic
                               If  F04 is clear & F03 is set , longitudes are geocentric but latitudes are geocentric

                               If  F04 is set, longitudes and latitudes are geocentric.

Subroutines: /
 

-Line 190 is a three-byte GTO 14
 
 

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

 
         ( 592 bytes / SIZE 017 )
 
 

      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 )

     • If CF 00  CF 03  CF 04  geodetic longitudes & latitudes

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.65756 km
                                           LASTX  d = 12138.65876 km = R06 = Great elliptic arc distance

     • If CF 00  SF 03  CF 04   geocentric longitudes & geodetic latitudes

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12138.70239 km
                                           LASTX  d = 12138.70359 km = R06 = Great elliptic arc distance

     • If CF 00  SF 04   geocentric longitudes & latitudes

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "TGD"   >>>>   D = 12157.22639 km
                                           LASTX  d = 12157.22759 km = R06 = Great elliptic arc distance
 

Notes:

-The latitudes given on the web are always geodetic, but I don't really know if the longitudes are geocentric or geodetic !
-It's the same for an ellipsoid of revolution, but not for a triaxial ellipsoid...
 
 

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 !