hp41programs

Geod&Lox3axial

Geodesics & Loxodromes on a Triaxial Ellipsoid for the HP-41


Overview
 

 0°)  2 M-Code Routines
 1°)  Geodesic Distance
 2°)  Loxodromes
 
 

0°)  2 M-Code Routines
 

-In the following programs, the inegrals are evaluated with Gauss-Legendre 7-point formula
-GX7 & GW7 return the abscissas & weights ( the 1st abscissa is zero, so I didn't code it ):

   1  XEQ "GX7"  ->  0.4058451514
   2  XEQ "GX7"  ->  0.7415311856
   3  XEQ "GX7"  ->  0.9491079123

   0  XEQ "GW7"  ->  0.4179591837
   1  XEQ "GW7"  ->  0.3818300505
   2  XEQ "GW7"  ->  0.2797053915
   3  XEQ "GW7"  ->  0.1294849662
 

0B7  "7"
018  "X"
007  "G"
2A0  SETDEC
0F8  C=X
2FC  RCR 13
27E  C=C-1 MS
10E  A=C ALL
04E  C=0 ALL
35C  PT=12
1BE  A=A-1 MS
063  JNC+12d
110
010
150
210
110  C=4.058451514
150
050
150
050
110
0C3  JNC+24d
1BE  A=A-1 MS
063  JNC+12d
1D0
110
050
150
0D0  C=7.415311856
050
050
210
150
190
05B  JNC+11d
250
110
250
050
010  C=9.491079123
1D0
250
050
090
0D0
266  C=C/10
0E8  X=C
3E0  RTN
0B7  "7"
017  "W"
007  "G"
2A0  SETDEC
0F8  C=X
2FC  RCR 13
10E  A=C ALL
04E  C=0 ALL
35C  PT=12
1BE  A=A-1 MS
063  JNC+12d
110
050
1D0
250
150  C=4.179591837
250
050
210
0D0
1D0
12B  JNC+37d
1BE  A=A-1 MS
063  JNC+12d
0D0
210
050
210
0D0  C=3.818300505
010
010
150
010
150
0C3  JNC+24d
1BE  A=A-1 MS
063  JNC+12d
090
1D0
250
1D0
010  C=2.797053915
150
0D0
250
050
150
05B  JNC+11d
050
090
250
110
210  C=1.294849662
110
250
190
190
090
266  C=C/10
0E8  X=C
3E0  RTN

Note:

-The content of X-register is overwritten.
 

1°)  Geodesic Distance
 

>>>  This program is just a slight modification of the last programs listed in "Geodesics on a Triaxial Ellipsoid for the HP41"

 "DGT3" uses an approximate method to evaluate the geodesic distance between 2 points M & N on the Earth
  The error is about a few centimeters for long geodesics, except for nearly antipodal points ( do not use this program in this case )
 

-Let a point P ( on the ellipsoid ) with ( OM , OP' ) = µ , ( OP' , OP ) = h    ,    P' = orthogonal projection of P on the plane (OMN)

    OP = x u + y v + z w       where  u = OM / || OM ||   ,  v = ON / || ON ||  ,  w = u x v          ( x = cross-product )

-We have, with  R = OP  and  W = ( OM , ON )

    x = R Cos h Sin ( W - µ ) / Sin W
    y = R Cos h Sin µ / Sin W
    z = R Sin h / Sin W

-Writing that P(X,Y,Z)  is on the ellipsoid  X2/a2 + Y2/b2 + Z2/c2 = 1,  we find that

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

      A = [ x1 Cos h Sin ( W - µ ) + x2 Cos h Sin µ + ( y1z2 - y2z1 ) Sin h ] / a
      B = [ y1 Cos h Sin ( W - µ ) + y2 Cos h Sin µ + ( z1x2 - z2x1 ) Sin h ] / b                       where  u ( x1 , y1 , z1 )  &  v ( x2 , y2 , z2 )
      C = [ z1 Cos h Sin ( W - µ ) + z2 Cos h Sin µ + ( x1y2 - x2y1 ) Sin h ] / c
 

-We search an approximation of the function  h(µ) under the form:

   h(µ) = h1 Sin ( 180° µ / W ) + h2 Sin ( 360° µ / W ) + h3 Sin ( 540° µ / W )                 ( the first terms of a Fourier series )

-And the corresponding integral  s = §0W [ R2 Cos2 h  + R2 ( dh/dµ )2 + ( dR/dµ )2  ] 1/2 dµ

      =  ( Sin W )    §0W [ ( Cos2 h ) / T + (1/T) ( dh/dµ )2 + ( dr/dµ )2  ] 1/2 dµ    is evaluated by a Gaussian quadrature ( Gauss-Legendre 7-point formula )

  where   T = ( A2 + B2 + C2 )

    r = 1/sqrt(T)                  dr/dµ = r/µ + ( r/h ) ( dh/dµ )

    r = 1 / ( A2 + B2 + C2 )1/2

  the partial derivatives

  r/µ = ( -1/sqrt(T3) ) [ (A/a) { -x1 Cos h Cos ( W - µ ) + x2 Cos h Cos µ }
                                   + (B/b) { -y1 Cos h Cos ( W - µ ) + y2 Cos h Cos µ }
                                   + (C/c) { -z1 Cos h Cos ( W - µ ) + z2 Cos h Cos µ } ]
  and

  r/h = ( -1/sqrt(T3) ) [ (A/a) { -x1 Sin h Sin ( W - µ ) - x2 Sin h Sin µ + ( y1z2 - y2z1 ) Cos h }
                                    + (B/b) { -y1 Sin h Sin ( W - µ ) - y2 Sin h Sin µ + ( z1x2 - z2x1 ) Cos h }
                                    + (C/c) { -z1 Sin h Sin ( W - µ ) - z2 Sin h Sin µ + ( x1y2 - x2y1 ) Cos h } ]
 

-With  H = W / 1000  , we calculate the following lengths

    A = f(0,0,0)
    B = f(H,0,0)
    C = f(-H,0,0)
    D = f(0,H,0)                     f = length s corresponding to h1 = -H , 0 , +H ;  h2 = -H , 0 , +H ;  h3 = -H , 0 , +H
    E = f(0,-H,0)
    F = f(0,0,H)
    G = f(0,0,-H)
 

-Then, the extremum is estimated by the formula:

   s ~ A - (B-C)2 / [ 8.(B+C-2.A) ] - (D-E)2 / [ 8.(D+E-2.A) ] - (F-G)2 / [ 8.(F+G-2.A) ]
 
 

Data Registers:           •  R00 = L0               ( If SF 00 , Registers R00 thru R03 are to be initialized before executing "DGT3" )

                                      •  R01 = a                  R04 to R34: temp
                                      •  R02 = b
                                      •  R03 = c

Flags:  F00-F01-F02-F04

   CF 00  6378.171274  6378.101946  6356.751868 & -14.92911 are stored  in registers R01-R02-R03-R00
   SF 00:  you have to initialize R00-R01-R02-R03

   CF 01  CF 02  3 terms in the Fourier series:
   CF 01  SF 02   2 terms in the Fourier series:
   SF 01               1 term in the Fourier series:

   SF 04   geocentric longitudes and latitudes
   CF 04  geodetic longitudes and latitudes

Subroutines:  GX7 & GW7

-Instead if the M-Code routines GX7 & GW7, you can store the coefficients in registers R40 thru R46 ( for example )

   R40 = 0.4179591837
   R41 = 0.3818300505
   R42 = 0.4058451514
   R43 = 0.2797053915
   R44 = 0.7415311856
   R45 = 0.1294849662
   R46 = 0.9491079123

-Then, replace

   lines 158-159 by  RCL IND 24
   lines 148-149 by  RCL IND 24  DSE 24
   lines 141-142 by  46.04
   lines 136 to 138 by  RCL 40
 
 

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

 
     ( 575 bytes / SIZE 035 )
 
 

      STACK        INPUTS    OUTPUTS1    OUTPUTS2    OUTPUTS3
           T       l1 ( ° ' " )       D0 ( km )       D0 ( km )       D0 ( km )
           Z       b1 ( ° ' " )       D1 ( km )       D0 ( km )       D0 ( km )
           Y       l2 ( ° ' " )       D2 ( km )       D1 ( km )       D0 ( km )
           X       b2 ( ° ' " )       D3 ( km )       D2 ( km )       D1 ( km )

 
Example1:      Geodesic distance between the U.S. Naval Observatory at Washington (D.C.)

   l1 = 77°03'56"0 W ,  b1 = 38°55'17"2 N   and Cape Town Observatory ( South Africa )
   l2 = 18°28'35"7 E  ,  b2 = 33°56'02"5 S

    CF 01  CF 02  CF 04

  -77.0356    ENTER^
   38.55172   ENTER^
   18.28357   ENTER^
  -33.56025  XEQ "DGT3"  >>>>   D3 =  12709.47906 km = R33                        ---Execution time = 5m46s---
                                            RDN    D2 = 12709.47906 km
                                            RDN    D1 = 12709.48047 km
                                            RDN    D0 = 12709.48074 km = R21

Example2:    With ( -40° , -40° ) ( +120° , +50° ) it yields:

  •  If  CF 01  CF 02  CF 04   we get:

                     D3 =  18095.49323 km
         RDN    D2 = 18095.49323 km
         RDN    D1 = 18095.49858 km
         RDN    D0 = 18095.55363 km

  •  If  CF 01  CF 02  SF 04   we get:

                     D3 =  18099.69052 km
         RDN    D2 = 18099.69052 km
         RDN    D1 = 18099.69587 km
         RDN    D0 = 18099.75037 km

Example3:     On the ellipsoid    x2 / 41 + y2 / 37 + z2 / 35 = 1              the geocentric coordinates of the 2 points are:

    l1 =  9°17'35"55660               l2 = 63°20'07"3683
    b1 = -8°11'55"79344              b2 = 57°46'42"9935

    0   STO 00

   41  SQRT  STO 01
   37  SQRT  STO 02
   35  SQRT  STO 03

    SF 00   CF 01  CF 02  SF 04

    9.173555660  ENTER^
  -8.115579344  ENTER^
   63.20073683  ENTER^
   57.46429935  XEQ "DGT3"  >>>>   D3 = 8.594823580
                                                 RDN    D2 = 8.594824326
                                                 RDN    D1 = 8.594849588
                                                 RDN    D0 = 8.594880192

-The exact value: about  8.594822580 whence an error of  E-6 with D3

Notes:

-On the Earth, the maximum error given by D0 is of the same order as Andoyer's first order formula for an ellipsoid of revolution i-e about 60 meters.

-For the Earth, the difference between D2 & D3 is often negligible
-So, you can set F02 without a great loss of precision.

    SF 02  ->  Exectution time = 4m05s
    SF 01  ->  Execution time =  2m29s
 

2°)  Loxodromes
 

-Let  r = semi-major axis of a parallel and  rm the radius of curvature of the meridian at the point P(L,B)  L = geodetic longitude B = geodetis latitude.

     r2 = b2 + ( a2 - b2 ) / [ 1 + (b2/a2) Tan2 L ]
   rm2 = c2 / { r [ 1 - ( 1 - b2/a2 ) Sin2 B ]3/2 }

-Let  a' & b'  the semi-axes of this meridian

    a' = ( a Cos B ) / [ 1 - ( 1 - c2/a2 ) Sin2 B ]1/2
    b' = ( b Cos B ) / [ 1 - ( 1 - c2/b2 ) Sin2 B ]1/2

-Let r'm = the radius of curvature of the parallel

    r'm2 = b'2 / { a' [ 1 - ( 1 - b'2/a'2 ) Sin2 L ]3/2 }

-Then  the azimuth µ of the loxodrome verifies

        rm Tan µ dB =  r'm dL
        ( Cos µ ) ds =  rm dB

-Since r & r'm depend on L and L is an unknown function of B, we could solve a differential equation to find Tan Az
-But the eccentricity of the equator is very small, and L only appears in a corrective term,
  so we can use the expression L(B) valid for an ellipsoid of revolution

     L - L1 =  ( Tan µ ).{ Ln [ ( Tan ( 45° + B/2 ) ) / (Tan  ( 45° + B1/2 ) ) ] - (e/2).Ln [ ( 1 + e.sin B ).( 1 - e.sin B1 )/( 1 - e.sin B )/( 1 + e.sin B1 ) ] }
 

Data Registers:           •  R00 = L0               ( If SF 00 , Registers R00 thru R03 are to be initialized before executing "DGT3" )

                                      •  R01 = a                  R04 to R24: temp
                                      •  R02 = b
                                      •  R03 = c

Flags:  F00-F07-F08-F09-F10

   CF 00  6378.171274  6378.101946  6356.751868 & -14.92911 are stored  in registers R01-R02-R03-R00
   SF 00:  you have to initialize R00-R01-R02-R03

Subroutines:  GX7 & GW7

-Instead if the M-Code routines GX7 & GW7, you can store the coefficients in registers R40 thru R46 ( for example )

   R40 = 0.4179591837
   R41 = 0.3818300505
   R42 = 0.4058451514
   R43 = 0.2797053915
   R44 = 0.7415311856
   R45 = 0.1294849662
   R46 = 0.9491079123

-Then, replace

   lines 420-421 by  RCL IND 15
   lines 410-411 by  RCL IND 15  DSE 15
   lines 401 to 403 by  RCL 40
   line  397 by  46.04

-Lines 42-159-216 are three-byte GTOs
 
 

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

 
     ( 640 bytes / SIZE 025 )
 
 

      STACK        INPUTS      OUTPUTS
           T       L1 ( ° ' " )             /
           Z       b1 ( ° ' " )             /
           Y       L2 ( ° ' " )        µ ( ° ' " )
           X       b2 ( ° ' " )        D ( km )

 
Example1:

   -77.03560  ENTER^
    38.55172  ENTER^
      2.20138  ENTER^
    48.50112  XEQ "LOX3"  >>>>   D = 6453.471591 km                          ---Execution time = 2m40s---
                                            X<>Y   µ  = 80°10'15"8102

Example2:

   0   ENTER^
   0   ENTER^
  70  ENTER^
   0      R/S       >>>>  D = 7792.380608 km                          ---Execution time = 14s---
                       X<>Y   µ = 90°

Example2:

   0  ENTER^
   0  ENTER^
   0  ENTER^
  80    R/S       >>>>  D = 8885.152534 km                          ---Execution time = 15s---
                       X<>Y  µ = 0°

Notes:

-Due to roundoff errors, there is a loss of significant digits if the 2 latitudes are very close but not equal.
-If the latitudes are equal , µ = +/- 90° and if the longitudes are equal , µ = 0 or 180° so the program is much faster
-The azimuths are measured clockwise from North.