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) :  x2 / a2 + y2 / b2 + z2 / c2 = 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 = f0           R16-R17-R18 = Vertex1         R25-R26-R27 = Vertex4
•  R02 = n      •  R05 = b               R10-R11-R12 = f1           R19-R20-R21 = Vertex2         R28 = t0
•  R03 = p      •  R06 = c               R13-R14-R15 = f2           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) :  x2 / 25 + y2 / 16 + z2 / 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 ) = f0 =  center of the ellipse
R10-R11-R12 = ( 3.7153 , -2.4769 ,  0 )        = f1
R13-R14-R15 = ( 1.8862 , 1.8107 , -2.3011 ) = f2

-So any point of the ellipse may be obtained by the formula:  f0 +  f1 Cos t + f2 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 ) = V1
R19-R20-R21 = ( 1.4021 , 3.0561 , -1.7432 )  = V2
R22-R23-R24 = ( -0.1135 , -1.8190 , 2.6710 ) = V3
R25-R26-R27 = ( -3.4528 , 2.4818 , 1.1150 )  = V4

-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 =  t0     R14 = 6378.172   R15 = 6378.102   R16 = 6356.752314

-> R01 = xM   R03 = xN      R05-R06-R07 = XM YM ZM   R11-R12-R13: temp    R23-R24-R25 = U = unit vector of the 1st vertex
-> R02 = yM   R04 = yN      R08-R09-R10 = XN YN ZN     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 L1 ( ° ' " ) / Z b1 ( ° ' " ) / Y L2 ( ° ' " ) B ( km ) X b2 ( ° ' " ) 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.) :

L1 = 77°03'56"0 W ;  b1 = 38°55'17"2 N  and the Observatoire de Paris:  L2 = 2°20'13"8 E ;  b2 = 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