1 /* 2 3 -Procedure inelpl_c ( Intersection of ellipse and plane ) 4 5 -Abstract 6 7 Find the intersection of an ellipse and a plane. 8 9 -Disclaimer 10 11 THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE 12 CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. 13 GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE 14 ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE 15 PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" 16 TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY 17 WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A 18 PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC 19 SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE 20 SOFTWARE AND RELATED MATERIALS, HOWEVER USED. 21 22 IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA 23 BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT 24 LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, 25 INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, 26 REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE 27 REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. 28 29 RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF 30 THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY 31 CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE 32 ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. 33 34 -Required_Reading 35 36 ELLIPSES 37 PLANES 38 39 -Keywords 40 41 ELLIPSE 42 GEOMETRY 43 MATH 44 45 */ 46 #include <math.h> 47 #include "SpiceUsr.h" 48 #include "SpiceZfc.h" 49 #include "SpiceZst.h" 50 #include "SpiceZim.h" 51 #undef inelpl_c 52 53 inelpl_c(ConstSpiceEllipse * ellips,ConstSpicePlane * plane,SpiceInt * nxpts,SpiceDouble xpt1[3],SpiceDouble xpt2[3])54 void inelpl_c ( ConstSpiceEllipse * ellips, 55 ConstSpicePlane * plane, 56 SpiceInt * nxpts, 57 SpiceDouble xpt1[3], 58 SpiceDouble xpt2[3] ) 59 /* 60 61 -Brief_I/O 62 63 Variable I/O Description 64 -------- --- -------------------------------------------------- 65 ellips I A CSPICE ellipse. 66 plane I A CSPICE plane. 67 nxpts O Number of intersection points of plane and ellipse. 68 xpt1, 69 xpt2 O Intersection points. 70 71 -Detailed_Input 72 73 ellips is a CSPICE ellipse. The ellipse is allowed to 74 be degenerate: one or both semi-axes may have 75 zero length. 76 77 plane is a CSPICE plane. The intersection of plane 78 and ellipse is sought. 79 80 -Detailed_Output 81 82 nxpts is the number of points of intersection of the 83 geometric plane and ellipse represented by `plane' and 84 `ellips'. `nxpts' may take the values 0, 1, 2 or -1. 85 The value -1 indicates that the ellipse consists of 86 more than one point lies in the plane, so the number 87 of intersection points is infinite. 88 89 When the ellipse consists of a single point and 90 lies in the plane, `nxpts' is set to 1. 91 92 xpt1, 93 xpt2 are the points of intersection of the input plane 94 and ellipse. If there is only one intersection 95 point, both xpt1 and xpt2 contain that point. If 96 the number of intersection points is zero or 97 infinite, the contents of xpt1 and xpt2 are 98 undefined. 99 100 -Parameters 101 102 None. 103 104 -Exceptions 105 106 1) The input plane must be a CSPICE plane: the normal vector must 107 be non-zero and the constant must be non-negative. 108 If the input plane is invalid, the error SPICE(INVALIDPLANE) 109 will be signaled. 110 111 2) If the input ellipse has non-orthogonal axes, the error 112 SPICE(INVALIDELLIPSE) will be signaled. 113 114 3) The input ellipse is allowed to be a line segment or a point; 115 these cases are not considered to be errors. If the ellipse 116 consists of a single point and lies in the plane, the number 117 of intersection points is set to 1 (rather than -1) and 118 the output arguments `xpt1' and `xpt2' are assigned the value 119 of the ellipse's center. 120 121 -Files 122 123 None. 124 125 -Particulars 126 127 This routine computes the intersection set of a non-degenerate 128 plane with a possibly degenerate ellipse. The ellipse is allowed 129 to consist of a line segment or a point. 130 131 A plane may intersect an ellipse in 0, 1, 2, or infinitely many 132 points. For there to be an infinite set of intersection points, 133 the ellipse must lie in the plane and consist of more than one 134 135 -Examples 136 137 1) If we want to find the angle of some ray above the limb of an 138 ellipsoid, where the angle is measured in a plane containing 139 the ray and a "down" vector, we can follow the procedure 140 given below. We assume the ray does not intersect the 141 ellipsoid. The result we seek is called angle, imaginatively 142 enough. 143 144 We assume that all vectors are given in body-fixed 145 coordinates. 146 147 #include "SpiceUsr.h" 148 . 149 . 150 . 151 /. 152 Find the limb of the ellipsoid as seen from the 153 point observ. Here a, b, and c are the lengths of 154 the semi-axes of the ellipsoid. The limb is 155 returned as a SpiceEllipse. 156 ./ 157 158 edlimb_c ( a, b, c, observ, &limb ); 159 160 /. 161 The ray direction vector is raydir, so the ray is the 162 set of points 163 164 observ + t * raydir 165 166 where t is any non-negative real number. 167 168 The `down' vector is just -observ. The vectors 169 observ and raydir are spanning vectors for the plane 170 we're interested in. We can use psv2pl_c to represent 171 this plane by a CSPICE plane. 172 ./ 173 psv2pl_c ( observ, observ, raydir, &plane ); 174 175 /. 176 Find the intersection of the plane defined by observ 177 and raydir with the limb. 178 ./ 179 inelpl_c ( limb, plane, nxpts, xpt1, xpt2 ); 180 181 /. 182 We always expect two intersection points, if the vector 183 down is valid. 184 ./ 185 if ( nxpts < 2 ) 186 { 187 [ do something about the error ] 188 } 189 190 /. 191 Form the vectors from observ to the intersection 192 points. Find the angular separation between the 193 boresight ray and each vector from observ to the 194 intersection points. 195 ./ 196 vsub_c ( xpt1, observ, vec1 ); 197 vsub_c ( xpt2, observ, vec2 ); 198 199 sep1 = vsep_c ( vec1, raydir ); 200 sep2 = vsep_c ( vec2, raydir ); 201 202 /. 203 The angular separation we're after is the minimum of 204 the two separations we've computed. 205 ./ 206 angle = mind_c ( 2, sep1, sep2 ); 207 208 209 -Restrictions 210 211 None. 212 213 -Literature_References 214 215 None. 216 217 -Author_and_Institution 218 219 N.J. Bachman (JPL) 220 221 -Version 222 223 -CSPICE Version 2.1.0, 07-OCT-2011 (NJB) 224 225 Relaxed ellipse semi-axes orthogonality test limit 226 SEPLIM from 1.D-12 TO 1.D-9 radians. The angular 227 separation of the axes of the input ellipse must not 228 differ from pi/2 radians by more than this limit. 229 230 -CSPICE Version 2.0.0, 14-JAN-2008 (NJB) 231 232 Bug fix: the routine's specification and behavior have been 233 updated so the routine now returns a meaningful result for the 234 case of an ellipse consisting of a single point. 235 236 Bug fix: in the degenerate case where the input ellipse is a 237 line segment of positive length, and this segment intersects 238 the plane, the number of intersection points is set to 1 239 rather than 2. 240 241 Invalid input planes and ellipses are now diagnosed. 242 243 -CSPICE Version 1.0.0, 28-AUG-2001 (NJB) 244 245 -Index_Entries 246 247 intersection of ellipse and plane 248 249 -& 250 */ 251 252 253 /* 254 -Revisions 255 256 -CSPICE Version 2.0.0, 14-JAN-2008 (NJB) 257 258 Bug fix: the routine's specification and behavior have been 259 updated so the routine now returns a meaningful result for the 260 case of an ellipse consisting of a single point. In this case, 261 if an intersection is found, the number of intersection points 262 is set to 1 and both intersection arguments are set equal to 263 the ellipse's center. 264 265 Bug fix: in the degenerate case where the input ellipse is a 266 line segment of positive length, and this segment intersects 267 the plane, the number of intersection points is set to 1 268 rather than 2. 269 270 Invalid input planes and ellipses are now diagnosed. 271 Error handling code has been added to trap errors that had 272 been erroneously passed off to lower level routines for 273 diagnosis. 274 -& 275 */ 276 277 278 279 { /* Begin inelpl_c */ 280 281 282 /* 283 Local constants 284 */ 285 #define SEPLIM ( 1.0e-9 ) 286 287 /* 288 Local variables 289 */ 290 SpiceDouble alpha; 291 SpiceDouble angle1; 292 SpiceDouble angle2; 293 SpiceDouble beta; 294 SpiceDouble center [3]; 295 SpiceDouble constant; 296 SpiceDouble inpcon; 297 SpiceDouble normal [3]; 298 SpiceDouble point [3]; 299 SpiceDouble sep; 300 SpiceDouble smajor [3]; 301 SpiceDouble sminor [3]; 302 SpiceDouble v [2]; 303 304 SpicePlane trans; 305 306 307 308 /* 309 Participate in error tracing. 310 */ 311 chkin_c ( "inelpl_c" ); 312 313 314 /* 315 Check the input plane. 316 */ 317 pl2nvc_c ( plane, normal, &inpcon ); 318 319 if ( vzero_c(normal) ) 320 { 321 setmsg_c ( "Input SPICE plane has zero normal vector." ); 322 sigerr_c ( "SPICE(INVALIDPLANE)" ); 323 chkout_c ( "inelpl_c" ); 324 return; 325 } 326 else if ( inpcon < 0.0 ) 327 { 328 setmsg_c ( "Input SPICE plane has non-positive " 329 "constant #. Properly constructed " 330 "SPICE planes always have non-negative " 331 "constants." ); 332 errdp_c ( "#", inpcon ); 333 sigerr_c ( "SPICE(INVALIDPLANE)" ); 334 chkout_c ( "inelpl_c" ); 335 return; 336 } 337 338 /* 339 Get the components of the input ellipse; check for 340 invalid semi-axes. The semi-axes may have zero length 341 but they must always be orthogonal. We require this 342 check only if both semi-axes have non-zero length. 343 */ 344 el2cgv_c ( ellips, center, smajor, sminor ); 345 346 if ( !vzero_c(sminor) ) 347 { 348 sep = vsep_c( smajor, sminor ); 349 350 if ( fabs( sep-halfpi_c() ) > SEPLIM ) 351 { 352 setmsg_c ( "Input SPICE ellipse has non-orthogonal " 353 "semi-axes: (#,#,#) and (#,#,#). Angular " 354 "separation of these vectors is # radians. " 355 "Properly constructed SPICE ellipses " 356 "always have orthogonal semi-axes." ); 357 errdp_c ( "#", smajor[0] ); 358 errdp_c ( "#", smajor[1] ); 359 errdp_c ( "#", smajor[2] ); 360 errdp_c ( "#", sminor[0] ); 361 errdp_c ( "#", sminor[1] ); 362 errdp_c ( "#", sminor[2] ); 363 errdp_c ( "#", sep ); 364 sigerr_c ( "SPICE(INVALIDELLIPSE)" ); 365 chkout_c ( "inelpl_c" ); 366 return; 367 } 368 } 369 370 /* 371 If the input ellipse is a single point, decide now 372 whether the ellipse lies in the plane. 373 */ 374 375 if ( vzero_c(smajor) ) 376 { 377 /* 378 The ellipse is a single point. If the ellipse's center 379 lies in the plane, the whole ellipse is the one 380 intersection point. Check the inner product of the 381 center and the plane's normal vector. 382 */ 383 384 if ( vdot_c(center, normal) == inpcon ) 385 { 386 /* 387 The center does in fact lie in the plane. 388 */ 389 390 *nxpts = 1; 391 392 vequ_c ( center, xpt1 ); 393 vequ_c ( center, xpt2 ); 394 } 395 else 396 { 397 /* 398 There's no intersection: the intersection arguments 399 are left undefined in this case. 400 */ 401 402 *nxpts = 0; 403 } 404 405 /* 406 Return now; this simplifies the logic to follow. 407 */ 408 chkout_c ( "inelpl_c" ); 409 return; 410 } 411 412 /* 413 At this point the ellipse may still be degenerate: it can be a 414 line segment. We'll need to compute the intersection point or 415 points if we have a positive, finite intersection set. 416 417 The first thing we want to do is translate the plane and the 418 ellipse so as to center the ellipse at the origin. To translate 419 the plane, just get a point and normal vector, and translate 420 the point. Find the plane constant of the translated plane. 421 */ 422 pl2nvp_c ( plane, normal, point ); 423 vsub_c ( point, center, point ); 424 nvp2pl_c ( normal, point, &trans ); 425 pl2nvc_c ( &trans, normal, &constant ); 426 427 /* 428 Ok, we can get to work. The locus of the ellipse is 429 430 cos(theta) smajor + sin(theta) sminor, 431 432 and any point x of the ellipse that intersects the input plane 433 satisfies 434 435 < x, normal > = constant. 436 437 Substituting our expression for points on the ellipse into the 438 second equation, we arrive at 439 440 cos(theta) < smajor, normal > 441 + sin(theta) < sminor, normal > = constant. (1) 442 443 This equation merits a little analysis. First, if `normal' 444 is orthogonal to `smajor' and `sminor, the plane and ellipse must 445 be parallel. Also, the left side of the equation is zero in 446 this case. If `constant' is non-zero, there are no solutions: 447 the ellipse and plane are parallel but do not intersect. If 448 `constant' is zero, the ellipse lies in the plane: all values of 449 theta are solutions. Let's get this case out of the way 450 right now, shall we? 451 */ 452 v[0] = vdot_c ( smajor, normal ); 453 v[1] = vdot_c ( sminor, normal ); 454 455 /* 456 Test whether the plane and ellipse are parallel. 457 */ 458 if ( vzerog_c( v, 2 ) ) 459 { 460 /* 461 The ellipse lies in the plane if and only if constant is zero. 462 In any case, we don't modify xpt1 or xpt2. 463 */ 464 if ( constant == 0.0 ) 465 { 466 *nxpts = -1; 467 } 468 else 469 { 470 *nxpts = 0; 471 } 472 473 chkout_c ( "inelpl_c" ); 474 return; 475 } 476 477 478 /* 479 Now if `normal' is not orthogonal to both `smajor' and `sminor', 480 the vector 481 482 v = ( < smajor, normal >, < sminor, normal > ) 483 484 is non-zero. We can re-write (1) as 485 486 < u, v > = constant, 487 488 where 489 490 u = ( cos(theta), sin(theta) ). 491 492 If alpha is the angle between u and v, we have 493 494 < u, v > = || u || * || v || * cos(alpha), 495 496 so 497 498 || v || * cos(alpha) = constant. (2) 499 500 `constant' is positive, since pl2nvc_c returns the distance 501 between its input plane and the origin as the output 502 plane constant. 503 504 Equation (2) has solutions if and only if 505 506 || v || > constant. (3) 507 - 508 509 510 Let's return right now if there are no solutions. 511 */ 512 if ( vnormg_c ( v, 2 ) < constant ) 513 { 514 *nxpts = 0; 515 516 chkout_c ( "inelpl_c" ); 517 return; 518 } 519 520 521 /* 522 Since (3) above is satisfied, the plane and ellipse intersect. 523 We can find alpha by the formula 524 525 alpha = + arccos ( constant / || v || ) 526 527 Since `alpha' is the angular separation between `u' and `v', we 528 can find `u' once we have the angular position of `v'; let's 529 call that `beta'. The angular position of `u'(which we called 530 `theta' earlier) will be 531 532 theta = beta + alpha. 533 - 534 535 The values of `theta' are the angles we seek. 536 */ 537 alpha = acos ( constant / vnormg_c ( v, 2 ) ); 538 539 beta = atan2 ( v[1], v[0] ); 540 541 angle1 = beta - alpha; 542 angle2 = beta + alpha; 543 544 /* 545 Determine the number of intersection points. We have a special 546 case if the semi-minor axis has length zero: in that case `beta' is 547 zero or Pi, and although `angle1' and `angle2' may differ, the 548 cosines of these angles are identical. Since in this case 549 the solutions corresponding to `angle1' and `angle2' have the 550 form 551 552 center + cos(angle1)*smajor 553 center + cos(angle2)*smajor 554 555 the solutions are identical. 556 */ 557 558 if ( vzero_c(sminor) ) 559 { 560 *nxpts = 1; 561 } 562 else 563 { 564 if ( angle1 == angle2 ) 565 { 566 /* 567 This case occurs when `alpha' is zero. 568 */ 569 570 *nxpts = 1; 571 } 572 else 573 { 574 *nxpts = 2; 575 } 576 } 577 578 /* 579 Compute the intersection points. 580 */ 581 vlcom3_c ( 1.0, center, 582 cos(angle1), smajor, 583 sin(angle1), sminor, xpt1 ); 584 585 vlcom3_c ( 1.0, center, 586 cos(angle2), smajor, 587 sin(angle2), sminor, xpt2 ); 588 589 chkout_c ( "inelpl_c" ); 590 591 } /* End inelpl_c */ 592 593