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)  x2 / a2 + y2 / b2 + z2 / c2 = 1

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

-The vector  N = ( x / a2 , y / b2 , z / c2 ) 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:

( b2 - a2 ) Sin u Cos u Cos v + a X Sin u - b Y Cos u = 0
( b2 - c2 ) 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 XY 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):   x2 / 49 + y2 / 36 + z2 / 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   x12 / a12 + ............... + xn2 / an2 = 1  and a point  P( X1 , ............... , Xn )

writing that  N( x1 / a12 , .............. , xn / an2 )  and  PQ ( X1 - x1 , ................ , Xn - xn )  are colinear gives - after squaring the equalities:

( x12 / a12 ) / [ a12 ( X1 - x1 )2 ] = ...................... = ( xn2 / an2 ) / [ an2 ( Xn - xn )2 ] = 1 / [ a12 ( X1 - x1 )2 + ........... + an2 ( Xn - xn )2 ]

-Let  A = SQRT [ a12 ( X1 - x1 )2 + ........... + an2 ( Xn - xn )2 ]

-We get the system:

x1 = a12 X1 / ( a12 + A ) , ................... , xn = an2 Xn / ( an2 + 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 = a1      •  R02 = a2    ..........................   •  Rnn = an
•  Rn+1 = X1   •  Rn+2 = X2 ..........................   •  R2n = Xn

•  R2n+1 = x1  •  R2n+2 = x2 .........................   •  R3n = xn
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( x1 , ......... , xn )

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

(E):   x2 / 49 + y2 / 36 + z2 / 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):   x2 / 49 + y2 / 36 + z2 / 25 + t2 / 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?