1 /* zzelvupy.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__10000 = 10000;
11 static doublereal c_b79 = 2.;
12 static doublereal c_b90 = .5;
13
14 /* $Procedure ZZELVUPY ( Is ellipse in polygonal field of view? ) */
zzelvupy_(doublereal * ellips,doublereal * vertex,doublereal * axis,integer * n,doublereal * bounds,logical * found)15 /* Subroutine */ int zzelvupy_(doublereal *ellips, doublereal *vertex,
16 doublereal *axis, integer *n, doublereal *bounds, logical *found)
17 {
18 /* Initialized data */
19
20 static doublereal origin[3] = { 0.,0.,0. };
21
22 /* System generated locals */
23 integer bounds_dim2, i__1, i__2, i__3;
24 doublereal d__1, d__2;
25
26 /* Builtin functions */
27 integer s_rnge(char *, integer, char *, integer);
28 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
29 double asin(doublereal), pow_dd(doublereal *, doublereal *);
30
31 /* Local variables */
32 doublereal asep, apex[3];
33 extern /* Subroutine */ int vhat_(doublereal *, doublereal *);
34 extern doublereal vdot_(doublereal *, doublereal *), vsep_(doublereal *,
35 doublereal *);
36 extern /* Subroutine */ int vsub_(doublereal *, doublereal *, doublereal *
37 );
38 static doublereal work[30000] /* was [3][10000] */;
39 doublereal edge1[3], edge2[3], a, b, d__;
40 integer i__, j;
41 doublereal vxpt1[3], vxpt2[3], scale, x, y;
42 extern /* Subroutine */ int chkin_(char *, ftnlen);
43 doublereal plane[4];
44 extern /* Subroutine */ int errdp_(char *, doublereal *, ftnlen), vlcom_(
45 doublereal *, doublereal *, doublereal *, doublereal *,
46 doublereal *);
47 extern doublereal vdist_(doublereal *, doublereal *);
48 doublereal vtemp[3];
49 extern /* Subroutine */ int ucrss_(doublereal *, doublereal *, doublereal
50 *);
51 extern doublereal vnorm_(doublereal *);
52 extern logical vzero_(doublereal *);
53 integer nxpts;
54 extern /* Subroutine */ int el2cgv_(doublereal *, doublereal *,
55 doublereal *, doublereal *), cgv2el_(doublereal *, doublereal *,
56 doublereal *, doublereal *);
57 doublereal hafedg;
58 extern /* Subroutine */ int nvp2pl_(doublereal *, doublereal *,
59 doublereal *);
60 doublereal cp[3];
61 extern /* Subroutine */ int psv2pl_(doublereal *, doublereal *,
62 doublereal *, doublereal *);
63 extern doublereal pi_(void);
64 doublereal hafsec, eplane[4], ellscl[9], center[3], easize, ebsctr[3];
65 extern /* Subroutine */ int saelgv_(doublereal *, doublereal *,
66 doublereal *, doublereal *), inelpl_(doublereal *, doublereal *,
67 integer *, doublereal *, doublereal *);
68 doublereal ctrvec[3], consep, offset[3], pasize, smajor[3];
69 char errmsg[1840];
70 extern /* Subroutine */ int sigerr_(char *, ftnlen), setmsg_(char *,
71 ftnlen);
72 doublereal fovpln[4], vbsctr[3];
73 extern /* Subroutine */ int chkout_(char *, ftnlen);
74 doublereal sminor[3];
75 extern /* Subroutine */ int errint_(char *, integer *, ftnlen), repmot_(
76 char *, char *, integer *, char *, char *, ftnlen, ftnlen, ftnlen,
77 ftnlen);
78 doublereal gv1[3];
79 extern logical return_(void);
80 doublereal gv2[3];
81 extern /* Subroutine */ int inrypl_(doublereal *, doublereal *,
82 doublereal *, integer *, doublereal *);
83 extern integer zzwind_(doublereal *, integer *, doublereal *, doublereal *
84 );
85 doublereal xpt[3], xpt1[3], xpt2[3];
86
87 /* $ Abstract */
88
89 /* Determine whether a specified ellipse intersects the pyramid */
90 /* defined by a polygonal field of view. */
91
92 /* $ Disclaimer */
93
94 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
95 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
96 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
97 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
98 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
99 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
100 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
101 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
102 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
103 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
104
105 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
106 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
107 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
108 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
109 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
110 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
111
112 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
113 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
114 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
115 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
116
117 /* $ Required_Reading */
118
119 /* ELLIPSES */
120 /* PLANES */
121
122 /* $ Keywords */
123
124 /* ELLIPSE */
125 /* GEOMETRY */
126 /* MATH */
127 /* PLANE */
128
129 /* $ Declarations */
130 /* $ Disclaimer */
131
132 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
133 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
134 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
135 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
136 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
137 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
138 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
139 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
140 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
141 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
142
143 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
144 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
145 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
146 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
147 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
148 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
149
150 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
151 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
152 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
153 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
154
155
156 /* Include File: SPICELIB Error Handling Parameters */
157
158 /* errhnd.inc Version 2 18-JUN-1997 (WLT) */
159
160 /* The size of the long error message was */
161 /* reduced from 25*80 to 23*80 so that it */
162 /* will be accepted by the Microsoft Power Station */
163 /* FORTRAN compiler which has an upper bound */
164 /* of 1900 for the length of a character string. */
165
166 /* errhnd.inc Version 1 29-JUL-1997 (NJB) */
167
168
169
170 /* Maximum length of the long error message: */
171
172
173 /* Maximum length of the short error message: */
174
175
176 /* End Include File: SPICELIB Error Handling Parameters */
177
178 /* $ Brief_I/O */
179
180 /* Variable I/O Description */
181 /* -------- --- -------------------------------------------------- */
182 /* ELLIPS I A SPICELIB ellipse. */
183 /* VERTEX I Vertex of a pyramid. */
184 /* AXIS I Axis of a pyramid. */
185 /* N I Number of boundary vectors of the pyramid. */
186 /* BOUNDS I Boundary vectors of the pyramid. */
187 /* FOUND O Flag indicating whether intersection was found. */
188 /* UBEL P Upper bound of SPICELIB ellipse array. */
189 /* UBPL P Upper bound of SPICELIB plane array. */
190 /* MAXFOV P Maximum number of boundary vectors. */
191
192 /* $ Detailed_Input */
193
194 /* ELLIPS is a SPICELIB ellipse having non-zero semi-axes. */
195
196 /* VERTEX is the single point of intersection of the vectors */
197 /* defining the edges of a pyramid. The vectors */
198 /* emanate from this point. The pyramid represents */
199 /* the spatial region viewed by a polygonal field of */
200 /* view (FOV). */
201
202 /* AXIS is a vector emanating from VERTEX that lies inside */
203 /* the pyramid defined by VERTEX, N, and BOUNDS. */
204 /* AXIS represents the boresight direction of the FOV. */
205
206 /* N, */
207 /* BOUNDS are, respectively, the number of boundary vectors */
208 /* defining the pyramid and the boundary vectors */
209 /* themselves. Each pair of consecutive vectors in */
210 /* the array BOUNDS, together with VERTEX, defines a */
211 /* face of the pyramid. */
212
213 /* Each boundary vector must have angular separation */
214 /* of less than pi/2 radians from AXIS. */
215
216 /* For any plane that doesn't contain VERTEX and that */
217 /* intersects AXIS at right angles, the intersections */
218 /* of the boundary vectors with that plane are the */
219 /* vertices of a polygon. The polygon need not be */
220 /* convex, but it must be non-self-intersecting. */
221
222
223 /* $ Detailed_Output */
224
225 /* FOUND is set to .TRUE. if the pyramid and ellipse */
226 /* intersect; otherwise FOUND is .FALSE. */
227
228 /* $ Parameters */
229
230 /* UBEL is the array upper bound for SPICELIB ellipses. */
231
232 /* UBPL is the array upper bound for SPICELIB planes. */
233
234 /* MAXFOV is the maximum number of boundary vectors that */
235 /* may be supplied in the input array argument */
236 /* BOUNDS. */
237
238 /* $ Exceptions */
239
240 /* If an error is found, the output argument FOUND will be set to */
241 /* .FALSE. */
242
243
244 /* 1) If either of the semi-axes of the input ellipse is the */
245 /* zero vector, the error SPICE(ZEROVECTOR) will be signaled. */
246
247 /* 2) If the norm of the input ellipse's semi-minor axis is */
248 /* zero after division by the maximum of the norms of the */
249 /* semi-major axis, the ellipse's center, and the vertex of */
250 /* the pyramid, the error SPICE(DEGENERATECASE) will be */
251 /* signaled. */
252
253 /* 3) If the vertex of the pyramid lies in the plane containing */
254 /* the ellipse, at most the edge of the ellipse can be "seen" */
255 /* from the vertex. This case is not considered to be an */
256 /* error. */
257
258 /* 4) If the number of boundary vectors N is not at least 3, */
259 /* or if the number exceeds MAXFOV, the error */
260 /* SPICE(INVALIDCOUNT) will be signaled. */
261
262 /* 5) If any boundary vector is the zero vector, the error */
263 /* SPICE(ZEROVECTOR) will be signaled. */
264
265 /* 6) If the axis is the zero vector, the error SPICE(ZEROVECTOR) */
266 /* will be signaled. */
267
268 /* 7) If any boundary vector has angular separation of at least */
269 /* pi/2 radians from AXIS, the error SPICE(INVALIDFOV) */
270 /* will be signaled. */
271
272 /* 8) If any boundary vector has angular separation of zero */
273 /* radians from one of its neighbors, the error SPICE(INVALIDFOV) */
274 /* will be signaled. */
275
276 /* 9) No test is done to ensure that the input boundary vectors */
277 /* define a non-self-intersecting polygon via their intersection */
278 /* with a plane normal to AXIS. If the boundary vectors don't */
279 /* meet this condition, the results of this routine are */
280 /* unreliable. */
281
282 /* 10) The pyramidal field of view and the input ellipse must not */
283 /* differ too radically in scale, or great loss of precision */
284 /* will result, making the results of this routine unreliable. */
285 /* For example, if the ratio of the norm of the semi-minor axis */
286 /* of the ellipse to the distance from VERTEX to the center of */
287 /* the ellipse is less than double precision epsilon on the host */
288 /* system, a meaningful result can't be computed. */
289
290 /* This routine does not attempt to judge the minimum */
291 /* acceptable level of accuracy. */
292
293 /* $ Files */
294
295 /* None. */
296
297 /* $ Particulars */
298
299 /* This routine is useful for determining whether an ellipsoidal */
300 /* body is in the field of view of a remote-sensing instrument */
301 /* with a field of view having polygonal cross section. */
302
303 /* $ Examples */
304
305 /* Test an ellipse for intersection with a square field */
306 /* of view. */
307
308
309 /* PROGRAM EX1 */
310 /* IMPLICIT NONE */
311
312 /* INTEGER MAXN */
313 /* PARAMETER ( MAXN = 4 ) */
314
315 /* INTEGER UBEL */
316 /* PARAMETER ( UBEL = 9 ) */
317
318 /* DOUBLE PRECISION AXIS ( 3 ) */
319 /* DOUBLE PRECISION CENTER ( 3 ) */
320 /* DOUBLE PRECISION ELLIPS ( UBEL ) */
321 /* DOUBLE PRECISION FOV ( 3, MAXN ) */
322 /* DOUBLE PRECISION SMAJOR ( 3 ) */
323 /* DOUBLE PRECISION SMINOR ( 3 ) */
324 /* DOUBLE PRECISION VERTEX ( 3 ) */
325
326 /* INTEGER N */
327
328 /* LOGICAL FOUND */
329
330 /* C */
331 /* C The FOV (field of view) "looks" in the -x direction: */
332 /* C the axis of the FOV is parallel to the x axis. */
333 /* C The FOV intersects the plane of the ellipse in a */
334 /* C square having height and width 4 units. The edges */
335 /* C of the square are parallel to the y and z axes. */
336 /* C */
337 /* N = 4 */
338
339 /* CALL VPACK ( -1.D0, -1.D0, -1.D0, FOV(1,1) ) */
340 /* CALL VPACK ( -1.D0, 1.D0, -1.D0, FOV(1,2) ) */
341 /* CALL VPACK ( -1.D0, 1.D0, 1.D0, FOV(1,3) ) */
342 /* CALL VPACK ( -1.D0, -1.D0, 1.D0, FOV(1,4) ) */
343
344 /* CALL VPACK ( -1.D0, 0.D0, 0.D0, AXIS ) */
345 /* CALL VPACK ( 1.D0, 0.D0, 0.D0, VERTEX ) */
346
347 /* C */
348 /* C The ellipse is oriented with the major axis */
349 /* C vertical and is parallel to the x-z plane. The ellipse */
350 /* C lies in the plane defined by x = -1. The ellipse */
351 /* C ever-so-slightly overlaps the bottom edge of the FOV. */
352 /* C */
353 /* CALL VPACK ( 0.D0, 0.D0, 1.D0, SMAJOR ) */
354 /* CALL VPACK ( 0.D0, 5.D-1, 0.D0, SMINOR ) */
355 /* CALL VPACK ( -1.D0, 0.D0, -3.D0 + 1.D-12, CENTER ) */
356
357 /* C */
358 /* C Create a SPICELIB ellipse from the center and semi-axes. */
359 /* C */
360 /* CALL CGV2EL ( CENTER, SMAJOR, SMINOR, ELLIPS ) */
361
362 /* C */
363 /* C Test for intersection. We expect an intersection to be */
364 /* C found. */
365 /* C */
366 /* CALL ZZELVUPY ( ELLIPS, VERTEX, AXIS, N, FOV, FOUND ) */
367
368 /* WRITE (*,*) 'Case 1: FOUND = ', FOUND */
369
370 /* C */
371 /* C Shift the ellipse center to move the ellipse outside of */
372 /* C the FOV, then repeat the test. We expect FOUND to be */
373 /* C .FALSE. */
374 /* C */
375 /* CALL VPACK ( -1.D0, 0.D0, -3.D0 - 1.D-12, CENTER ) */
376
377 /* CALL CGV2EL ( CENTER, SMAJOR, SMINOR, ELLIPS ) */
378
379 /* CALL ZZELVUPY ( ELLIPS, VERTEX, AXIS, N, FOV, FOUND ) */
380
381 /* WRITE (*,*) 'Case 2: FOUND = ', FOUND */
382
383 /* END */
384
385
386 /* $ Restrictions */
387
388 /* None. */
389
390 /* $ Literature_References */
391
392 /* [1] `Calculus and Analytic Geometry', Thomas and Finney. */
393
394 /* $ Author_and_Institution */
395
396 /* N.J. Bachman (JPL) */
397 /* B.V. Semenov (JPL) */
398
399 /* $ Version */
400
401 /* - SPICELIB Version 1.0.1, 28-FEB-2008 (BVS) */
402
403 /* Corrected the contents of the Required_Reading section. */
404
405 /* - SPICELIB Version 1.0.0, 10-AUG-2005 (NJB) */
406
407 /* -& */
408 /* $ Index_Entries */
409
410 /* test whether pyramid intersects ellipse */
411 /* test whether ellipse is in pyramidal field of view */
412
413 /* -& */
414
415 /* SPICELIB functions */
416
417
418 /* Local parameters */
419
420
421 /* Local variables */
422
423
424 /* Saved variables */
425
426
427 /* Initial values */
428
429 /* Parameter adjustments */
430 bounds_dim2 = *n;
431
432 /* Function Body */
433 if (return_()) {
434 return 0;
435 }
436 chkin_("ZZELVUPY", (ftnlen)8);
437
438 /* We start out by checking the inputs. */
439
440 /* The next step will be to look for an intersection of the ellipse */
441 /* and pyramid. There are three intersection cases: */
442
443 /* 1) The ellipse is completely contained in the pyramid. */
444
445 /* 2) The ellipse "contains" the field of view in the sense */
446 /* that the intersection of the pyramid and the plane of the */
447 /* ellipse is contained in the region bounded by the ellipse. */
448
449 /* 3) One or more sides of the pyramid intersect the ellipse. */
450
451 /* There is also a non-intersection case: this is when cones */
452 /* bounding the ellipse and pyramid and having their apexes in */
453 /* common with that of the pyramid intersect only in that common */
454 /* apex. Before test (1), we perform this non-intersection test, */
455 /* since it can be done quickly. */
456
457 /* No intersection has been found so far. Set the default value */
458 /* of the FOUND flag here so it won't have to be set in every error */
459 /* checking block below. */
460
461 *found = FALSE_;
462
463 /* Validate the ellipse. First find the center and the semi-axes */
464 /* of the ellipse. */
465
466 el2cgv_(ellips, center, gv1, gv2);
467 saelgv_(gv1, gv2, smajor, sminor);
468
469 /* Check the semi-axis lengths. */
470
471 /* If the semi-major axis is the zero vector, we'd expect */
472 /* the semi-minor axis to be the zero vector as well. But */
473 /* round-off error could conceivably violate this assumption. */
474
475 if (vzero_(smajor) || vzero_(sminor)) {
476 setmsg_("Input ellipse has semi-major axis length # and semi-minor a"
477 "xis length #. Both vectors are required to be non-zero.", (
478 ftnlen)115);
479 d__1 = vnorm_(smajor);
480 errdp_("#", &d__1, (ftnlen)1);
481 d__1 = vnorm_(sminor);
482 errdp_("#", &d__1, (ftnlen)1);
483 sigerr_("SPICE(ZEROVECTOR)", (ftnlen)17);
484 chkout_("ZZELVUPY", (ftnlen)8);
485 return 0;
486 }
487
488 /* Scale the vectors defining the ellipse and the vertex of the */
489 /* pyramid so that the largest of these vectors has unit length. */
490
491 /* Computing MAX */
492 d__1 = vnorm_(center), d__2 = vnorm_(smajor), d__1 = max(d__1,d__2), d__2
493 = vnorm_(vertex);
494 scale = 1. / max(d__1,d__2);
495 for (i__ = 1; i__ <= 3; ++i__) {
496 center[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("center",
497 i__1, "zzelvupy_", (ftnlen)452)] = scale * center[(i__2 = i__
498 - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("center", i__2, "zzelv"
499 "upy_", (ftnlen)452)];
500 smajor[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("smajor",
501 i__1, "zzelvupy_", (ftnlen)453)] = scale * smajor[(i__2 = i__
502 - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("smajor", i__2, "zzelv"
503 "upy_", (ftnlen)453)];
504 sminor[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("sminor",
505 i__1, "zzelvupy_", (ftnlen)454)] = scale * sminor[(i__2 = i__
506 - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("sminor", i__2, "zzelv"
507 "upy_", (ftnlen)454)];
508 apex[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("apex", i__1,
509 "zzelvupy_", (ftnlen)455)] = scale * vertex[(i__2 = i__ - 1) <
510 3 && 0 <= i__2 ? i__2 : s_rnge("vertex", i__2, "zzelvupy_", (
511 ftnlen)455)];
512 }
513
514 /* Create a scaled ellipse. We'll perform the FOV side-ellipse */
515 /* intersection computations using this ellipse. */
516
517 cgv2el_(center, smajor, sminor, ellscl);
518
519 /* After scaling, make sure the semi-axes have sufficient length to */
520 /* prevent numerical problems. Let A and B be the scaled semi-axis */
521 /* lengths of the ellipse. */
522
523 a = vnorm_(smajor);
524 b = vnorm_(sminor);
525 if (b == 0.) {
526 setmsg_("Scaled ellipse's semi-minor axis length = 0.", (ftnlen)44);
527 sigerr_("SPICE(DEGENERATECASE)", (ftnlen)21);
528 chkout_("ZZELVUPY", (ftnlen)8);
529 return 0;
530 }
531
532 /* Validate the input pyramid. */
533
534 /* The axis must not be the zero vector. */
535
536 if (vzero_(axis)) {
537 setmsg_("The pyramid's axis the zero vector.", (ftnlen)35);
538 sigerr_("SPICE(ZEROVECTOR)", (ftnlen)17);
539 chkout_("ZZELVUPY", (ftnlen)8);
540 return 0;
541 }
542
543 /* There must be at least three boundary vectors. */
544
545 if (*n < 3) {
546 setmsg_("The number of boundary vectors was #; this number must be a"
547 "t least 3.", (ftnlen)69);
548 errint_("#", n, (ftnlen)1);
549 sigerr_("SPICE(INVALIDCOUNT)", (ftnlen)19);
550 chkout_("ZZELVUPY", (ftnlen)8);
551 return 0;
552 }
553
554 /* There must be no more than MAXFOV boundary vectors. */
555
556 if (*n > 10000) {
557 setmsg_("The number of boundary vectors was #; this number must not "
558 "exceed #.", (ftnlen)68);
559 errint_("#", n, (ftnlen)1);
560 errint_("#", &c__10000, (ftnlen)1);
561 sigerr_("SPICE(INVALIDCOUNT)", (ftnlen)19);
562 chkout_("ZZELVUPY", (ftnlen)8);
563 return 0;
564 }
565
566 /* We must initialize certain variables before continuing with */
567 /* the checks. */
568
569 /* Let CTRVEC be the vector from the apex to the center of the */
570 /* ellipse. This vector will be used in several places later; */
571 /* it's convenient to compute it here. */
572
573 vsub_(center, apex, ctrvec);
574
575 /* Compute PASIZE: an upper bound on the angular radius of a */
576 /* circular cone whose axis is the input central axis. While */
577 /* we're at it, check the angular separation of the boundary */
578 /* vectors from the central axis and from each other. */
579
580 pasize = 0.;
581 i__1 = *n;
582 for (i__ = 1; i__ <= i__1; ++i__) {
583
584 /* Each boundary vector must have angular separation from the */
585 /* axis of less than pi/2 radians. Keep track of the maximum */
586 /* angular separation PASIZE as we go. We'll use this variable */
587 /* later in a non-intersection test. */
588
589 asep = vsep_(axis, &bounds[(i__2 = i__ * 3 - 3) < bounds_dim2 * 3 &&
590 0 <= i__2 ? i__2 : s_rnge("bounds", i__2, "zzelvupy_", (
591 ftnlen)550)]);
592 if (asep >= pi_() / 2) {
593 setmsg_("The angular separation of boundary vector # from the ax"
594 "is is #. This number must less than pi/2.", (ftnlen)96);
595 errint_("#", &i__, (ftnlen)1);
596 errdp_("#", &asep, (ftnlen)1);
597 sigerr_("SPICE(INVALIDFOV)", (ftnlen)17);
598 chkout_("ZZELVUPY", (ftnlen)8);
599 return 0;
600 }
601 pasize = max(pasize,asep);
602
603 /* Each boundary vector must have non-zero angular separation */
604 /* from its neighbors. */
605
606 if (i__ < *n) {
607 j = i__ + 1;
608 } else {
609 j = 1;
610 }
611 ucrss_(&bounds[(i__2 = i__ * 3 - 3) < bounds_dim2 * 3 && 0 <= i__2 ?
612 i__2 : s_rnge("bounds", i__2, "zzelvupy_", (ftnlen)577)], &
613 bounds[(i__3 = j * 3 - 3) < bounds_dim2 * 3 && 0 <= i__3 ?
614 i__3 : s_rnge("bounds", i__3, "zzelvupy_", (ftnlen)577)], cp);
615 if (vzero_(cp)) {
616
617 /* The cross product may be zero because one of the */
618 /* boundary vectors is zero. Check this first. */
619
620 if (vzero_(&bounds[(i__2 = j * 3 - 3) < bounds_dim2 * 3 && 0 <=
621 i__2 ? i__2 : s_rnge("bounds", i__2, "zzelvupy_", (ftnlen)
622 584)]) || vzero_(&bounds[(i__3 = i__ * 3 - 3) <
623 bounds_dim2 * 3 && 0 <= i__3 ? i__3 : s_rnge("bounds",
624 i__3, "zzelvupy_", (ftnlen)584)])) {
625 s_copy(errmsg, "The # boundary vector is the zero vector.", (
626 ftnlen)1840, (ftnlen)41);
627 if (vzero_(&bounds[(i__2 = i__ * 3 - 3) < bounds_dim2 * 3 &&
628 0 <= i__2 ? i__2 : s_rnge("bounds", i__2, "zzelvupy_",
629 (ftnlen)588)])) {
630 j = i__;
631 }
632 repmot_(errmsg, "#", &j, "L", errmsg, (ftnlen)1840, (ftnlen)1,
633 (ftnlen)1, (ftnlen)1840);
634 setmsg_(errmsg, (ftnlen)1840);
635 sigerr_("SPICE(ZEROVECTOR)", (ftnlen)17);
636 } else {
637 setmsg_("The angular separation of boundary vector # from ve"
638 "ctor # is 0.This number must be positive.", (ftnlen)
639 92);
640 errint_("#", &i__, (ftnlen)1);
641 errint_("#", &j, (ftnlen)1);
642 sigerr_("SPICE(INVALIDFOV)", (ftnlen)17);
643 }
644 chkout_("ZZELVUPY", (ftnlen)8);
645 return 0;
646 }
647 }
648
649 /* That's it for the error checks. We'll now answer the question */
650 /* this routine is meant to answer: does the ellipse or the region */
651 /* it bounds intersect the pyramid? */
652
653 /* We'll start out with a simple check to rule out intersection */
654 /* when the ellipse and pyramid are contained in disjoint right */
655 /* circular cones with a common apex. */
656
657 /* Find the angular radius (that is, one-half of the angular extent) */
658 /* of a bounding cone of the ellipse as seen from the apex. The */
659 /* cone circumscribes a sphere of radius A centered at the ellipse's */
660 /* center, where A is the length of the semi-major axis. Note that */
661 /* the cone does not in general circumscribe the ellipse itself. */
662
663 /* The test can be performed only if the apex of the FOV is outside */
664 /* of the sphere of radius A centered at the ellipse center. */
665
666 d__ = vdist_(center, apex);
667 if (a < d__) {
668 easize = asin(a / d__);
669
670 /* The variable PASIZE already contains the angular radius of a */
671 /* bounding cone of the pyramid as seen from the pyramid's apex. */
672 /* The angular radius is the maximum of the angular separations */
673 /* of each pyramid edge from the pyramid's axis. Check whether */
674 /* the bounding cones of ellipse and pyramid are disjoint. Recall */
675 /* CTRVEC is the vector from the apex to the center of the */
676 /* ellipse. If the angular separation of CTRVEC and AXIS exceeds */
677 /* the sum of the angular radii of the ellipse's and pyramid's */
678 /* bounding cones, there can be no intersection. */
679
680 consep = vsep_(ctrvec, axis) - (easize + pasize);
681 if (consep > 0.) {
682 chkout_("ZZELVUPY", (ftnlen)8);
683 return 0;
684 }
685 }
686
687 /* At this point, we have to take a more detailed look at the */
688 /* possible intersection of ellipse and pyramid. First check */
689 /* whether the center of the ellipse is contained in the pyramid. */
690 /* If the ellipse is completely contained in the pyramid, this */
691 /* check will yield a positive result. */
692
693 /* The center of the ellipse is inside the pyramid if a plane */
694 /* containing this point and normal to the axis vector */
695 /* chops the pyramid in a polygon that has non-zero winding */
696 /* number about the center. */
697
698 /* The center of the ellipse must lie in the correct half-space */
699 /* for this test to be applicable. */
700
701 if (vdot_(axis, ctrvec) > 0.) {
702
703 /* Construct the plane and find the polygon. */
704
705 nvp2pl_(axis, ctrvec, fovpln);
706
707 /* Create the planar FOV boundary using the intersections */
708 /* of the FOV boundary vectors with FOVPLN. */
709
710 i__1 = *n;
711 for (i__ = 1; i__ <= i__1; ++i__) {
712 inrypl_(origin, &bounds[(i__2 = i__ * 3 - 3) < bounds_dim2 * 3 &&
713 0 <= i__2 ? i__2 : s_rnge("bounds", i__2, "zzelvupy_", (
714 ftnlen)686)], fovpln, &nxpts, &work[(i__3 = i__ * 3 - 3) <
715 30000 && 0 <= i__3 ? i__3 : s_rnge("work", i__3, "zzelv"
716 "upy_", (ftnlen)686)]);
717
718 /* We expect to have a single point of intersection for each */
719 /* boundary vector. */
720
721 if (nxpts != 1) {
722 setmsg_("NXPTS = # for boundary vector #/FOV plane intersect"
723 "ion.", (ftnlen)55);
724 errint_("#", &nxpts, (ftnlen)1);
725 errint_("#", &i__, (ftnlen)1);
726 sigerr_("SPICE(BUG)", (ftnlen)10);
727 chkout_("ZZELVUPY", (ftnlen)8);
728 return 0;
729 }
730 }
731
732 /* Now WORK contains the polygon representing the intersection of */
733 /* the pyramid with the plane FOVPLN. If the winding number of */
734 /* the polygon about the ellipse center is non-zero, we conclude */
735 /* the center is in the pyramid. */
736
737 if (zzwind_(fovpln, n, work, ctrvec) != 0) {
738
739 /* The center of the ellipse is inside the pyramid. We're */
740 /* done. */
741
742 *found = TRUE_;
743 chkout_("ZZELVUPY", (ftnlen)8);
744 return 0;
745 }
746 }
747
748 /* Check whether the ray defined by APEX and the first boundary */
749 /* vector of the pyramid (the "boundary ray") intersects the plane */
750 /* region bounded by the ellipse. If the intersection of the */
751 /* pyramid and the plane of the ellipse is completely contained in */
752 /* the region bounded by the ellipse, this check will yield a */
753 /* positive result. */
754
755 /* First find the intersection of the boundary ray and the plane */
756 /* containing the ellipse; represent this plane using the SPICELIB */
757 /* plane EPLANE. */
758 /* We don't check FAILED() here because the spanning vectors */
759 /* are orthogonal, and because PSV2PL (via a call to UCRSS) */
760 /* does scaling to prevent underflow. */
761
762 psv2pl_(center, smajor, sminor, eplane);
763 inrypl_(apex, &bounds[(i__1 = 0) < bounds_dim2 * 3 ? i__1 : s_rnge("boun"
764 "ds", i__1, "zzelvupy_", (ftnlen)745)], eplane, &nxpts, xpt);
765
766 /* The routine INRYPL can return the NXPTS values 1, 0, or INF---a */
767 /* code indicating an infinite number of intersection points of ray */
768 /* and plane. If the value is 1, the boundary ray may intersect */
769 /* the region bounded by the ellipse. */
770
771 if (nxpts == 1) {
772
773 /* The boundary ray intersects the plane of the ellipse in a */
774 /* single point. Decide whether this point is inside the ellipse. */
775 /* To test for containment, find the "coordinates" of the */
776 /* center-to-point vector relative to the two-dimensional basis */
777 /* formed by the semi-axes of the ellipse. Call this */
778 /* center-to-point vector OFFSET. Recall A and B are the */
779 /* semi-axis lengths of the ellipse. Let X and Y be the */
780 /* coordinates of OFFSET in the two-dimensional reference frame */
781 /* whose basis consists of normalized versions of SMAJOR and */
782 /* SMINOR. */
783
784 /* Note that we could have the special case in which the vertex */
785 /* of the pyramid lies in the plane of the ellipse, in which case */
786 /* the FOV "sees" the ellipse edge-on. However, since NXPTS is */
787 /* not INF, the boundary vector does not lie in the plane of the */
788 /* ellipse. So in this special case, APEX would be in the region */
789 /* bounded by the ellipse. */
790
791 vsub_(xpt, center, offset);
792 x = vdot_(offset, smajor) / a;
793 y = vdot_(offset, sminor) / b;
794 d__1 = x / a;
795 d__2 = y / b;
796 if (pow_dd(&d__1, &c_b79) + pow_dd(&d__2, &c_b79) <= 1.) {
797
798 /* The boundary-vector-plane intercept lies in the */
799 /* topologically closed region bounded by the ellipse. */
800
801 *found = TRUE_;
802 chkout_("ZZELVUPY", (ftnlen)8);
803 return 0;
804 }
805 }
806
807 /* Check whether one of the pyramid's sides intersects the ellipse. */
808 /* For each side, we first test whether the plane containing that */
809 /* side intersects the ellipse. If it does, the intersection is */
810 /* a (possibly degenerate) line segment with endpoints on the */
811 /* ellipse. The triangle (or segment) defined by the pyramid's */
812 /* apex and this segment (point) is then checked for intersection */
813 /* with the currently considered side of the pyramid. */
814
815 i__ = 1;
816 while(i__ <= *n && ! (*found)) {
817
818 /* Create a SPICELIB plane containing the Ith side of the */
819 /* pyramid. */
820
821 if (i__ < *n) {
822 j = i__ + 1;
823 } else {
824 j = 1;
825 }
826
827 /* Although PSV2PL can signal an error if the spanning */
828 /* vectors are linearly dependent, it won't do so here */
829 /* because we've already ensured the cross product of */
830 /* these vectors is non-zero. */
831
832 psv2pl_(apex, &bounds[(i__1 = i__ * 3 - 3) < bounds_dim2 * 3 && 0 <=
833 i__1 ? i__1 : s_rnge("bounds", i__1, "zzelvupy_", (ftnlen)820)
834 ], &bounds[(i__2 = j * 3 - 3) < bounds_dim2 * 3 && 0 <= i__2 ?
835 i__2 : s_rnge("bounds", i__2, "zzelvupy_", (ftnlen)820)],
836 plane);
837
838 /* Find the intersection of the plane and the ellipse. */
839
840 inelpl_(ellscl, plane, &nxpts, xpt1, xpt2);
841
842 /* If the ellipse-plane intersection is non-empty, test it to see */
843 /* whether it has non-empty intersection with the current side of */
844 /* the pyramid. */
845
846 if (nxpts > 0) {
847
848 /* Let EDGE1 and EDGE2 be the unit length boundary vectors */
849 /* forming the edges of the currently considered side of the */
850 /* pyramid. */
851
852 vhat_(&bounds[(i__1 = i__ * 3 - 3) < bounds_dim2 * 3 && 0 <= i__1
853 ? i__1 : s_rnge("bounds", i__1, "zzelvupy_", (ftnlen)837)]
854 , edge1);
855 vhat_(&bounds[(i__1 = j * 3 - 3) < bounds_dim2 * 3 && 0 <= i__1 ?
856 i__1 : s_rnge("bounds", i__1, "zzelvupy_", (ftnlen)838)],
857 edge2);
858
859 /* Let EBSCTR ("pyramid edge bisector") be a bisector of the */
860 /* sector bounded by EDGE1 and EDGE2. */
861
862 vlcom_(&c_b90, edge1, &c_b90, edge2, ebsctr);
863
864 /* Let HAFEDG be half of the angular measure of this sector. */
865
866 hafedg = vsep_(edge1, edge2) / 2.;
867
868 /* Let VXPT1 and VXPT2 be the unit vectors pointing from the */
869 /* pyramid's apex to the points of intersection of the ellipse */
870 /* and the plane containing the currently considered side of */
871 /* the pyramid. */
872
873 vsub_(xpt1, apex, vtemp);
874 vhat_(vtemp, vxpt1);
875 vsub_(xpt2, apex, vtemp);
876 vhat_(vtemp, vxpt2);
877
878 /* At this point we'll introduce a bit of terminology. We're */
879 /* going to work with plane regions defined by pairs of */
880 /* vectors with a common endpoint. We'll abuse standard */
881 /* terminology a bit and call the region bounded by such a */
882 /* vector pair a "sector." Strictly speaking, sectors refer */
883 /* only to subsets of a disc. */
884
885 /* When it's convenient, we'll also identify "sectors" with */
886 /* regions of the unit circle. This will make it possible */
887 /* to talk about intersections of sectors in terms of */
888 /* intersections of the associated arcs on the unit circle. */
889 /* By the "endpoints" of a sector we mean the endpoints */
890 /* of the arc associated with the sector on the unit circle. */
891
892 /* Let VBSCTR ("VXPT bisector") be a bisector of the sector */
893 /* bounded by VXPT1 and VXPT2. */
894
895 vlcom_(&c_b90, vxpt1, &c_b90, vxpt2, vbsctr);
896
897 /* Let HAFSEC be half of the angular measure of the sector */
898 /* bounded by VXPT1 and VXPT2. */
899
900 hafsec = vsep_(vxpt1, vxpt2) / 2.;
901
902 /* EDGE1, EDGE2, VXPT1, and VXPT2 are four co-planar vectors */
903 /* emanating from APEX. We want to find out whether the */
904 /* sector bounded by EDGE1 and EDGE2 intersects the sector */
905 /* bounded by VXPT1 and VXPT2. If there's an intersection, at */
906 /* least one endpoint of one sector is contained in the other */
907 /* sector. */
908
909 /* Because of potential round-off problems when the sectors */
910 /* are nearly coincident, we perform the precautionary check */
911 /* (case 3) on the angle bisector of the sector defined by */
912 /* VXPT1 and VXPT2. */
913
914 /* If the sector defined by VXPT1 and VXPT2 has no endpoint */
915 /* contained in the other sector, it's possible that the */
916 /* former sector contains the latter. In that case the */
917 /* angular bisector of the latter sector is contained in the */
918 /* former (case 4). */
919
920 /* We test a vector's containment in a sector by comparing the */
921 /* vector's angular separation from the sector's angle */
922 /* bisector to one-half of the angular measure of the sector. */
923
924 /* Case 1: VXPT1 lies between EDGE1 and EDGE2. */
925 /* Case 2: VXPT2 lies between EDGE1 and EDGE2. */
926 /* Case 3: VBSCTR lies between EDGE1 and EDGE2. */
927 /* Case 4: EBSCTR lies between VXPT1 and VXPT2. */
928
929 if (vsep_(vxpt1, ebsctr) <= hafedg) {
930 *found = TRUE_;
931 } else if (vsep_(vxpt2, ebsctr) <= hafedg) {
932 *found = TRUE_;
933 } else if (vsep_(vbsctr, ebsctr) <= hafedg) {
934 *found = TRUE_;
935 } else if (vsep_(ebsctr, vbsctr) <= hafsec) {
936 *found = TRUE_;
937 }
938 if (*found) {
939
940 /* We've found an intersection. We're done. */
941
942 chkout_("ZZELVUPY", (ftnlen)8);
943 return 0;
944 }
945 }
946
947 /* If no intersection was found, look at the next side of the */
948 /* pyramid. */
949
950 ++i__;
951 }
952
953 /* If we got this far, the ellipse is not in view. FOUND has */
954 /* already been set to .FALSE. */
955
956 chkout_("ZZELVUPY", (ftnlen)8);
957 return 0;
958 } /* zzelvupy_ */
959
960