Great Elliptic Arc Length for the HP-41

Overview
 

 1°)  Ellipsoid of Revolution

   a)  Computing Parametric Latitudes
   b)  With Geocentric Latitudes ( without computing parametric latitudes )

 2°)  Triaxial Ellipsoid

   a)  With Andoyer's Formula of order 2
   b)  With a Formula of order 3 ( 2 programs )
   c)  Andoyer's Formula of order 3
   d)  Numerical Integration ( 2 programs )

Latest Update:    paragraph 2°)b)

 

1°)  Ellipsoid of Revolution

       a)  Computing Parametric Latitudes

-This program employs WGS84 ellipsoid:  a = 6378.137 km  &  flattening  f = 1/298.25722  ( lines 166 & 18 )

-The latitudes are geodetic.

 "GEL" employs a formula of order 3 given in reference 1 ( page 69 )
 

Formula:

  With      L , L' = longitudes  &  p , p' = parametric latitudes  ,  f = earth flattening  ,  e²  = f ( 2 - f )

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

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

              S = sin² G  cos² H  +  cos² F  sin² H                               µ = 2 Arc Sin S^1/2

 we define  A = ( sin² G  cos² F )/S + ( sin² F  cos² G )/( 1-S )
      and      B = ( sin² G  cos² F )/S - ( sin² F  cos² G )/( 1-S )

  Dist / a  =  µ - (e²/4) ( A µ + B Sin µ ) - (e⁴/128) [ A² ( 6 µ - Sin 2µ ) + 8 A B Sin µ + 2 B² Sin 2µ ]
                         - (e⁶/1536) [ 3 A³ ( 10 µ - 3 Sin 2µ ) + 3 A² B ( 15 Sin µ - Sin 3µ ) + 18 A B² Sin 2µ  + 4 B³ Sin 3µ ]

 

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to  R07:  temp
Flags: /
Subroutines: /
 
 
 

01 LBL "GEL"  02 DEG  03 HR  04 STO 01  05 X<> T  06 HMS-  07 HR  08 2  09 /  10 COS  11 X^2  12 STO 00  13 X<>Y  14 HR  15 1  16 P-R  17 LASTX  18 298.25722  19 1/X  20 STO 05            21 -  22 STO 06  23 /  24 R-P

25 X<> 01  26 1  27 P-R  28 RCL 06  29 /  30 R-P  31 RDN  32 +  33 2  34 /  35 ST- Y  36 1  37 P-R  38 X^2  39 STO 02  40 STO T  41 RDN  42 X^2  43 STO 01            44 SIGN  45 P-R  46 X^2  47 ST* 01  48 RDN

49 X^2  50 ST* 02  51 -  52 RCL 00  53 *  54 -  55 ST/ 02  56 RCL 02  57 RCL 01  58 1  59 R^  60 STO 01  61 -  62 /  63 ST- 02  64 +  65 X<> 01  66 SQRT  67 ASIN  68 ST+ X  69 STO 04            70 ST+ X  71 SIN  72 STO 03

73 3  74 *  75 RCL 04  76 D-R  77 STO 00  78 10  79 *  80 -  81 RCL 01  82 *  83 RCL 04  84 3  85 *  86 SIN  87 STO 06  88 RCL 04  89 SIN  90 STO 04            91 15  92 *  93 -  94 RCL 02  95 *  96 +

97 RCL 01  98 *  99 RCL 03 100 6 101 * 102 RCL 02 103 X^2 104 STO 07 105 * 106 - 107 3 108 * 109 RCL 01 110 * 111 RCL 06 112 4 113 * 114 RCL 02           115 * 116 RCL 07 117 * 118 - 119 2 120 RCL 05

121 ST- Y 122 * 123 4 124 / 125 STO 05 126 * 127 3 128 / 129 RCL 03 130 RCL 00 131 6 132 * 133 - 134 RCL 01 135 * 136 RCL 04 137 8 138 * 139 RCL 02           140 * 141 - 142 RCL 01 143 * 144 +

145 RCL 03 146 ST+ X 147 RCL 07 148 * 149 - 150 RCL 05 151 * 152 8 153 / 154 RCL 00 155 RCL 01 156 * 157 - 158 RCL 02 159 RCL 04 160 * 161 - 162 RCL 05           163 * 164 RCL 00 165 + 166 6378.137 167 * 168 END

 
      ( 206 bytes / SIZE 008 )
 
 

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

 b & b' = geodetic latitudes 

Example:    Calculate the great elliptic arc 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 "GEL"   >>>>   D = 12138.68553 km                                                      ---Execution time = 13s---

Notes:

-If b & b' are geocentric latitudes, simply replace lines 23 & 29 ( / ) with *

-If µ is almost 180° , we have a better accuracy if µ is computed with:

    S  = Sin² (µ/2)  = Sin² G Cos² H  + Cos² F Sin² H
    C = Cos² (µ/2) = Cos² G Cos² H + Sin² F Sin² H

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to  R07:  temp
Flags: /
Subroutines: /
 
 
 

01 LBL "GEL"  02 DEG  03 STO 01  04 X<> T  05 -  06 STO 02  07 CLX  08 SIGN  09 STO 04  10 P-R  11 LASTX  12 298.25722  13 1/X  14 STO 07  15 -  16 STO 05  17 *  18 R-P  19 X<> 01  20 RCL 04            21 P-R  22 RCL 05  23 *  24 R-P  25 RDN  26 X<>Y

27 -  28 2  29 ST/ 02  30 /  31 ST+ Y  32 RCL 02  33 X<>Y  34 RCL 04  35 P-R  36 STO 03  37 RDN  38 STO 05  39 STO 06  40 X<> L  41 P-R  42 ST* 06  43 R^  44 *  45 X^2  46 X<>Y  47 STO 00            48 R^  49 RCL 04  50 P-R  51 ST* 05  52 ST* 00

53 RDN  54 ST* 03  55 *  56 X^2  57 +  58 RCL 05  59 X^2  60 RCL 06  61 X^2  62 RCL 00  63 X^2  64 +  65 STO 01  66 /  67 STO 02  68 RCL 03            69 X^2  70 R^  71 /  72 ST- 02  73 +  74 X<> 01  75 SQRT  76 X<>Y  77 SQRT  78 R-P

79 X<>Y  80 ST+ X  81 STO 04  82 ST+ X  83 SIN  84 STO 03  85 3  86 *  87 RCL 04  88 D-R  89 STO 00  90 10  91 *  92 -  93 RCL 01  94 *  95 RCL 04            96 3  97 *  98 SIN  99 STO 06 100 RCL 04 101 SIN 102 STO 04 103 15 104 *

105 - 106 RCL 02 107 * 108 + 109 RCL 01 110 * 111 RCL 03 112 6 113 * 114 RCL 02 115 X^2 116 STO 05 117 * 118 - 119 3 120 * 121 RCL 01         122 * 123 RCL 06 124 4 125 * 126 RCL 02 127 * 128 RCL 05 129 * 130 -

131 2 132 RCL 07 133 ST- Y 134 * 135 4 136 / 137 STO 07 138 * 139 3 140 / 141 RCL 03 142 RCL 00 143 6 144 * 145 - 146 RCL 01         147 * 148 RCL 04 149 8 150 * 151 RCL 02 152 * 153 - 154 RCL 01 155 *

156 + 157 RCL 03 158 ST+ X 159 RCL 05 160 * 161 - 162 RCL 07 163 * 164 8 165 / 166 RCL 00 167 RCL 01 168 * 169 - 170 RCL 02 171 RCL 04         172 * 173 - 174 RCL 07 175 * 176 RCL 00 177 + 178 6378.137 179 * 180 END

 
      ( 220 bytes / SIZE 008 )
 
 

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

    b & b' geocentric latitudes

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

     -30   ENTER^
     -40   ENTER^
     120  ENTER^
      50   XEQ "GEL"   >>>>   D = 17433.98161 km                                                      ---Execution time = 13.3s---

Note:

-If b & b' are geodetic latitudes, simply replace lines 17 & 23 ( * ) with  /

       b)  With Geocentric Latitudes ( without computing parametric latitudes )

-If we need less accurate distances, we can use the following formula:


   D / ( a.µ ) ~  1 + (f/2) ( -A + B ( Sin µ ) / µ )
                          + (f²/4) [ -3 A + ( A²/4 ) ( 13 + ( Sin µ ) ( Cos µ ) / µ ) - ( B²/2 ) ( Sin µ ) ( Cos µ ) / µ  + 3 B ( 1 - A ) ( Sin µ ) / µ ]
                          + (f³/8) { B³ ( -7 ( Sin µ ) / µ + 10 ( Sin³ µ ) / µ ) + A² B ( 7 ( Sin µ ) / µ - 7 ( Sin³ µ ) / µ )  + ( 1 - A ) {  A ( 14 A - 9 ) + B ( ( Sin µ ) / µ ) [ - 3 B ( Cos µ ) + 4 ( 1 - A ) ] } }

Data Registers:   R00 = µ  R01 = A  R02 = B  R03 to R07:  temp
Flags: /
Subroutines: /
 
 
 

01 LBL "GEL"  02 DEG  03 X<> T  04 -  05 X<> Z  06 -  07 2  08 ST/ Z  09 /  10 ST+ Z  11 1  12 STO 04            13 P-R  14 STO 03  15 RDN  16 STO 05  17 STO 06  18 X<> L  19 P-R  20 ST* 06  21 R^  22 *  23 X^2

24 X<>Y  25 STO 00  26 R^  27 RCL 04  28 P-R  29 ST* 00  30 ST* 05  31 RDN  32 ST* 03  33 *  34 X^2  35 +  36 RCL 05            37 X^2  38 RCL 00  39 X^2  40 RCL 06  41 X^2  42 +  43 STO 01  44 /  45 STO 02  46 RCL 03

47 X^2  48 R^  49 /  50 ST- 02  51 +  52 X<> 01  53 SQRT  54 X<>Y  55 SQRT  56 R-P  57 X<>Y  58 ST+ X  59 D-R  60 STO 00  61 LASTX  62 RCL 04            63 P-R  64 STO 03  65 RDN  66 STO 07  67 ST* 07  68 /  69 STO 06

70 RCL 04  71 RCL 01  72 -  73 STO 04  74 4  75 *  76 RCL 02  77 RCL 03  78 *  79 3  80 *  81 -  82 X<>Y  83 ST/ 03  84 /  85 RCL 02            86 *  87 RCL 01  88 14  89 *  90 9  91 -  92 RCL 01

93 *  94 +  95 RCL 04  96 *  97 RCL 07            98 10  99 * 100 7 101 ST* 07 102 ST- 07 103 - 104 RCL 02 105 X^2 106 * 107 RCL 07 108 RCL 01 109 X^2 110 * 111 - 112 RCL 02 113 * 114 RCL 06 115 /

116 + 117 596.51444 118 STO 05 119 / 120 RCL 03         121 13 122 + 123 RCL 01 124 * 125 4 126 / 127 3 128 ST* 04 129 - 130 RCL 01 131 * 132 + 133 RCL 04 134 RCL 06 135 / 136 RCL 02 137 RCL 03 138 *

139 2 140 / 141 - 142 RCL 02         143 * 144 + 145 RCL 05 146 / 147 RCL 02 148 RCL 06 149 / 150 + 151 RCL 01 152 - 153 RCL 05 154 / 155 RCL 00 156 ST* Y 157 + 158 6378.137 159 * 160 END

 
      ( 204 bytes / SIZE 008 )
 
 

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

    b & b' geocentric latitudes

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

     -30   ENTER^
     -40   ENTER^
     120  ENTER^
      50   XEQ "GEL"   >>>>   D = 17433.98161 km                                                      ---Execution time = 9.9s---


-Here are a few root-mean-square & maximal errors in centimeters.

| λ'-λ |	1°	45°	90°	120°	150°	160°	170°	175°	179°
rms	0.9	0.9	0.9	0.9	1.1	1.2	1.4	1.6	1.8
max	2.4	2.4	2.4	2.4	2.4	2.4	2.7	2.9	2.9

2°)  Triaxial Ellipsoid

       a)  With Andoyer's Formula of order 2
 

-The latitudes L and the longitudes B are geodetic.

-The semi-axes are:

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

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

-The center O of the ellipsoid and 2 points M & N define a plane ( at least if they are not aligned ) that intersects the ellipsoid.
  and we can compute the distance of the arc of ellipse MN thus defined.

-A point P on the arc MN may be determined by the vectorial equality:

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

   OP = ( Sin W ) / ( A² + B² + C² )^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' ) / a² + ( y Sin W ) ( y Cos W - y' ) / b² + ( z Sin W ) ( z Cos W - z' ) / c² ]
    Q =  [ x² Sin² W - ( x Cos W - x' )² ] / a² + [ y² Sin² W - ( y Cos W - y' )² ] / b² + [ z² Sin² W - ( z Cos W - z' )² ] / c²

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

    A² + B² + C² = T Sin² ( W - µ ) + S Sin² µ + 2 U Sin µ Sin ( W - µ )

   where  S = x'² / a² + y'² / b² + z'² / c²  ,   T =  x² / a² + y² / b² + z² / c²  ,   U =  xx' / a² + yy' / b² + zz' / c²

    and   P = 2 U - 2 T Cos W  ,  Q = T Sin W + ( 2 U Cos W - S - T Cos² W ) / Sin W

-Then, the distance on the great ellipse of the triaxial ellipsoid - between M & N is computed by Andoyer's formula of order 2, re-written with geocentric "latitudes":

   D/a' = W + (f/2) (-W+(Cos2F)(SinW)) +(f²/4) [W/4+(SinW)(CosW)/4-(Cos²(2F))(SinW)(CosW)/2]

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

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

 
 
 

01 LBL "GEL"  02 DEG  03 HR  04 STO 14  05 RDN  06 HR  07 14.92911  08 STO 11  09 +  10 STO 13  11 RDN  12 HR  13 STO 12  14 X<>Y  15 HR  16 ST+ 11  17 6378.172  18 STO 01  19 .07  20 -  21 STO 02            22 21.349686  23 -  24 STO 03  25 RCL 11  26 RCL 12  27 XEQ 01  28 STO 05

29 RCL 09  30 STO 06  31 RCL 10  32 STO 07  33 RCL 13  34 RCL 14  35 XEQ 01  36 XEQ 00  37 STO 04  38 RCL 01  39 ST/ 05  40 ST/ 08  41 RCL 02  42 ST/ 06  43 ST/ 09  44 RCL 03  45 ST/ 07  46 ST/ 10  47 RCL 08  48 X^2  49 RCL 09            50 X^2  51 RCL 10  52 X^2  53 +  54 +  55 STO 13  56 RCL 05

57 X^2  58 RCL 06  59 X^2  60 RCL 07  61 X^2  62 +  63 +  64 STO 00  65 -  66 RCL 04  67 ACOS  68 STO 02  69 SIN  70 STO 15  71 *  72 RCL 08  73 XEQ 00  74 ST+ X  75 STO 09  76 RCL 00            77 RCL 13  78 +  79 RCL 04  80 *  81 -  82 R-P  83 X<>Y  84 STO 08

85 2  86 /  87 STO 05  88 XEQ 02  89 STO 06  90 RCL 05  91 90  92 +  93 XEQ 02  94 STO 07  95 RCL 02  96 D-R  97 RCL 08            98 COS  99 STO 05 100 X^2 101 ST+ X 102 RCL 04 103 RCL 15 104 ST* 05 105 * 106 ST* Y 107 - 108 - 109 4 110 / 111 RCL 06 112 ST/ 15

113 RCL 07 114 ST- Y 115 ST+ X 116 / 117 STO 07 118 * 119 RCL 05 120 - 121 + 122 RCL 07 123 * 124 + 125 RCL 15         126 * 127 GTO 03 128 LBL 00 129 RCL 05 130 * 131 RCL 09 132 RCL 06 133 * 134 + 135 RCL 10 136 RCL 07 137 * 138 + 139 RTN 140 LBL 01

141 1 142 P-R 143 X<>Y 144 RCL 03 145 X^2 146 * 147 STO 10 148 RDN 149 P-R 150 RCL 01 151 X^2 152 * 153 STO Z 154 X^2 155 X<>Y 156 RCL 02         157 X^2 158 * 159 STO 09 160 X^2 161 RCL 10 162 X^2 163 + 164 + 165 SQRT 166 ST/ 09 167 ST/ 10 168 / 169 STO 08

170 RTN 171 LBL 02 172 RCL 02 173 2 174 / 175 + 176 ENTER 177 SIN 178 ENTER 179 X^2 180 RCL 13 181 * 182 RCL 02         183 R^ 184 - 185 SIN 186 ST* Z 187 X^2 188 RCL 00 189 * 190 + 191 X<>Y 192 RCL 09 193 * 194 + 195 SQRT 196 RTN 197 LBL 03 198 END

 
      ( 262 bytes / SIZE 016 )
 
 

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

 
Example:    Calculate the great elliptic arc 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 "GEL"   >>>>   D = 12138.65875 km                                                      ---Execution time = 16s---


       b)  With a Formula of Order 3


-We can calculate D more directly:


               m = [ y² (a²/b² - 1)/2 + z² (a²/c² - 1)/2 ] / ( Sin² W )
    With    n  = [ y'² (a²/b² - 1)/2 + z'² (a²/c² - 1)/2 ] / ( Sin² W )      
               p  = 2 [ y.y' (a²/b² - 1)/2 + z.z' (a²/c² - 1)/2 ] / ( Sin² W )

 R/a = 1 / SQRT [ 1 + m Sin² (W-µ) + n Sin² µ + p Sin µ Sin (W-µ) ]  = 1 / SQRT ( 1 + q )     where q is a small number

-After some calculations, it yields:  

     1 + q = ( 1 + m Sin² W ) [ 1 + e Sin² µ + g Sin µ Cos µ ]  =  ( 1 + m Sin² W ) ( 1 + e/2 )  [ 1 + E Cos 2.µ + F Sin 2.µ ] =  ( 1 + m Sin² W ) ( 1 + e/2 ) [ 1 + G Sin ( 2.µ + D ) ]

              = ( 1 + m Sin² W ) ( 1 + e/2 ) ( 1 + h )


-The formula of order 3 to compute the arc length  =  §₀^W [ R² + ( dR/dµ )² ] ^1/2 dµ   is

    [ R² + ( dR/dµ )² ] ^1/2  ~   a  ( 1 + m Sin² W ) ^-1/2 (1 + e/2 )^-1/2 [ 1 - h/2 + (3/8) h² - (5/16) h³ + (1/8) h'² - (5/16) h h'² ]             where  h' = dh/dµ


      D / a   =   ( 1 + m Sin² W )^-1/2 ( 1 + e/2 )^-1/2 {  W + ( G / 4 )    [ - Cos D + Cos ( 2.µ + D ) ]  + ( G² / 32 ) [ 14 W - ( Sin D ) ( Cos D ) + Sin ( 2.W + D ) Cos ( 2.W + D ) ]
                                                                                   + ( G³ / 32 ) [ -5 Cos D + 5 Cos ( 2.W + D ) - 5 Cos³ D + 5 Cos³ ( 2.W + D ) ] }



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

 
   

01 LBL "GEL"  02 DEG  03 HR  04 STO 14  05 RDN  06 HR  07 14.92911  08 STO 11  09 +  10 STO 13  11 RDN  12 HR  13 STO 12  14 X<>Y  15 HR  16 ST+ 11  17 6378.172  18 STO 01  19 .07  20 -  21 STO 02  22 21.349686  23 -  24 STO 03  25 RCL 11  26 RCL 12            27 XEQ 01  28 STO 05  29 RCL 09  30 STO 06  31 RCL 10  32 STO 07

33 RCL 13  34 RCL 14  35 XEQ 01  36 RCL 05  37 *  38 RCL 09  39 RCL 06  40 *  41 STO 05  42 +  43 RCL 10  44 RCL 07  45 *  46 STO 08  47 +  48 STO 00  49 RCL 09  50 X^2  51 RCL 01  52 RCL 02  53 /  54 X^2  55 1  56 -  57 ST* 05  58 STO 09            59 *  60 RCL 01  61 RCL 03  62 /  63 X^2  64 1

65 -  66 STO 04  67 ST* 08  68 RCL 10  69 X^2  70 *  71 +  72 STO 10  73 RCL 06  74 X^2  75 RCL 09  76 *  77 RCL 07  78 X^2  79 RCL 04  80 *  81 +  82 STO 06  83 RCL 05  84 RCL 08  85 +  86 ST+ X  87 RCL 00  88 ACOS  89 STO 08            90 SIN  91 STO 07  92 X^2  93 ST/ 06  94 ST/ 10  95 /  96 STO 09

97 RCL 00  98 *  99 RCL 10 100 - 101 RCL 08 102 ST+ X 103 STO 05 104 COS 105 RCL 06 106 * 107 - 108 RCL 09 109 RCL 00 110 ST+ X 111 RCL 06 112 * 113 - 114 RCL 06 115 RCL 07 116 ST* Z 117 X^2 118 * 119 1 120 + 121 STO 03         122 ST+ X 123 ST/ Z 124 / 125 1 126 RCL Z 127 - 128 ST* 03

129 ST/ Z 130 / 131 R-P 132 STO 10 133 X<>Y 134 STO 02 135 1 136 P-R 137 STO 00 138 X<>Y 139 STO 06 140 RCL 02 141 RCL 05 142 + 143 1 144 P-R 145 STO 07 146 ST* Y 147 X^2 148 LASTX 149 ST* Y 150 + 151 RCL 00         152 ST* 06 153 ST- 07 154 X^2 155 LASTX 156 ST* Y 157 + 158 - 159 RCL 10

160 STO T 161 * 162 5 163 * 164 + 165 RCL 06 166 - 167 RCL 08 168 D-R 169 STO 00 170 14 171 * 172 + 173 * 174 8 175 / 176 RCL 07 177 + 178 * 179 4 180 / 181 RCL 00         182 + 183 GTO 02 184 LBL 01 185 1 186 P-R 187 X<>Y 188 RCL 03 189 X^2 190 *

191 STO 10 192 RDN 193 P-R 194 RCL 01 195 X^2 196 * 197 STO Z 198 X^2 199 X<>Y 200 RCL 02 201 X^2 202 * 203 STO 09 204 X^2 205 RCL 10 206 X^2 207 + 208 + 209 SQRT 210 ST/ 09 211 ST/ 10 212 / 213 STO 08         214 RTN 215 LBL 02 216 RCL 03 217 SQRT 218 / 219 RCL 01 220 * 221 END

 
      ( 282 bytes / SIZE 016 )
 
 

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

 
Example:    Calculate the great elliptic arc 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 "GEL"   >>>>   D = 12138.65875 km                                                      ---Execution time = 15s---


Note:

-We can simplify a little this program if the coordinates are geocentric.
-In the following version, the geocentric longitudes are measured from the equatorial major-axis


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

 
   

01 LBL "GEL"  02 DEG  03 HR  04 STO 02  05 RDN  06 HR  07 STO 06  08 X<> Z  09 HR  10 X<>Y  11 HR  12 1  13 STO 04            14 P-R  15 X<>Y  16 STO 07  17 RDN  18 P-R  19 STO 05  20 X<>Y  21 X<> 06  22 RCL 02  23 RCL 04  24 P-R  25 X<>Y

26 STO 01  27 RDN  28 P-R  29 RCL 05  30 *  31 X<>Y  32 STO 02  33 RCL 06  34 *  35 STO 05  36 +  37 RCL 01            38 RCL 07  39 *  40 STO 08  41 +  42 STO 00  43 RCL 02  44 X^2  45 2195 E-8  46 ST* 05  47 STO 02  48 *  49 6750546 E-9  50 STO 03

51 ST* 08  52 RCL 01  53 X^2  54 *  55 +  56 STO 01  57 RCL 06  58 X^2  59 RCL 02  60 *  61 RCL 07            62 X^2  63 RCL 03  64 *  65 +  66 STO 06  67 RCL 05  68 RCL 08  69 +  70 ST+ X  71 RCL 00  72 ACOS  73 STO 08  74 SIN  75 STO 07

76 X^2  77 ST/ 01  78 ST/ 06  79 /  80 STO 02  81 RCL 00  82 *  83 RCL 01            84 -  85 RCL 08  86 ST+ X  87 STO 05  88 COS  89 RCL 06  90 *  91 -  92 RCL 02  93 RCL 00  94 ST+ X  95 RCL 06  96 *  97 -  98 RCL 06  99 RCL 07 100 ST* Z

101 X^2 102 * 103 RCL 04 104 + 105 STO 03 106 ST+ X 107 ST/ Z 108 / 109 RCL 04         110 RCL Z 111 - 112 ST* 03 113 ST/ Z 114 / 115 R-P 116 STO 01 117 X<>Y 118 STO 02 119 RCL 04 120 P-R 121 STO 00 122 X<>Y 123 STO 06 124 RCL 02

125 RCL 05 126 + 127 RCL 04 128 P-R 129 STO 07 130 ST* Y 131 X^2 132 LASTX 133 ST* Y 134 + 135 RCL 00 136 ST* 06 137 ST- 07 138 X^2 139 LASTX 140 ST* Y 141 + 142 - 143 RCL 01         144 STO T 145 * 146 5 147 * 148 +

149 RCL 06 150 - 151 RCL 08 152 D-R 153 STO 00 154 14 155 * 156 + 157 * 158 8 159 / 160 RCL 07         161 + 162 * 163 4 164 / 165 RCL 00 166 + 167 RCL 03 168 SQRT 169 / 170 6378.172 171 * 172 END

 
      ( 224 bytes / SIZE 009 )
 
 

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

 
Example:   

 L = 166°08'02"2      L' = -101°56'06"0
 b = -33°41'00"9       b' =  33°10'47"0

     166.08022  ENTER^
    -33.41009   ENTER^
   -101.56060  ENTER^
     33.10470   XEQ "GEL"   >>>>   D = 12138.70371 km                                                      ---Execution time = 12s---


Notes:

-With free42 , replace line 45 by  2195022 E-11  and line 49 by  675054581 E-11
-The maximum error ( except with ~ antipodal points ) is only  1 mm (!)


       c)  Andoyer's Formula of order 3


-Though it seems impossible to find a formula of order 3 without calculating parametric latitudes in general case, I've found such a formula if the longitudes are the same.
-So, to compute the great elliptic arc distance with the geocentric coordinates, here is the formula employed by the following program ( lines 134... ):


  D/a = W + (f/2) ( -W + B.SinW )  + (f ² /4) [ W/4 + ( 1-2B² ) ( SinW ) ( CosW )/4 ] + (f ³/8) [ W/4 + ( 1 - 2B² + 30.B.CosW ) ( SinW ) ( CosW )/4 - (5/2) B³ Sin (3W) ]

     a = major axis of the great ellipse , f = eccentricity
   2F = b1+b2 where b1 & b2 = geocentric "latitudes"   B = Cos(2F)
    W = | b1 - b2 |


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

   

01 LBL "GEL"  02 DEG  03 HR  04 STO 14  05 RDN  06 HR  07 STO 13  08 RDN  09 HR  10 STO 12  11 6378.172  12 STO 01  13 .07  14 -  15 STO 02  16 21.349686  17 -  18 STO 03  19 R^  20 HR  21 14.92911  22 ST+ 13  23 +  24 STO 11  25 RCL 12  26 XEQ 01  27 STO 07            28 RCL 09  29 STO 06  30 RCL 08  31 STO 05  32 RCL 13

33 RCL 14  34 XEQ 01  35 RCL 07  36 -  37 X^2  38 RCL 09  39 RCL 06  40 -  41 X^2  42 RCL 08  43 RCL 05  44 -  45 X^2  46 +  47 +  48 SQRT  49 2  50 /  51 STO 00  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 08  62 X^2  63 RCL 09  64 X^2

65 RCL 10  66 X^2  67 +  68 +  69 STO 13  70 RCL 05  71 ST* 08  72 X^2  73 RCL 06  74 ST* 09  75 X^2  76 RCL 07  77 ST* 10  78 X^2  79 +  80 +  81 ST+ 13  82 -  83 RCL 00  84 ASIN  85 ST+ X  86 STO 02  87 1  88 P-R  89 STO 04            90 RDN  91 STO 11  92 *  93 RCL 08  94 RCL 09  95 RCL 10  96 +

97 +  98 ST+ X  99 STO 09 100 RCL 13 101 ST+ 09 102 RCL 04 103 * 104 - 105 R-P 106 X<>Y 107 STO 08 108 2 109 / 110 X<>Y 111 SQRT 112 P-R 113 X^2 114 CHS 115 X<>Y 116 X^2 117 RCL 00 118 X^2 119 RCL 09         120 * 121 ST+ Z 122 + 123 STO 05 124 SQRT 125 STO 06 126 X<>Y 127 X<=0? 128 X<> 05

129 SQRT 130 ST- Y 131 ST+ X 132 / 133 STO 05 134 3 135 RCL 11 136 ST+ X 137 X^2 138 - 139 10 140 * 141 RCL 08 142 COS 143 STO 09 144 X^2 145 ST* Y 146 ST+ X 147 1 148 - 149 STO 13         150 RCL 09 151 ST* Z 152 RCL 04 153 * 154 30 155 * 156 - 157 RCL 04 158 * 159 + 160 RCL 11

161 ST* 09 162 ST* 13 163 * 164 RCL 02 165 D-R 166 STO 00 167 STO Z 168 - 169 RCL 05 170 * 171 RCL 04 172 RCL 13 173 * 174 - 175 + 176 RCL 05 177 * 178 4 179 / 180 RCL 09         181 - 182 + 183 RCL 05 184 * 185 + 186 RCL 06 187 / 188 GTO 02 189 LBL 01 190 1 191 P-R 192 X<>Y

193 RCL 03 194 X^2 195 * 196 STO 10 197 RDN 198 P-R 199 RCL 01 200 X^2 201 * 202 STO 08 203 X^2 204 X<>Y 205 RCL 02 206 X^2 207 * 208 STO 09 209 X^2 210 RCL 10 211 X^2 212 + 213 + 214 SQRT 215 ST/ 08 216 ST/ 09 217 ST/ 10 218 RCL 10         219 RTN 220 LBL 02 221 RCL 11 222 * 223 END

 
      ( 290 bytes / SIZE 015 )
 
 

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

 
Example:    Calculate the great elliptic arc 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 "GEL"   >>>>   D = 12138.65875 km                                                      ---Execution time = 15s---


       d)  Numerical Integration

-The arc length:    §₀^W [ R² + ( dR/dµ )² ] ^1/2 dµ    is calculated by Newton-Cotes formula of order 9

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

-Line 48 is a three-byte GTO 03
 
 

01 LBL "GEL"  02 DEG  03 HR  04 STO 04  05 RDN  06 HR  07 14.929  08 STO 01  09 +  10 STO 03  11 RDN  12 HR  13 STO 02  14 X<>Y  15 HR  16 ST+ 01  17 6378.172  18 STO 05  19 .07  20 -  21 STO 06  22 21.349686  23 -  24 STO 07  25 RCL 01  26 RCL 02  27 XEQ 02  28 STO 08            29 RCL 12  30 STO 09  31 RCL 13  32 STO 10  33 RCL 03  34 RCL 04  35 XEQ 02

36 RCL 08  37 *  38 RCL 09  39 RCL 12  40 *  41 +  42 RCL 10  43 RCL 13  44 *  45 +  46 ACOS  47 STO 14  48 GTO 03  49 LBL 01  50 STO 01  51 1  52 P-R  53 STO 03  54 X<>Y  55 STO 04  56 RCL 14  57 RCL 01  58 -  59 1  60 P-R  61 STO 01  62 X<>Y  63 STO 02  64 RCL 08            65 *  66 RCL 04  67 RCL 11  68 *  69 +  70 RCL 05

71 /  72 STO 16  73 X^2  74 RCL 02  75 RCL 09  76 *  77 RCL 04  78 RCL 12  79 *  80 +  81 RCL 06  82 /  83 STO 17  84 X^2  85 +  86 RCL 02  87 RCL 10  88 *  89 RCL 04  90 RCL 13  91 *  92 +  93 RCL 07  94 /  95 STO 18  96 X^2  97 +  98 STO 15            99 RCL 03 100 RCL 11 101 * 102 RCL 01 103 RCL 08 104 * 105 -

106 RCL 16 107 * 108 RCL 05 109 / 110 RCL 03 111 RCL 12 112 * 113 RCL 01 114 RCL 09 115 * 116 - 117 RCL 17 118 * 119 RCL 06 120 / 121 + 122 RCL 03 123 RCL 13 124 * 125 RCL 01 126 RCL 10 127 * 128 - 129 RCL 18 130 * 131 RCL 07           132 / 133 + 134 RCL 15 135 / 136 X^2 137 1 138 + 139 RCL 15 140 /

141 SQRT 142 RTN 143 LBL 02 144 1 145 P-R 146 X<>Y 147 RCL 07 148 X^2 149 * 150 STO 13 151 X^2 152 X<> Z 153 X<>Y 154 P-R 155 RCL 05 156 X^2 157 * 158 STO 11 159 X^2 160 X<>Y 161 RCL 06 162 X^2 163 * 164 STO 12           165 X^2 166 + 167 + 168 SQRT 169 ST/ 11 170 ST/ 12 171 ST/ 13 172 RCL 11 173 RTN 174 LBL 03 175 CLX

176 XEQ 01 177 STO 19 178 RCL 14 179 XEQ 01 180 ST+ 19 181 RCL 14 182 9 183 / 184 STO 00 185 XEQ 01 186 STO 20 187 RCL 00 188 8 189 * 190 XEQ 01 191 ST+ 20 192 RCL 00 193 ST+ X 194 XEQ 01 195 STO 21 196 RCL 00           197 7 198 * 199 XEQ 01 200 ST+ 21 201 RCL 00 202 3 203 * 204 XEQ 01 205 STO 22 206 RCL 00 207 6 208 * 209 XEQ 01 210 ST+ 22

211 RCL 00 212 4 213 * 214 XEQ 01 215 X<> 00 216 5 217 * 218 XEQ 01 219 RCL 00 220 + 221 5778 222 * 223 RCL 22 224 19344 225 * 226 + 227 RCL 21 228 1080 229 * 230 + 231 RCL 20           232 15741 233 * 234 + 235 RCL 19 236 2857 237 * 238 + 239 RCL 14 240 ST* Y 241 SIN 242 * 243 D-R 244 89600 245 / 246 END

 
      ( 349 bytes / SIZE 023 )
 
 

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

 
Example:    Calculate the great elliptic arc 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 "GEL"   >>>>   D = 12138.65876 km                                                      ---Execution time = 56s---


Notes:

  W is computed from  Cos W = x x' + y y' + z z'  so the precision may be small if W is very small

-With the formula:   W = 2 Arc Sin ( || u - v || / 2 )  the accuracy is better with small W-values.
-But the precision will be even better if we use

   || u x v || = Sin W  &  u.v  = Cos W

-Gauss-Legendre 10-point formula is employed.
-The interval  [0,W]  is divided in n subintervals


Data Registers:         •  R00 = n                                          ( Register R00 is to be initialized before executing "GEL" )

                                       R05 = a  R06 = b  R07 = c              R01 to R34:  temp
Flags: /
Subroutines: /
 

-Line 76 is a three-byte GTO 03

 

01 LBL "GEL"  02 DEG  03 HR  04 STO 04  05 RDN  06 HR  07 14.929  08 STO 01  09 +  10 STO 03  11 RDN  12 HR  13 STO 02  14 X<>Y  15 HR  16 ST+ 01  17 6378.172  18 STO 05  19 .07  20 -  21 STO 06  22 21.349686  23 -  24 STO 07  25 RCL 01  26 RCL 02  27 XEQ 02  28 STO 08  29 RCL 12  30 STO 09                  31 RCL 13  32 STO 10            33 RCL 03  34 RCL 04  35 XEQ 02  36 RCL 10  37 *  38 RCL 08

39 RCL 13  40 *  41 -  42 X^2  43 RCL 09  44 RCL 13  45 *  46 RCL 10  47 RCL 12  48 *  49 -  50 X^2  51 +  52 RCL 08  53 RCL 12  54 *  55 RCL 09  56 RCL 11  57 *  58 -  59 X^2  60 +  61 SQRT  62 RCL 08  63 RCL 11  64 *  65 RCL 09  66 RCL 12                  67 *  68 +  69 RCL 10            70 RCL 13  71 *  72 +  73 R-P  74 X<>Y  75 STO 14  76 GTO 03

77 LBL 01  78 STO 01  79 1  80 P-R  81 STO 03  82 X<>Y  83 STO 04  84 RCL 14  85 RCL 01  86 -  87 1  88 P-R  89 STO 01  90 X<>Y  91 STO 02  92 RCL 08  93 *  94 RCL 04  95 RCL 11  96 *  97 +  98 RCL 05  99 / 100 STO 16 101 X^2 102 RCL 02 103 RCL 09 104 * 105 RCL 04                 106 RCL 12 107 * 108 + 109 RCL 06 110 / 111 STO 17 112 X^2 113 + 114 RCL 02

115 RCL 10 116 * 117 RCL 04 118 RCL 13 119 * 120 + 121 RCL 07 122 / 123 STO 18 124 X^2 125 + 126 STO 15 127 RCL 03 128 RCL 11 129 * 130 RCL 01 131 RCL 08 132 * 133 - 134 RCL 16 135 * 136 RCL 05 137 / 138 RCL 03 139 RCL 12 140 * 141 RCL 01 142 RCL 09                 143 * 144 - 145 RCL 17 146 * 147 RCL 06 148 / 149 + 150 RCL 03 151 RCL 13 152 *

153 RCL 01 154 RCL 10 155 * 156 - 157 RCL 18 158 * 159 RCL 07 160 / 161 + 162 RCL 15 163 / 164 X^2 165 1 166 + 167 RCL 15 168 / 169 SQRT 170 RTN 171 LBL 02 172 1 173 P-R 174 X<>Y 175 RCL 07 176 X^2 177 * 178 STO 13 179 X^2 180 X<> Z 181 X<>Y 182 P-R 183 RCL 05                 184 X^2 185 * 186 STO 11 187 X^2 188 X<>Y 189 RCL 06 190 X^2

191 * 192 STO 12 193 X^2 194 + 195 + 196 SQRT 197 ST/ 11 198 ST/ 12 199 ST/ 13 200 RCL 11 201 RTN 202 LBL 03 203 .2955242247 204 STO 21 205 .148874339 206 STO 22 207 .2692667193 208 STO 23 209 .4333953941 210 STO 24 211 .2190863625 212 STO 25 213 .6794095683 214 STO 26 215 .1494513492 216 STO 27 217 .8650633667 218 STO 28 219 6.667134431 E-2 220 STO 29                 221 .9739065285 222 STO 30 223 CLX 224 STO 19 225 STO 20 226 RCL 14 227 RCL 00 228 STO 31

229 ST+ X 230 / 231 STO 32 232 LBL 04 233 ST+ 20 234 30.02 235 STO 33 236 LBL 05 237 RCL 20 238 STO 34 239 RCL 32 240 RCL IND 33 241 * 242 ST+ 34 243 - 244 XEQ 01 245 X<> 34 246 XEQ 01 247 RCL 34 248 + 249 DSE 33 250 RCL IND 33 251 * 252 ST+ 19 253 DSE 33 254 GTO 05 255 RCL 32 256 ST+ 20 257 DSE 31 258 GTO 04 259 RCL 19                 260 * 261 RCL 14 262 SIN 263 * 264 D-R 265 END

 
      ( 436 bytes / SIZE 035 )
 
 

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

 
Example:    Calculate the great elliptic arc 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 )

 •   With  n = 1  STO 00

     151.12178  ENTER^
    -33.51411   ENTER^
   -116.51504  ENTER^
     33.21224   XEQ "GEL"   >>>>   D = 12138.65876 km                                                      ---Execution time = 63s---

-Same result with  n = 2


Notes:

-For the Earth, n = 1 is enough for a good precision.
-If the ellipsoid has larger flattenings, test this program with different n-values to obtain the best accuracy.
 ( modify lines 07 & 17 to 24 to store the different axis-values )

 

Reference:

[1]  Paul D. Thomas - "MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS" - US Naval Oceanographic Office.