hp41programs

Geodesic2

Great Elliptic Arc Length for the HP-41


Overview
 

 1°)  Ellipsoid of Revolution
 2°)  Triaxial Ellipsoid

   a)  With Andoyer's Formula
   b)  Numerical Integration ( 2 programs )


 
 

1°)  Ellipsoid of Revolution


-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  ,  e2  = 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 = 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 )


  Dist / a  =  µ - (e2/4) ( A µ + B Sin µ ) - (e4/128) [ A2 ( 6 µ - Sin 2µ ) + 8 A B Sin µ + 2 B2 Sin 2µ ]
                         - (e6/1536) [ 3 A3 ( 10 µ - 3 Sin 2µ ) + 3 A2 B ( 15 Sin µ - Sin 3µ ) + 18 A B2 Sin 2µ  + 4 B3 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 )

 
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 "GEL"   >>>>   D = 12138.68553 km                                                      ---Execution time = 13s---

 

2°)  Triaxial Ellipsoid


       a)  With Andoyer's Formula

 

-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°929 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 ) / ( 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 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)) +(f2/4) [W/4+(SinW)(CosW)/4-(Cos2(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:     R01 = a  R02 = b  R03 = c                                                                                     R00 to R14:  temp

                                 R06 & R07 ( or R07 & R06 ) = major & minor semi-axes of the great ellipse.
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 06
 39 RCL 09
 40 *
 41 RCL 07
 42 RCL 10
 43 *
 44 +
 45 +
 46 STO 04
 47 RCL 01
 48 ST/ 05
 49 ST/ 08
 50 RCL 02
 51 ST+ 01
 52 ST/ 06
 53 ST/ 09
 54 RCL 03          
 55 ST/ 07
 56 ST/ 10
 57 RCL 05
 58 X^2
 59 RCL 06
 60 X^2
 61 RCL 07
 62 X^2
 63 +
 64 +
 65 STO 00
 66 RCL 08
 67 X^2
 68 RCL 09
 69 X^2
 70 RCL 10
 71 X^2
 72 +
 73 +
 74 X<> 10
 75 RCL 07
 76 *
 77 RCL 09
 78 RCL 06
 79 *
 80 RCL 08
 81 RCL 05
 82 *
 83 +
 84 +
 85 ST+ X
 86 STO 09          
 87 RCL 00
 88 RCL 04
 89 *
 90 ST+ X
 91 -
 92 RCL 04
 93 RCL 09
 94 *
 95 RCL 10
 96 -
 97 RCL 00
 98 RCL 04
 99 X^2
100 *
101 -
102 RCL 00
103 RCL 04
104 ACOS
105 STO 02
106 SIN
107 STO 08
108 ST/ Z
109 *
110 +
111 X#0?
112 /
113 ATAN
114 2
115 /
116 90
117 STO 06
118 X<>Y
119 X<0?
120 +
121 STO 05
122 XEQ 02
123 X<> 06
124 RCL 05          
125 +
126 XEQ 02
127 STO 07
128 RCL 02
129 D-R
130 LASTX
131 RCL 05
132 ST+ X
133 -
134 COS
135 STO 05
136 X^2
137 ST+ X
138 RCL 04
139 RCL 08
140 ST* 05
141 *
142 ST* Y
143 -
144 -
145 4
146 /
147 RCL 07          
148 RCL 06
149 /
150 1
151 -
152 2
153 ST/ 01
154 /
155 STO Z
156 *
157 RCL 05
158 -
159 R^
160 +
161 *
162 +
163 RCL 06
164 *
165 GTO 03
166 LBL 01
167 1
168 P-R
169 X<>Y
170 RCL 03
171 X^2
172 *
173 STO 10
174 X^2
175 X<> Z
176 X<>Y
177 P-R
178 RCL 01          
179 X^2
180 *
181 STO 08
182 X^2
183 X<>Y
184 RCL 02
185 X^2
186 *
187 STO 09
188 X^2
189 +
190 +
191 SQRT
192 ST/ 08
193 ST/ 09
194 ST/ 10
195 RCL 08
196 RTN
197 LBL 02
198 ENTER
199 SIN
200 STO 07
201 X^2
202 RCL 10
203 *
204 X<>Y
205 RCL 02
206 -
207 SIN
208 ST* 07
209 X^2
210 RCL 00          
211 *
212 +
213 RCL 07
214 RCL 09
215 *
216 -
217 SQRT
218 RCL 08
219 X<>Y
220 /
221 RTN
222 LBL 03
223 END

 
      ( 289 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 "GEL"   >>>>   D = 12138.68576 km                                                      ---Execution time = 17s---


       b)  Numerical Integration


-The arc length:    §0W [ R2 + ( dR/dµ )2 ] 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 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 "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 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 )

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