Point-Ellipsoid

Distance Point-Ellipsoid for the HP-41

Overview

1°) Distance from a Point to an Ellipsoid
2°) Distance from a Point to a HyperEllipsoid

1°) Distance from a Point P(X,Y,Z) to an Ellipsoid

"DPE" uses the parametric coordinates (u,v) of a point Q(x,y,z) on the ellipsoid (E) x² / a² + y² / b² + z² / c² = 1

    x/a = Cos u Cos v
    y/b = Sin u Cos v
    z/c = Sin v

-The vector N = ( x / a² , y / b² , z / c² ) is normal to the ellipsoid (E)
-Writing that N and PQ are colinear yields:

Cos u Cos v / [ a ( a Cos u Cos v - X ] = Sin u Cos v / [ b ( b Sin u Cos v - Y ] = Sin v / [ c ( c Sin v - Z ) ]

-From which we deduce:

( b² - a² ) Sin u Cos u Cos v + a X Sin u - b Y Cos u = 0
( b² - c² ) Sin u Sin v Cos v + c Z Sin u Cos v - b Y Sin v = 0

-Then "DPE" solves this non-linear system by Newton's method

Data Registers: R00 temp ( Registers R01 thru R03 are to be initialized before executing "DPE" )

                                      • R01 = a     R04 = X     R07 = x    R10 = d     R11 = u        R13 to R20: temp
                                      • R02 = b     R05 = Y     R08 = y                      R12 = v
                                      • R03 = c     R06 = Z      R09 = z
Flags: /
Subroutines: /

01 LBL "DPE"
02 DEG
03 STO 04
04 RDN
05 STO 05
06 X<>Y
07 STO 06
08 RCL 04
09 RCL 01
10 /
11 X^2
12 RCL 05
13 RCL 02
14 /
15 X^2
16 +
17 SQRT
18 RCL 03
19 *
20 R-P
21 X<>Y
22 STO 12
23 RCL 02
24 RCL 05
25 *
26 RCL 01
27 RCL 04
28 *
29 R-P
30 X<>Y
31 STO 11
32 RCL 02
33 RCL 02

34 RCL 01
35 ST- Z
36 +
37 *
38 STO 13
39 RCL 02
40 RCL 02
41 RCL 03
42 ST- Z
43 +
44 *
45 STO 14
46 LBL 01
47 VIEW 11
48 RCL 11
49 1
50 P-R
51 STO 08
52 X<>Y
53 STO 07
54 RCL 12
55 LASTX
56 P-R
57 STO 10
58 X<>Y
59 STO 09
60 RCL 14
61 *
62 RCL 03
63 RCL 06
64 *
65 +
66 RCL 07

67 RCL 10
68 *
69 *
70 RCL 02
71 RCL 05
72 *
73 RCL 09
74 *
75 -
76 STO 16
77 RCL 13
78 RCL 08
79 *
80 RCL 10
81 *
82 RCL 01
83 RCL 04
84 *
85 +
86 RCL 07
87 *
88 RCL 02
89 RCL 05
90 *
91 RCL 08
92 *
93 -
94 STO 15
95 RCL 13
96 RCL 08
97 X^2
98 RCL 07
99 X^2

100 -
101 *
102 RCL 10
103 *
104 RCL 01
105 RCL 04
106 *
107 RCL 08
108 *
109 +
110 RCL 02
111 RCL 05
112 *
113 RCL 07
114 *
115 +
116 STO 17
117 RCL 13
118 RCL 07
119 *
120 RCL 08
121 *
122 RCL 09
123 *
124 CHS
125 STO 18
126 RCL 14
127 RCL 09
128 *
129 RCL 03
130 RCL 06
131 *
132 +

133 RCL 08
134 RCL 10
135 *
136 *
137 STO 19
138 RCL 14
139 RCL 10
140 X^2
141 RCL 09
142 X^2
143 -
144 *
145 RCL 03
146 RCL 06
147 *
148 RCL 09
149 *
150 -
151 RCL 07
152 *
153 RCL 02
154 RCL 05
155 *
156 RCL 10
157 *
158 -
159 STO 20
160 RCL 17
161 *
162 RCL 18
163 RCL 19
164 *
165 -

166 STO 00
167 X=0?
168 GTO 02
169 RCL 16
170 RCL 18
171 *
172 RCL 15
173 RCL 20
174 *
175 -
176 X<>Y
177 /
178 R-D
179 ST+ 11
180 RCL 15
181 RCL 19
182 *
183 RCL 16
184 RCL 17
185 *
186 -
187 RCL 00
188 /
189 R-D
190 ST+ 12
191 R-P
192 2 E-7
193 X<Y?
194 GTO 01
195 LBL 02
196 RCL 12
197 1
198 P-R
199 RCL 11

200 X<>Y
201 P-R
202 RCL 01
203 *
204 STO 07
205 CLX
206 RCL 02
207 *
208 STO 08
209 CLX
210 RCL 03
211 *
212 STO 09
213 RCL 06
214 -
215 X^2
216 RCL 05
217 RCL 08
218 -
219 X^2
220 +
221 RCL 04
222 RCL 07
223 -
224 X^2
225 +
226 SQRT
227 RCL 09
228 RCL 08
229 RCL 07
230 R^
231 STO 10
232 CLD
233 END

( 266 bytes / SIZE 021 )

STACK	INPUTS	OUTPUTS
T	/	z
Z	Z	y
Y	Y	x
X	X	d

Example1: (E): x² / 49 + y² / 36 + z² / 25 = 1 P(10,11,12)

    7 STO 01
    6 STO 02
    5 STO 03

    12 ENTER^
    11 ENTER^
    10 XEQ "DPE"   >>>>    d = 13.23891299   = R10                   ---Execution time = 32s---
                                RDN     x = 3.896191634 = R07
                                RDN     y = 3.511764224 = R08
                                RDN     z = 2.948002121 = R09

Example2: P(0.1,0.2,0.3)

     0.3   ENTER^
     0.2   ENTER^
     0.1      R/S    >>>>        d = 4.691015000                        ---Execution time = 53s---
                              RDN     x = 0.192100176
                              RDN     y = 0.575653468
                              RDN     z = 4.975042648

Notes:

"DPE" works well if P is outside the ellipsoid
Otherwise, the results may be doubtful, for instance, with P(0,0,0) it yields d = 7 which is obviously wrong ( d = 5 )

-You can store other initial values for u & v in R11 & R12 and press XEQ 01 again

2°) Distance from a Point to a HyperEllipsoid

-With a hyperellipsoid x₁² / a₁² + ............... + x_n² / a_n² = 1 and a point P( X₁ , ............... , X_n )

writing that N( x₁ / a₁² , .............. , x_n / a_n² ) and PQ ( X₁ - x₁ , ................ , X_n - x_n ) are colinear gives - after squaring the equalities:

( x₁² / a₁² ) / [ a₁² ( X₁ - x₁ )² ] = ...................... = ( x_n² / a_n² ) / [ a_n² ( X_n - x_n )² ] = 1 / [ a₁² ( X₁ - x₁ )² + ........... + a_n² ( X_n - x_n )² ]

-Let A = SQRT [ a₁² ( X₁ - x₁ )² + ........... + a_n² ( X_n - x_n )² ]

-We get the system:

x₁ = a₁² X₁ / ( a₁² + A ) , ................... , x_n = a_n² X_n / ( a_n² + A )

-"DPHE" solves this non-linear system by a successive approximation method ( fixed-point theorem ).

Data Registers: • R00 = n ( Registers R00 thru R2n are to be initialized before executing "DPHE" )

                                      • R01 = a₁      • R02 = a₂    ..........................   • Rnn = a_n
                                      • Rn+1 = X₁   • Rn+2 = X₂ ..........................   • R2n = X_n

                                      • R2n+1 = x₁ • R2n+2 = x₂ .........................   • R3n = x_n
Flags: /
Subroutines: /

01 LBL "DPHE"
02 RCL 00
03 3
04 *
05 .1
06 %
07 +
08 RCL 00
09 -
10 ISG X
11 CLRGX
12 LBL 01
13 RCL 00
14 STO M
15 ENTER
16 ST+ X
17 STO N
18 +

19 STO O
20 CLX
21 LBL 02
22 RCL IND N
23 RCL IND O
24 -
25 RCL IND M
26 *
27 X^2
28 +
29 DSE O
30 DSE N
31 DSE M
32 GTO 02
33 SQRT
34 RCL 00
35 STO M
36 ST+ N

37 ST+ O
38 CLX
39 LBL 03
40 RCL Y
41 RCL IND M
42 X^2
43 ST+ Y
44 X<>Y
45 /
46 RCL IND N
47 *
48 ENTER
49 X<> IND O
50 -
51 ABS
52 +
53 DSE O
54 DSE N

55 DSE M
56 GTO 03
57 VIEW X
58 X#0?
59 GTO 01
60 RCL 00
61 STO M
62 ST+ N
63 ST+ O
64 CLX
65 CLD
66 LBL 04
67 RCL IND N
68 RCL IND O
69 -
70 X^2
71 +
72 DSE O

73 DSE N
74 DSE M
75 GTO 04
76 SQRT
77 RCL 00
78 3
79 *
80 .1
81 %
82 +
83 RCL 00
84 -
85 ISG X
86 X<>Y
87 CLA
88 END

( 136 bytes / SIZE 3n+1 )

STACK	INPUTS	OUTPUTS
Y	/	bbb.eee
X	/	/

Where bbb.eee = 2n+1.3n = control number of Q( x₁ , ......... , x_n )

Example1: Let's take the same example as above ( 3-dimensional )

(E): x² / 49 + y² / 36 + z² / 25 = 1 P(10,11,12)

    3 STO 00

    7 STO 01    10 STO 04
    6 STO 02    11 STO 05
    5 STO 03    12 STO 06

XEQ "DPHE" >>>> d = 13.23891300 ---Execution time = 74s---
X<>Y 7.009

-And you get x, y , z in R07-R08-R09:

      x = 3.896191632 = R07
      y = 3.511764223 = R08
      z = 2.948002118 = R09

Example2: (E): x² / 49 + y² / 36 + z² / 25 + t² / 16 = 1 P(10,11,12,13)

     4 STO 00

     7 STO 01     10 STO 05
     6 STO 02     11 STO 06
     5 STO 03     12 STO 07
     4 STO 04     13 STO 08

XEQ "DPHE" >>>> d = 17.81606768 ---Execution time = 85s---
X<>Y 9.012

-And

      x = 3.464440290 = R09
      y = 3.083223967 = R10
      z = 2.554561428 = R11
      t = 1.918164674 = R12

Note:

-With the same hyperellipsoid and P(4,4,3,3), the result

d = 1.528630348 is returned in 6mn32s but in less than 1 second with V41.

Warning:

-The HP41 displays a series of positive numbers which tend to zero.
-This program returns zero if P is inside (E)

-"DPHE" works well if P is far enough from the surface of the (hyper)ellipsoid.
-But if P is very close to (E), the execution time becomes prohibitive, even with V41 !
-If you know a method to accelerate the convergence...

-The test line 58 X#0? could perhaps lead to an infinite loop
-So, line 58 could be replaced by something like E-8 X<Y?

hp41programs

Distance Point-Ellipsoid for the HP-41