Ellipsoid-Plane

Ellipsoid-Plane Intersection for the HP-41

Overview

1°) General Case
2°) Great Ellipse defined by 2 points on the Earth.

-The intersection of an ellipsoid and a plane is an ellipse ( or a point or the empty set )
-"EPI" calculates the elements of this ellipse, assuming the plane & the ellipsoid are given by

(P) : m x + n y + p z = q
(E) : x² / a² + y² / b² + z² / c² = 1

-The second program computes the elements of the great ellipse defined by 2 points on the Earth.

1°) General Case

"EPI" employs the method described in reference [1]

Data Registers: • R00 = q ( Registers R00 thru R06 are to be initialized before executing "EPI" )

                                      • R01 = m     • R04 = a               R07-R08-R09 = f₀           R16-R17-R18 = Vertex1         R25-R26-R27 = Vertex4
                                      • R02 = n      • R05 = b               R10-R11-R12 = f₁           R19-R20-R21 = Vertex2         R28 = t₀
                                      • R03 = p      • R06 = c               R13-R14-R15 = f₂           R22-R23-R24 = Vertex3         R29 & R30 = semi-axes
Flags: /
Subroutines: /

01 LBL "EPI"
02 DEG
03 CLX
04 STO 12
05 RCL 00
06 RCL 01
07 RCL 04
08 *
09 STO 07
10 STO 13
11 CHS
12 STO 11
13 X^2
14 RCL 02
15 RCL 05
16 *
17 STO 08
18 STO 10
19 STO 14
20 X^2
21 +
22 RCL 03
23 RCL 06
24 *
25 STO 09
26 STO 15
27 X^2
28 +
29 SQRT
30 ST/ 07
31 ST/ 08

32 ST/ 09
33 ST/ 10
34 ST/ 11
35 ST/ 13
36 ST/ 14
37 ST/ 15
38 /
39 ST* 07
40 ST* 08
41 ST* 09
42 ASIN
43 COS
44 ST* 10
45 ST* 11
46 RCL 13
47 RCL 14
48 R-P
49 X#0?
50 GTO 00
51 STO 11
52 STO 13
53 STO 15
54 X<> Z
55 STO 10
56 STO 14
57 GTO 01
58 LBL 00
59 ST/ 10
60 ST/ 11
61 RCL 11
62 RCL 15

63 *
64 CHS
65 X<> 13
66 RCL 10
67 LASTX
68 *
69 X<> 14
70 RCL 10
71 *
72 X<>Y
73 RCL 11
74 *
75 X<>Y
76 -
77 STO 15
78 LBL 01
79 RCL 04
80 ST* 07
81 ST* 10
82 ST* 13
83 RCL 05
84 ST* 08
85 ST* 11
86 ST* 14
87 RCL 06
88 ST* 09
89 ST* 15
90 RCL 10
91 RCL 13
92 *
93 RCL 11

94 RCL 14
95 *
96 +
97 ST+ X
98 RCL 10
99 X^2
100 RCL 13
101 X^2
102 -
103 RCL 11
104 X^2
105 +
106 RCL 14
107 X^2
108 -
109 RCL 15
110 X^2
111 -
112 R-P
113 CLX
114 2
115 /
116 STO 28
117 XEQ 10
118 STO 16
119 RDN
120 STO 17
121 X<>Y
122 STO 18
123 RCL 28
124 90

125 +
126 XEQ 10
127 STO 19
128 RDN
129 STO 20
130 X<>Y
131 STO 21
132 RCL 28
133 90
134 -
135 XEQ 10
136 STO 22
137 RDN
138 STO 23
139 X<>Y
140 STO 24
141 RCL 28
142 180
143 +
144 XEQ 10
145 GTO 14
146 LBL 10
147 1
148 P-R
149 RCL 10
150 X<>Y
151 *
152 STO M
153 CLX
154 RCL 11
155 LASTX

156 *
157 STO N
158 RCL 13
159 RCL Z
160 *
161 RCL 07
162 +
163 ST+ M
164 CLX
165 RCL 14
166 R^
167 *
168 RCL 08
169 +
170 ST+ N
171 RCL 15
172 R^
173 *
174 RCL 09
175 +
176 RCL N
177 RCL M
178 CLA
179 RTN
180 LBL 14
181 STO 25
182 RCL 07
183 -
184 X^2
185 X<>Y
186 STO 26

187 RCL 08
188 -
189 X^2
190 +
191 X<>Y
192 STO 27
193 RCL 09
194 -
195 X^2
196 +
197 SQRT
198 RCL 07
199 RCL 19
200 -
201 X^2
202 RCL 08
203 RCL 20
204 -
205 X^2
206 +
207 RCL 09
208 RCL 21
209 -
210 X^2
211 +
212 SQRT
213 X>Y?
214 X<>Y
215 STO 30
216 X<>Y
217 STO 29
218 END

( 296 bytes / SIZE 031 )

STACK	INPUTS	OUTPUTS
Y	/	B
X	/	A

Where A & B are the semi-axes of the ellipse.

Example: (P) : 2 x + 3 y + 4 z = 5 (E) : x² / 25 + y² / 16 + z² / 9 = 1

5 STO 00

      2 STO 01       5 STO 04
      3 STO 02       4 STO 05
      4 STO 03       3 STO 06

XEQ "EPI" >>>> A = 4.547809219 = R29 ---Execution time = 15s---
X<>Y B = 3.374481476 = R30

And we find in registers R07 to R28 ( rounded to 4 decimals )

     R07-R08-R09 = ( 0.6443 , 0.6186 , 0.4639 )= f₀ = center of the ellipse
     R10-R11-R12 = ( 3.7153 , -2.4769 , 0 )        = f₁
     R13-R14-R15 = ( 1.8862 , 1.8107 , -2.3011 ) = f₂

-So any point of the ellipse may be obtained by the formula: f₀ + f₁ Cos t + f₂ Sin t where t belongs to [0°,360°[

And the 4 vertices are in registers R16 to R27:

     R16-R17-R18 = ( 4.7415 , -1.2447 , -0.1872 ) = V₁
     R19-R20-R21 = ( 1.4021 , 3.0561 , -1.7432 ) = V₂
     R22-R23-R24 = ( -0.1135 , -1.8190 , 2.6710 ) = V₃
     R25-R26-R27 = ( -3.4528 , 2.4818 , 1.1150 ) = V₄

-In register R28 we have the value of t which corresponds to the 1st vertex, i-e t = 16°43698070

Notes:

-You can use LBL 10 to compute the coordinates of a point, for a given t-value.
-For example, with t = 41° ( t in decimal degree )

    41 XEQ 10 returns x = +4.6857
                          RDN   y = -0.0628
                          RDN   z = -1.0457

-There will be a data error message if the intersection is empty.

2°) Great Ellipse defined by 2 points on the Earth

-Given 2 points M & N at the surface of the Earth by their geodetic coordinates, "GEL" returns the semi-major and semi-minor axes of the great ellipse
passing by these 2 points and the center of the Earth.

-The coordinates of M & N on the great ellipse ( the new axes = the axes of the ellipse ) are also computed

Data Registers: R00 = t₀ R14 = 6378.172 R15 = 6378.102 R16 = 6356.752314

-> R01 = x_M R03 = x_N R05-R06-R07 = X_M Y_M Z_M R11-R12-R13: temp R23-R24-R25 = U = unit vector of the 1st vertex
-> R02 = y_M R04 = y_N R08-R09-R10 = X_N Y_N Z_N R17 to R22: temp R26-R27-R28 = V = unit vector of the 2nd vertex

R29 = A = semi-major axis R30 = B = semi-minor axis

Flags: /
Subroutines: /

01 LBL "GEL"
02 DEG
03 HR
04 STO 04
05 RDN
06 HR
07 14.929
08 STO 00
09 +
10 STO 03
11 RDN
12 HR
13 X<>Y
14 HR
15 RCL 00
16 +
17 XEQ 00
18 STO 05
19 RDN
20 STO 06
21 X<>Y
22 STO 07
23 RCL 04
24 RCL 03
25 XEQ 00
26 STO 08
27 RDN
28 STO 09
29 X<>Y
30 STO 10
31 GTO 01
32 LBL 00
33 STO 11
34 SIN
35 X<>Y
36 STO 12
37 COS
38 *
39 X^2
40 6356.752314
41 STO 16
42 X^2
43 X<>Y
44 6378.102
45 STO 15
46 X^2
47 6378.172
48 STO 14
49 X^2
50 ST- Y

51 ST- T
52 ST/ T
53 /
54 STO 13
55 *
56 X<>Y
57 RCL 12
58 SIN
59 X^2
60 *
61 +
62 1
63 ST+ 13
64 ST+ Z
65 +
66 SQRT
67 RCL 14
68 X<>Y
69 /
70 RCL 11
71 RCL 12
72 X<> Z
73 X<>Y
74 RDN
75 P-R
76 R^
77 X<>Y
78 P-R
79 R^
80 ST* T
81 X<> 13
82 ST* Z
83 RDN
84 RTN
85 LBL 01
86 RCL 06
87 RCL 10
88 *
89 RCL 07
90 RCL 09
91 *
92 -
93 RCL 14
94 *
95 STO 11
96 RCL 07
97 RCL 08
98 *
99 RCL 05
100 RCL 10

101 *
102 -
103 RCL 15
104 *
105 STO 12
106 RCL 05
107 RCL 09
108 *
109 RCL 06
110 RCL 08
111 *
112 -
113 RCL 16
114 *
115 STO 13
116 X^2
117 RCL 12
118 X^2
119 +
120 RCL 11
121 X^2
122 +
123 SQRT
124 ST/ 11
125 ST/ 12
126 ST/ 13
127 CLX
128 STO 19
129 STO 22
130 STO 18
131 STO 20
132 SIGN
133 STO 17
134 STO 21
135 RCL 11
136 RCL 12
137 R-P
138 X=0?
139 GTO 00
140 RCL 12
141 STO 17
142 RCL 11
143 CHS
144 STO 18
145 R-P
146 ST/ 17
147 ST/ 18
148 RCL 13
149 RCL 18
150 *

151 CHS
152 STO 20
153 RCL 13
154 RCL 17
155 *
156 STO 21
157 RCL 11
158 RCL 18
159 *
160 RCL 12
161 RCL 17
162 *
163 -
164 STO 22
165 LBL 00
166 RCL 14
167 ST* 17
168 ST* 20
169 RCL 15
170 ST* 18
171 ST* 21
172 RCL 16
173 ST* 19
174 ST* 22
175 RCL 17
176 RCL 20
177 *
178 RCL 18
179 RCL 21
180 *
181 +
182 RCL 19
183 RCL 22
184 *
185 +
186 ST+ X
187 RCL 17
188 X^2
189 RCL 18
190 X^2
191 +
192 RCL 19
193 X^2
194 +
195 RCL 20
196 X^2
197 -
198 RCL 21
199 X^2
200 -

201 RCL 22
202 X^2
203 -
204 R-P
205 CLX
206 2
207 /
208 STO 00
209 XEQ 02
210 STO 23
211 RDN
212 STO 24
213 X<>Y
214 STO 25
215 RCL 00
216 90
217 +
218 XEQ 02
219 STO 26
220 RDN
221 STO 27
222 X<>Y
223 STO 28
224 GTO 01
225 LBL 02
226 1
227 P-R
228 RCL 17
229 X<>Y
230 *
231 STO 26
232 CLX
233 RCL 18
234 LASTX
235 *
236 STO 27
237 CLX
238 RCL 19
239 LASTX
240 *
241 STO 28
242 CLX
243 RCL 20
244 X<>Y
245 *
246 ST+ 26
247 RCL 21
248 LASTX
249 *
250 ST+ 27

251 RCL 22
252 LASTX
253 *
254 ST+ 28
255 RCL 28
256 RCL 27
257 RCL 26
258 RTN
259 LBL 01
260 RCL 26
261 X^2
262 RCL 27
263 X^2
264 +
265 RCL 28
266 X^2
267 +
268 SQRT
269 STO 30
270 RCL 23
271 X^2
272 RCL 24
273 X^2
274 +
275 RCL 25
276 X^2
277 +
278 SQRT
279 STO 29
280 X>Y?
281 GTO 00
282 X<> 30
283 STO 29
284 RCL 23
285 X<> 26
286 STO 23
287 RCL 24
288 X<> 27
289 STO 24
290 RCL 25
291 X<> 28
292 STO 25
293 LBL 00
294 RCL 29
295 ST/ 23
296 ST/ 24
297 ST/ 25
298 RCL 30
299 ST/ 26
300 ST/ 27
301 ST/ 28

302 RCL 05
303 RCL 23
304 *
305 RCL 06
306 RCL 24
307 *
308 +
309 RCL 07
310 RCL 25
311 *
312 +
313 STO 01
314 RCL 05
315 RCL 26
316 *
317 RCL 06
318 RCL 27
319 *
320 +
321 RCL 07
322 RCL 28
323 *
324 +
325 STO 02
326 RCL 08
327 RCL 23
328 *
329 RCL 09
330 RCL 24
331 *
332 +
333 RCL 10
334 RCL 25
335 *
336 +
337 STO 03
338 RCL 08
339 RCL 26
340 *
341 RCL 09
342 RCL 27
343 *
344 +
345 RCL 10
346 RCL 28
347 *
348 +
349 STO 04
350 RCL 30
351 RCL 29
352 END

( 515 bytes / SIZE 031 )

STACK	INPUTS	OUTPUTS
T	L₁ ( ° ' " )	/
Z	b₁ ( ° ' " )	/
Y	L₂ ( ° ' " )	B ( km )
X	b₂ ( ° ' " )	A ( km )

Where A & B are the semi-axes of the ellipse.

Example: Calculate the elements of the great ellipse defined by the U.S. Naval Observatory at Washington (D.C.) :

L₁ = 77°03'56"0 W ; b₁ = 38°55'17"2 N and the Observatoire de Paris: L₂ = 2°20'13"8 E ; b₂ = 48°50'11"2 N

   -77.0356    ENTER^
    38.55172 ENTER^
     2.20138   ENTER^
    48.50112 XEQ "GEL"   >>>> a' = 6378.104423 km                            ---Execution time = 24s---
                                          X<>Y b' = 6364.794129 km

-The coordinates of the U.S. Naval Observatory at Washington (D.C.) are in R05-R06-R07: ( 2322.282431 , -4392.706835 , 3985.543384 )

and with respect to ( U , V ) are in R01 & R02: ( -3881.641327 , -5050.374807 )

-The coordinates of the Observatoire de Paris Observatoire de Paris are in R08-R09-R10: ( 4016.636531 , 1248.414112 , 4778.596908 )

and with respect to ( U , V ) are in R03-R04: ( 1976.680334 , -6051.415542 )

-The unit vector corresponding to the major axis is R23-R24-R25 = U = ( 0.186195902 , 0.982512530 , 0.000463853 )
-The unit vector corresponding to the minor axis is R26-R27-R28 = V = ( -0.602931120 , 0.114634184 , -0.789514451 )

Notes:

-This program employs a triaxial ellipsoidal model of the Earth with a = 6378.172 km, b = 6378.102 km , c = 6356.752314 km
and longitude of the major axis = -14.929°
-Change the corresponding lines ( 47-44-40-07 ) if you want to use other constants.

-The longitudes are positive East.

Reference:

[1] https://en.wikipedia.org/wiki/Ellipsoid

hp41programs

Ellipsoid-Plane Intersection for the HP-41