1 /* inrypl.f -- translated by f2c (version 19980913).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
8 /* Table of constant values */
9
10 static integer c__3 = 3;
11 static doublereal c_b16 = 1.;
12
13 /* $Procedure INRYPL ( Intersection of ray and plane ) */
inrypl_(doublereal * vertex,doublereal * dir,doublereal * plane,integer * nxpts,doublereal * xpt)14 /* Subroutine */ int inrypl_(doublereal *vertex, doublereal *dir, doublereal *
15 plane, integer *nxpts, doublereal *xpt)
16 {
17 /* System generated locals */
18 doublereal d__1, d__2;
19
20 /* Local variables */
21 doublereal udir[3];
22 extern /* Subroutine */ int vhat_(doublereal *, doublereal *), vscl_(
23 doublereal *, doublereal *, doublereal *);
24 extern doublereal vdot_(doublereal *, doublereal *);
25 extern /* Subroutine */ int vequ_(doublereal *, doublereal *);
26 doublereal scale;
27 extern /* Subroutine */ int chkin_(char *, ftnlen);
28 extern doublereal dpmax_(void);
29 extern /* Subroutine */ int vlcom_(doublereal *, doublereal *, doublereal
30 *, doublereal *, doublereal *);
31 doublereal const__, prjvn;
32 extern doublereal vnorm_(doublereal *);
33 extern logical vzero_(doublereal *);
34 extern /* Subroutine */ int pl2nvc_(doublereal *, doublereal *,
35 doublereal *), cleard_(integer *, doublereal *);
36 doublereal mscale, prjdif, sclcon, toobig, normal[3], prjdir;
37 extern logical smsgnd_(doublereal *, doublereal *);
38 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
39 ftnlen), vsclip_(doublereal *, doublereal *), setmsg_(char *,
40 ftnlen);
41 extern logical return_(void);
42 doublereal sclvtx[3];
43
44 /* $ Abstract */
45
46 /* Find the intersection of a ray and a plane. */
47
48 /* $ Disclaimer */
49
50 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
51 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
52 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
53 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
54 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
55 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
56 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
57 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
58 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
59 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
60
61 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
62 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
63 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
64 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
65 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
66 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
67
68 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
69 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
70 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
71 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
72
73 /* $ Required_Reading */
74
75 /* PLANES */
76
77 /* $ Keywords */
78
79 /* GEOMETRY */
80
81 /* $ Declarations */
82 /* $ Brief_I/O */
83
84 /* Variable I/O Description */
85 /* -------- --- -------------------------------------------------- */
86 /* VERTEX, */
87 /* DIR I Vertex and direction vector of ray. */
88 /* PLANE I A SPICELIB plane. */
89 /* NXPTS O Number of intersection points of ray and plane. */
90 /* XPT O Intersection point, if NXPTS = 1. */
91
92 /* $ Detailed_Input */
93
94 /* VERTEX, */
95 /* DIR are a point and direction vector that define a */
96 /* ray in three-dimensional space. */
97
98 /* PLANE is a SPICELIB plane. */
99
100 /* $ Detailed_Output */
101
102 /* NXPTS is the number of points of intersection of the */
103 /* input ray and plane. Values and meanings of */
104 /* NXPTS are: */
105
106 /* 0 No intersection. */
107
108 /* 1 One point of intersection. Note that */
109 /* this case may occur when the ray's */
110 /* vertex is in the plane. */
111
112 /* -1 An infinite number of points of */
113 /* intersection; the ray lies in the plane. */
114
115
116 /* XPT is the point of intersection of the input ray */
117 /* and plane, when there is exactly one point of */
118 /* intersection. Otherwise, XPT is the zero vector. */
119
120 /* $ Parameters */
121
122 /* None. */
123
124 /* $ Exceptions */
125
126 /* 1) If the ray's direction vector is the zero vector, the error */
127 /* SPICE(ZEROVECTOR) is signaled. NXPTS and XPT are not */
128 /* modified. */
129
130
131 /* 2) If the ray's vertex is further than DPMAX() / 3 from the */
132 /* origin, the error SPICE(VECTORTOOBIG) is signaled. NXPTS */
133 /* and XPT are not modified. */
134
135
136 /* 3) If the input plane is further than DPMAX() / 3 from the */
137 /* origin, the error SPICE(VECTORTOOBIG) is signaled. NXPTS */
138 /* and XPT are not modified. */
139
140
141 /* 4) The input plane should be created by one of the SPICELIB */
142 /* routines */
143
144 /* NVC2PL */
145 /* NVP2PL */
146 /* PSV2PL */
147
148 /* Invalid input planes will cause unpredictable results. */
149
150
151 /* 5) In the interest of good numerical behavior, in the case */
152 /* where the ray's vertex is not in the plane, this routine */
153 /* considers that an intersection of the ray and plane occurs */
154 /* only if the distance between the ray's vertex and the */
155 /* intersection point is less than DPMAX() / 3. */
156
157 /* If VERTEX is not in the plane and this condition is not */
158 /* met, then NXPTS is set to 0 and XPT is set to the zero */
159 /* vector. */
160
161 /* $ Files */
162
163 /* None. */
164
165 /* $ Particulars */
166
167 /* The intersection of a ray and plane in three-dimensional space */
168 /* can be a the empty set, a single point, or the ray itself. */
169
170 /* $ Examples */
171
172 /* 1) Find the camera projection of the center of an extended */
173 /* body. For simplicity, we assume: */
174
175 /* -- The camera has no distortion; the image of a point */
176 /* is determined by the intersection of the focal plane */
177 /* and the line determined by the point and the camera's */
178 /* focal point. */
179
180 /* -- The camera's pointing matrix (C-matrix) is available */
181 /* in a C-kernel. */
182
183
184 /* C */
185 /* C Load Leapseconds and SCLK kernels to support time */
186 /* C conversion. */
187 /* C */
188 /* CALL FURNSH ( 'LEAP.KER' ) */
189 /* CALL FURNSH ( 'SCLK.KER' ) */
190
191 /* C */
192 /* C Load an SPK file containing ephemeris data for */
193 /* C observer (a spacecraft, whose NAIF integer code */
194 /* C is SC) and target at the UTC epoch of observation. */
195 /* C */
196 /* CALL FURNSH ( 'SPK.BSP' ) */
197
198 /* C */
199 /* C Load a C-kernel containing camera pointing for */
200 /* C the UTC epoch of observation. */
201 /* C */
202 /* CALL FURNSH ( 'CK.BC' ) */
203
204 /* C */
205 /* C Find the ephemeris time (barycentric dynamical time) */
206 /* C and encoded spacecraft clock times corresponding to */
207 /* C the UTC epoch of observation. */
208 /* C */
209 /* CALL UTC2ET ( UTC, ET ) */
210 /* CALL SCE2C ( SC, ET, SCLKDP ) */
211
212 /* C */
213 /* C Encode the pointing lookup tolerance. */
214 /* C */
215 /* CALL SCTIKS ( SC, TOLCH, TOLDP ) */
216
217 /* C */
218 /* C Find the observer-target vector at the observation */
219 /* C epoch. In this example, we'll use a light-time */
220 /* C corrected state vector. */
221 /* C */
222 /* CALL SPKEZ ( TARGET, ET, 'J2000', 'LT', SC, */
223 /* . STATE, LT ) */
224
225 /* C */
226 /* C Look up camera pointing. */
227 /* C */
228 /* CALL CKGP ( CAMERA, SCLKDP, TOLDP, 'J2000', CMAT, */
229 /* . CLKOUT, FOUND ) */
230
231 /* IF ( .NOT. FOUND ) THEN */
232
233 /* [Handle this case...] */
234
235 /* END IF */
236
237 /* C */
238 /* C Negate the spacecraft-to-target body vector and */
239 /* C convert it to camera coordinates. */
240 /* C */
241 /* CALL VMINUS ( STATE, DIR ) */
242 /* CALL MXV ( CMAT, DIR, DIR ) */
243
244 /* C */
245 /* C If FL is the camera's focal length, the effective */
246 /* C focal point is */
247 /* C */
248 /* C FL * ( 0, 0, 1 ) */
249 /* C */
250 /* CALL VSCL ( FL, ZVEC, FOCUS ) */
251
252 /* C */
253 /* C The camera's focal plane contains the origin in */
254 /* C camera coordinates, and the z-vector is orthogonal */
255 /* C to the plane. Make a SPICELIB plane representing */
256 /* C the focal plane. */
257 /* C */
258 /* CALL NVC2PL ( ZVEC, 0.D0, FPLANE ) */
259
260 /* C */
261 /* C The image of the target body's center in the focal */
262 /* C plane is defined by the intersection with the focal */
263 /* C plane of the ray whose vertex is the focal point and */
264 /* C whose direction is DIR. */
265 /* C */
266 /* CALL INRYPL ( FOCUS, DIR, FPLANE, NXPTS, IMAGE ) */
267
268 /* IF ( NXPTS .EQ. 1 ) THEN */
269 /* C */
270 /* C The body center does project to the focal plane. */
271 /* C Check whether the image is actually in the */
272 /* C camera's field of view... */
273 /* C */
274 /* . */
275 /* . */
276 /* . */
277 /* ELSE */
278
279 /* C */
280 /* C The body center does not map to the focal plane. */
281 /* C Handle this case... */
282 /* C */
283 /* . */
284 /* . */
285 /* . */
286 /* END IF */
287
288
289
290 /* 2) Find the Saturn ring plane intercept of a spacecraft-mounted */
291 /* instrument's boresight vector. We want the find the point */
292 /* in the ring plane that will be observed by an instrument */
293 /* with a give boresight direction at a specified time. We */
294 /* must account for light time and stellar aberration in order */
295 /* to find this point. The intercept point will be expressed */
296 /* in Saturn body-fixed coordinates. */
297
298 /* In this example, we assume */
299
300 /* -- The ring plane is equatorial. */
301
302 /* -- Light travels in a straight line. */
303
304 /* -- The light time correction for the ring plane intercept */
305 /* can be obtained by performing three light-time */
306 /* correction iterations. If this assumption does not */
307 /* lead to a sufficiently accurate result, additional */
308 /* iterations can be performed. */
309
310 /* -- A Newtonian approximation of stellar aberration */
311 /* suffices. */
312
313 /* -- The boresight vector is given in J2000 coordinates. */
314
315 /* -- The observation epoch is ET ephemeris seconds past */
316 /* J2000. */
317
318 /* -- The boresight vector, spacecraft and planetary */
319 /* ephemerides, and ring plane orientation are all known */
320 /* with sufficient accuracy for the application. */
321
322 /* -- All necessary kernels are loaded by the caller of */
323 /* this example routine. */
324
325
326 /* SUBROUTINE RING_XPT ( SC, ET, BORVEC, SBFXPT, FOUND ) */
327 /* IMPLICIT NONE */
328
329 /* CHARACTER*(*) SC */
330 /* DOUBLE PRECISION ET */
331 /* DOUBLE PRECISION BORVEC ( 3 ) */
332 /* DOUBLE PRECISION SBFXPT ( 3 ) */
333 /* LOGICAL FOUND */
334
335 /* C */
336 /* C SPICELIB functions */
337 /* C */
338 /* DOUBLE PRECISION CLIGHT */
339 /* DOUBLE PRECISION VDIST */
340
341 /* C */
342 /* C Local parameters */
343 /* C */
344 /* INTEGER UBPL */
345 /* PARAMETER ( UBPL = 4 ) */
346
347 /* INTEGER SATURN */
348 /* PARAMETER ( SATURN = 699 ) */
349
350 /* C */
351 /* C Local variables */
352 /* C */
353 /* DOUBLE PRECISION BORV2 ( 3 ) */
354 /* DOUBLE PRECISION CORVEC ( 3 ) */
355 /* DOUBLE PRECISION LT */
356 /* DOUBLE PRECISION PLANE ( UBPL ) */
357 /* DOUBLE PRECISION SATSSB ( 6 ) */
358 /* DOUBLE PRECISION SCPOS ( 3 ) */
359 /* DOUBLE PRECISION SCSSB ( 6 ) */
360 /* DOUBLE PRECISION STATE ( 6 ) */
361 /* DOUBLE PRECISION STCORR ( 3 ) */
362 /* DOUBLE PRECISION TAU */
363 /* DOUBLE PRECISION TPMI ( 3, 3 ) */
364 /* DOUBLE PRECISION XPT ( 3 ) */
365 /* DOUBLE PRECISION ZVEC ( 3 ) */
366
367 /* INTEGER I */
368 /* INTEGER NXPTS */
369 /* INTEGER SCID */
370
371 /* LOGICAL FND */
372
373 /* C */
374 /* C First step: account for stellar aberration. Since the */
375 /* C instrument pointing is given, we need to find the intercept */
376 /* C point such that, when the stellar aberration correction is */
377 /* C applied to the vector from the spacecraft to that point, */
378 /* C the resulting vector is parallel to BORVEC. An easy */
379 /* C solution is to apply the inverse of the normal stellar */
380 /* C aberration correction to BORVEC, and then solve the */
381 /* C intercept problem with this corrected boresight vector. */
382 /* C */
383 /* C Find the position of the observer relative */
384 /* C to the solar system barycenter at ET. */
385 /* C */
386 /* CALL BODN2C ( SC, SCID, FND ) */
387
388 /* IF ( .NOT. FND ) THEN */
389
390 /* CALL SETMSG ( 'ID code for body # was not found.' ) */
391 /* CALL ERRCH ( '#', SC ) */
392 /* CALL SIGERR ( 'SPICE(NOTRANSLATION' ) */
393 /* RETURN */
394
395 /* END IF */
396
397 /* CALL SPKSSB ( SCID, ET, 'J2000', SCSSB ) */
398
399 /* C */
400 /* C We now wish to find the vector CORVEC that, when */
401 /* C corrected for stellar aberration, yields BORVEC. */
402 /* C A good first approximation is obtained by applying */
403 /* C the stellar aberration correction for transmission */
404 /* C to BORVEC. */
405 /* C */
406 /* CALL STLABX ( BORVEC, SCSSB(4), CORVEC ) */
407
408 /* C */
409 /* C The inverse of the stellar aberration correction */
410 /* C applicable to CORVEC should be a very good estimate of */
411 /* C the correction we need to apply to BORVEC. Apply */
412 /* C this correction to BORVEC to obtain an improved estimate */
413 /* C of CORVEC. */
414 /* C */
415 /* CALL STELAB ( CORVEC, SCSSB(4), BORV2 ) */
416 /* CALL VSUB ( BORV2, CORVEC, STCORR ) */
417 /* CALL VSUB ( BORVEC, STCORR, CORVEC ) */
418
419 /* C */
420 /* C Because the ring plane intercept may be quite far from */
421 /* C Saturn's center, we cannot assume light time from the */
422 /* C intercept to the observer is well approximated by */
423 /* C light time from Saturn's center to the observer. */
424 /* C We compute the light time explicitly using an iterative */
425 /* C approach. */
426 /* C */
427 /* C We can however use the light time from Saturn's center to */
428 /* C the observer to obtain a first estimate of the actual light */
429 /* C time. */
430 /* C */
431 /* CALL SPKEZR ( 'SATURN', ET, 'J2000', 'LT', SC, */
432 /* . STATE, LT ) */
433 /* TAU = LT */
434
435 /* C */
436 /* C Find the ring plane intercept and calculate the */
437 /* C light time from it to the spacecraft. */
438 /* C Perform three iterations. */
439 /* C */
440 /* I = 1 */
441 /* FOUND = .TRUE. */
442
443 /* DO WHILE ( ( I .LE. 3 ) .AND. ( FOUND ) ) */
444 /* C */
445 /* C Find the position of Saturn relative */
446 /* C to the solar system barycenter at ET-TAU. */
447 /* C */
448 /* CALL SPKSSB ( SATURN, ET-TAU, 'J2000', SATSSB ) */
449
450 /* C */
451 /* C Find the Saturn-to-observer vector defined by these */
452 /* C two position vectors. */
453 /* C */
454 /* CALL VSUB ( SCSSB, SATSSB, SCPOS ) */
455
456 /* C */
457 /* C Look up Saturn's pole at ET-TAU; this is the third */
458 /* C column of the matrix that transforms Saturn body-fixed */
459 /* C coordinates to J2000 coordinates. */
460 /* C */
461 /* CALL PXFORM ( 'IAU_SATURN', 'J2000', ET-TAU, TPMI ) */
462
463 /* CALL MOVED ( TPMI(1,3), 3, ZVEC ) */
464
465 /* C */
466 /* C Make a SPICELIB plane representing the ring plane. */
467 /* C We're treating Saturn's center as the origin, so */
468 /* C the plane constant is 0. */
469 /* C */
470 /* CALL NVC2PL ( ZVEC, 0.D0, PLANE ) */
471
472 /* C */
473 /* C Find the intersection of the ring plane and the */
474 /* C ray having vertex SCPOS and direction vector */
475 /* C CORVEC. */
476 /* C */
477 /* CALL INRYPL ( SCPOS, CORVEC, PLANE, NXPTS, XPT ) */
478
479 /* C */
480 /* C If the number of intersection points is 1, */
481 /* C find the next light time estimate. */
482 /* C */
483 /* IF ( NXPTS .EQ. 1 ) THEN */
484 /* C */
485 /* C Find the light time (zero-order) from the */
486 /* C intercept point to the spacecraft. */
487 /* C */
488 /* TAU = VDIST ( SCPOS, XPT ) / CLIGHT() */
489 /* I = I + 1 */
490
491 /* ELSE */
492
493 /* FOUND = .FALSE. */
494
495 /* END IF */
496
497 /* END DO */
498
499 /* C */
500 /* C At this point, if FOUND is .TRUE., we iterated */
501 /* C 3 times, and XPT is our estimate of the */
502 /* C position of the ring plane intercept point */
503 /* C relative to Saturn in the J2000 frame. This is the */
504 /* C point observed by an instrument pointed in direction */
505 /* C BORVEC at ET at mounted on the spacecraft SC. */
506 /* C */
507 /* C If FOUND is .FALSE., the boresight ray does not */
508 /* C intersect the ring plane. */
509 /* C */
510 /* C As a final step, transform XPT to Saturn body-fixed */
511 /* C coordinates. */
512 /* C */
513 /* IF ( FOUND ) THEN */
514
515 /* CALL MTXV ( TPMI, XPT, SBFXPT ) */
516
517 /* END IF */
518
519 /* END */
520
521
522
523 /* $ Restrictions */
524
525 /* None. */
526
527 /* $ Literature_References */
528
529 /* None. */
530
531 /* $ Author_and_Institution */
532
533 /* N.J. Bachman (JPL) */
534 /* W.L. Taber (JPL) */
535
536 /* $ Version */
537
538 /* - SPICELIB Version 1.2.0, 29-SEP-2016 (NJB) */
539
540 /* Changed from standard to discovery check-in. Fixed typo */
541 /* in header. */
542
543 /* - SPICELIB Version 1.1.1, 07-FEB-2008 (BVS) */
544
545 /* Fixed a few typos in the header. */
546
547 /* - SPICELIB Version 1.1.0, 02-SEP-2005 (NJB) */
548
549 /* Updated to remove non-standard use of duplicate arguments */
550 /* in VSCL call. */
551
552 /* - SPICELIB Version 1.0.3, 12-DEC-2002 (NJB) */
553
554 /* Header fix: ring plane intercept algorithm was corrected. */
555 /* Now light time is computed accurately, and stellar aberration */
556 /* is accounted for. Example was turned into a complete */
557 /* subroutine. */
558
559 /* - SPICELIB Version 1.0.2, 09-MAR-1999 (NJB) */
560
561 /* Reference to SCE2T replaced by reference to SCE2C. An */
562 /* occurrence of ENDIF was replaced by END IF. */
563
564 /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */
565
566 /* Comment section for permuted index source lines was added */
567 /* following the header. */
568
569 /* - SPICELIB Version 1.0.0, 01-APR-1991 (NJB) (WLT) */
570
571 /* -& */
572 /* $ Index_Entries */
573
574 /* intersection of ray and plane */
575
576 /* -& */
577 /* $ Revisions */
578
579 /* - SPICELIB Version 1.1.0, 02-SEP-2005 (NJB) */
580
581 /* Updated to remove non-standard use of duplicate arguments */
582 /* in VSCL call. */
583
584 /* -& */
585
586 /* SPICELIB functions */
587
588
589 /* Local parameters */
590
591
592 /* Local variables */
593
594
595 /* Standard SPICE error handling. */
596
597 if (return_()) {
598 return 0;
599 }
600
601 /* We'll give the name TOOBIG to the bound DPMAX() / MARGIN. */
602 /* If we let VTXPRJ be the orthogonal projection of VERTEX onto */
603 /* PLANE, and let DIFF be the vector VTXPRJ - VERTEX, then */
604 /* we know that */
605
606 /* || DIFF || < 2 * TOOBIG */
607
608 /* Check the distance of the ray's vertex from the origin. */
609
610 toobig = dpmax_() / 3.;
611 if (vnorm_(vertex) >= toobig) {
612 chkin_("INRYPL", (ftnlen)6);
613 setmsg_("Ray's vertex is too far from the origin.", (ftnlen)40);
614 sigerr_("SPICE(VECTORTOOBIG)", (ftnlen)19);
615 chkout_("INRYPL", (ftnlen)6);
616 return 0;
617 }
618
619 /* Check the distance of the plane from the origin. (The returned */
620 /* plane constant IS this distance.) */
621
622 pl2nvc_(plane, normal, &const__);
623 if (const__ >= toobig) {
624 chkin_("INRYPL", (ftnlen)6);
625 setmsg_("Plane is too far from the origin.", (ftnlen)33);
626 sigerr_("SPICE(VECTORTOOBIG)", (ftnlen)19);
627 chkout_("INRYPL", (ftnlen)6);
628 return 0;
629 }
630
631 /* Check the ray's direction vector. */
632
633 vhat_(dir, udir);
634 if (vzero_(udir)) {
635 chkin_("INRYPL", (ftnlen)6);
636 setmsg_("Ray's direction vector is the zero vector.", (ftnlen)42);
637 sigerr_("SPICE(ZEROVECTOR)", (ftnlen)17);
638 chkout_("INRYPL", (ftnlen)6);
639 return 0;
640 }
641
642 /* That takes care of the error cases. Now scale the input vertex */
643 /* and plane to improve numerical behavior. */
644
645 /* Computing MAX */
646 d__1 = const__, d__2 = vnorm_(vertex);
647 mscale = max(d__1,d__2);
648 if (mscale != 0.) {
649 d__1 = 1. / mscale;
650 vscl_(&d__1, vertex, sclvtx);
651 sclcon = const__ / mscale;
652 } else {
653 vequ_(vertex, sclvtx);
654 sclcon = const__;
655 }
656 if (mscale > 1.) {
657 toobig /= mscale;
658 }
659 /* Find the projection (coefficient) of the ray's vertex along the */
660 /* plane's normal direction. */
661
662 prjvn = vdot_(sclvtx, normal);
663
664 /* If this projection is the plane constant, the ray's vertex lies in */
665 /* the plane. We have one intersection or an infinite number of */
666 /* intersections. It all depends on whether the ray actually lies */
667 /* in the plane. */
668
669 /* The absolute value of PRJDIF is the distance of the ray's vertex */
670 /* from the plane. */
671
672 prjdif = sclcon - prjvn;
673 if (prjdif == 0.) {
674
675 /* XPT is the original, unscaled vertex. */
676
677 vequ_(vertex, xpt);
678 if (vdot_(normal, udir) == 0.) {
679
680 /* The ray's in the plane. */
681
682 *nxpts = -1;
683 } else {
684 *nxpts = 1;
685 }
686 return 0;
687 }
688
689 /* Ok, the ray's vertex is not in the plane. The ray may still be */
690 /* parallel to or may point away from the plane. If the ray does */
691 /* point towards the plane, mathematicians would say that the */
692 /* ray does intersect the plane, but the computer may disagree. */
693
694 /* For this routine to find an intersection, both of the following */
695 /* conditions must be met: */
696
697 /* -- The ray must point toward the plane; this happens when */
698 /* PRJDIF has the same sign as < UDIR, NORMAL >. */
699
700 /* -- The vector difference XPT - SCLVTX must not overflow. */
701
702 /* Qualitatively, the case of interest looks something like the */
703 /* picture below: */
704
705
706 /* * SCLVTX */
707 /* |\ */
708 /* | \ <-- UDIR */
709 /* | \ */
710 /* length of this | \| */
711 /* segment is | -* */
712 /* | */
713 /* | PRJDIF | --> | ___________________________ */
714 /* |/ / */
715 /* | * / <-- PLANE */
716 /* /| XPT / */
717 /* / ^ / */
718 /* / | NORMAL / */
719 /* / | . / */
720 /* / |/| / */
721 /* / .---| / / */
722 /* / | |/ / */
723 /* / `---* / */
724 /* / Projection of SCLVTX onto the plane */
725 /* / / */
726 /* / / */
727 /* ---------------------------- */
728
729
730
731
732 /* Find the projection of the direction vector along the plane's */
733 /* normal vector. */
734
735 prjdir = vdot_(udir, normal);
736
737 /* We're done if the ray doesn't point toward the plane. PRJDIF */
738 /* has already been found to be non-zero at this point; PRJDIR is */
739 /* zero if the ray and plane are parallel. The SPICELIB routine */
740 /* SMSGND will return a value of .FALSE. if PRJDIR is zero. */
741
742 if (! smsgnd_(&prjdir, &prjdif)) {
743
744 /* The ray is parallel to or points away from the plane. */
745
746 *nxpts = 0;
747 cleard_(&c__3, xpt);
748 return 0;
749 }
750
751 /* The difference XPT - SCLVTX is the hypotenuse of a right triangle */
752 /* formed by SCLVTX, XPT, and the orthogonal projection of SCLVTX */
753 /* onto the plane. We'll obtain the hypotenuse by scaling UDIR. */
754 /* We must make sure that this hypotenuse does not overflow. The */
755 /* scale factor has magnitude */
756
757 /* | PRJDIF | */
758 /* -------------- */
759 /* | PRJDIR | */
760
761 /* and UDIR is a unit vector, so as long as */
762
763 /* | PRJDIF | < | PRJDIR | * TOOBIG */
764
765 /* the hypotenuse is no longer than TOOBIG. The product can be */
766 /* computed safely since PRJDIR has magnitude 1 or less. */
767
768 if (abs(prjdif) >= abs(prjdir) * toobig) {
769
770 /* If the hypotenuse is too long, we say that no intersection */
771 /* exists. */
772
773 *nxpts = 0;
774 cleard_(&c__3, xpt);
775 return 0;
776 }
777
778 /* We conclude that it's safe to compute XPT. Scale UDIR and add */
779 /* the result to SCLVTX. The addition is safe because both addends */
780 /* have magnitude no larger than TOOBIG. The vector thus obtained */
781 /* is the intersection point. */
782
783 *nxpts = 1;
784 scale = abs(prjdif) / abs(prjdir);
785 vlcom_(&c_b16, sclvtx, &scale, udir, xpt);
786
787 /* Re-scale XPT. This is safe, since TOOBIG has already been */
788 /* scaled to allow for any growth of XPT at this step. */
789
790 vsclip_(&mscale, xpt);
791 return 0;
792 } /* inrypl_ */
793
794