1 /* xdda.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__1 = 1;
11 static doublereal c_b67 = 0.;
12 static doublereal c_b68 = 1.;
13
14 /* $Procedure XDDA ( list voxels intersected by a ray ) */
xdda_(doublereal * vertex,doublereal * raydir,integer * grdext,integer * maxnvx,integer * nvx,integer * voxlst)15 /* Subroutine */ int xdda_(doublereal *vertex, doublereal *raydir, integer *
16 grdext, integer *maxnvx, integer *nvx, integer *voxlst)
17 {
18 /* Initialized data */
19
20 static integer next[3] = { 2,3,1 };
21
22 /* System generated locals */
23 integer i__1, i__2, i__3, i__4;
24 doublereal d__1;
25
26 /* Builtin functions */
27 integer s_rnge(char *, integer, char *, integer);
28
29 /* Local variables */
30 integer step[3], i__;
31 extern /* Subroutine */ int chkin_(char *, ftnlen);
32 extern doublereal dpmax_(void);
33 extern /* Subroutine */ int errdp_(char *, doublereal *, ftnlen);
34 integer iaxis[3];
35 doublereal limit;
36 extern logical vzero_(doublereal *);
37 doublereal ax2err, ax3err, s12, s13;
38 extern doublereal brcktd_(doublereal *, doublereal *, doublereal *);
39 extern integer brckti_(integer *, integer *, integer *);
40 integer icoord[3];
41 doublereal maxcmp;
42 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
43 ftnlen), setmsg_(char *, ftnlen), errint_(char *, integer *,
44 ftnlen);
45 doublereal vtxoff[3];
46 extern logical return_(void);
47 integer intvtx[3];
48 extern logical zzingrd_(integer *, integer *);
49
50 /* $ Abstract */
51
52 /* Given a ray and a voxel grid, return a list of voxels the ray */
53 /* intersects. */
54
55 /* $ Disclaimer */
56
57 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
58 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
59 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
60 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
61 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
62 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
63 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
64 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
65 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
66 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
67
68 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
69 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
70 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
71 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
72 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
73 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
74
75 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
76 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
77 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
78 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
79
80 /* $ Required_Reading */
81
82 /* None. */
83
84 /* $ Keywords */
85
86 /* GRID */
87 /* INTERSECTION */
88 /* PLATE */
89 /* VOXEL */
90
91 /* $ Declarations */
92 /* $ Brief_I/O */
93
94 /* VARIABLE I/O DESCRIPTION */
95 /* -------- --- -------------------------------------------------- */
96 /* GRDTOL P Tolerance for vertex distance from grid. */
97 /* VERTEX I Voxel grid coordinates of ray's vertex. */
98 /* RAYDIR I Direction vector of ray. */
99 /* GRDEXT I Dimensions of grid in voxel units. */
100 /* MAXNVX I Maximum value of VOXLST. */
101 /* NVX O Number of voxels in the VOXLST list. */
102 /* VOXLST O List of voxels intersected by ray. */
103
104 /* $ Detailed_Input */
105
106 /* VERTEX Voxel grid coordinates of ray's vertex. These */
107 /* coordinates are zero-based, double precision */
108 /* offsets from the grid's origin. The units of */
109 /* the coordinates are voxels, that is, */
110 /* voxel edge lengths. */
111
112 /* RAYDIR Direction vector of ray from VERTEX. */
113
114 /* GRDEXT The integer 3-vector containing the voxel grid */
115 /* extents. These are the dimensions of the voxel grid in */
116 /* voxel units, in the X, Y, and Z directions */
117 /* respectively. */
118
119 /* MAXNVX The maximum number of voxel coordinate sets that */
120 /* can be stored in VOXLST. */
121
122 /* $ Detailed_Output */
123
124 /* NVX Number of voxel coordinate sets contained in VOXLST. */
125
126 /* VOXLST List of coordinate sets of voxels intersected by ray. */
127 /* Elements */
128
129 /* VOXLST(J,I), J = 1, 3 */
130
131 /* are the coordinates of the Ith voxel in the list. */
132 /* These coordinates are 1-based integer values. */
133
134 /* The voxels in the output list are ordered by */
135 /* increasing distance from the ray's vertex. */
136
137 /* $ Parameters */
138
139 /* GRDTOL is a tolerance value used to determine whether */
140 /* VERTEX is too far from the voxel grid. The Ith */
141 /* component of VERTEX must not differ from the */
142 /* Ith coordinate of the nearest grid point by more */
143 /* than */
144
145 /* GRDTOL * EXTENT(I) */
146
147 /* $ Exceptions */
148
149 /* 1) The error SPICE(ZEROVECTOR) is signaled if the input RAYDIR */
150 /* has all zero components. */
151
152 /* 2) The error SPICE(INVALIDSIZE) is signaled if the maximum */
153 /* output list size MAXNVX is non-positive. */
154
155 /* 3) The error SPICE(BADDIMENSIONS) is signaled if any element of */
156 /* the grid extents array GRDEXT is non-positive. */
157
158 /* 4) The error SPICE(VERTEXNOTINGRID) is signaled if the ray's */
159 /* vertex is inside, or within a small distance from, the voxel */
160 /* grid. See the description of the parameter GRDTOL. */
161
162 /* 5) The error SPICE(ARRAYTOOSMALL) is signaled if the value of */
163 /* the NVX counter (number of intersected voxels) exceeds the */
164 /* size of the VOXLST input vector. */
165
166 /* $ Files */
167
168 /* None. */
169
170 /* $ Particulars */
171
172 /* This routine supports use of a spatial index for rapid */
173 /* selection of plates that could be hit by a specified ray. */
174
175 /* $ Examples */
176
177 /* See the routine DSKX02 for a usage example. */
178
179 /* $ Restrictions */
180
181 /* None. */
182
183 /* $ Literature_References */
184
185 /* None. */
186
187 /* $ Author_and_Institution */
188
189 /* N.J. Bachman (JPL) */
190 /* J.A. Bytof (JPL) */
191
192 /* $ Version */
193
194 /* - SPICELIB Version 3.1.0, 02-FEB-2016 (NJB) */
195
196 /* Updated to call ZZINGRD rather than INGRD. */
197 /* Minor updates were made to header I/O sections. */
198
199 /* - DSKLIB Version 3.0.0, 11-JUL-2014 (NJB) */
200
201 /* Bug fix: a correction was made to the computation of */
202 /* the vertex offset from the bounding planes of the */
203 /* voxel containing the vertex. */
204
205 /* Minor edits were made to comments. */
206
207 /* Last update was 05-JUN-2014 (NJB) */
208
209 /* Bug fix: the use of the MOD function led to a 1-voxel */
210 /* size error when the input ray's vertex was on the */
211 /* voxel grid boundary. */
212
213 /* An error check for invalid grid dimensions was added. */
214
215 /* Code to prevent arithmetic overflow was added. */
216
217 /* Code was added to prevent the values AX2ERR and AX3ERR from */
218 /* ever becoming negative when the components of the ray's */
219 /* direction vector in the corresponding directions are zero or */
220 /* too small for a voxel step in those directions to occur. */
221
222 /* Renamed the routine's arguments, except for NVX. */
223
224 /* Detailed output descriptions were updated to refer to */
225 /* voxel coordinates rather than IDs. References to sorting */
226 /* were deleted. */
227
228 /* In-line comments now explain the routine's algorithm. */
229 /* Old comments that are no longer applicable were deleted. */
230
231 /* - DSKLIB Version 2.1.0, 26-JUL-2010 (NJB) */
232
233 /* Bug fix: voxel space coordinates of input */
234 /* vertex are now bracketed within the voxel */
235 /* grid. */
236
237 /* This prevents round-off errors from occurring */
238 /* when the vertex is slightly outside the grid, */
239 /* but may not be appropriate for all applications. */
240 /* Therefore it may make sense to make this a */
241 /* private routine. */
242
243 /* - DSKLIB Version 2.0.0, 20-APR-2010 (NJB) */
244
245 /* Removed commented out lines declaring and calling VOX2ID. */
246
247 /* - DSKLIB Version 1.1.0, 08-OCT-2009 (NJB) */
248
249 /* Updated header. */
250
251 /* Bug fix: driving axis for intercept computation is */
252 /* now determined by largest component of ray direction vector. */
253 /* This fix was made long before this header update. */
254
255 /* - DSKLIB Version 1.1.0, 19-OCT-2004 (EDW) */
256
257 /* Added logic to remove duplicate voxel IDs from */
258 /* the return list. Extended programing comments. */
259
260 /* - DSKLIB Version 1.0.1, 26-AUG-2002 (BVS) */
261
262 /* Replaced WRITE with normal error reporting calls. */
263
264 /* - DSKLIB Version 1.0.0, 03-FEB-1999 (JAB) */
265
266 /* -& */
267 /* $ Index_Entries */
268
269 /* list voxels intersected by a ray */
270
271 /* -& */
272
273 /* SPICELIB functions */
274
275
276 /* Local parameters */
277
278
279 /* Local variables */
280
281
282 /* Saved variables */
283
284
285 /* Initial values */
286
287
288 /* Use discovery check-in. */
289
290 if (return_()) {
291 return 0;
292 }
293
294 /* The algorithm below efficiently determines the set of voxels */
295 /* intersected by the input ray. */
296
297 /* This algorithm doesn't compute the intersections of the ray */
298 /* with the boundaries of the voxels, nor does it ever compute */
299 /* the coordinates of any point on the ray. Instead, it keeps */
300 /* track of the voxel boundary planes that the ray passes through. */
301
302 /* The algorithm starts out by determining which voxel contains the */
303 /* ray's vertex. It computes the distances from the vertex to the */
304 /* "next" voxel boundary planes---those that the ray is headed */
305 /* towards. It maintains measurements that enable it to determine */
306 /* which boundary plane is hit next. The voxel on the other side of */
307 /* that intersection point is the next voxel the ray goes through. */
308 /* In the case of ties, any of the candidate "next" voxels may be */
309 /* selected. The "next" voxel is added to the output voxel list, the */
310 /* measurements of relative distances to the next boundaries are */
311 /* updated, and the algorithm continues in this fashion until a */
312 /* voxel outside the grid is detected. */
313
314 /* The relative distance measurements from the ray's vertex to */
315 /* the "next" boundary planes are defined as follows: */
316
317 /* - For the primary ray direction---this is the direction */
318 /* corresponding to the component of largest magnitude of the */
319 /* ray's direction vector---the distance is just the */
320 /* difference of the primary coordinates of the next plane and */
321 /* of the ray's vertex. */
322
323 /* - For each axis orthogonal to the primary one, the algorithm */
324 /* computes the length of the projection onto the primary axis */
325 /* of the portion of the ray extending from the vertex to the */
326 /* "next" voxel boundary plane orthogonal to that non-primary */
327 /* axis. From that projection length the distance from the */
328 /* vertex to the boundary in the primary direction is */
329 /* subtracted. */
330
331 /* For the non-primary axes, these differences are stored in */
332 /* the respective variables */
333
334 /* AX2ERR */
335 /* AX3ERR */
336
337 /* When AX2ERR is negative, the ray will hit the next voxel */
338 /* boundary orthogonal to the "second" axis (having its */
339 /* index stored in the variable IAXIS(2)) before it hits */
340 /* the next boundary orthogonal to the primary axis. The */
341 /* quantity AX3ERR behaves similarly. */
342
343 /* If both AX2ERR and AX3ERR are negative, the more negative */
344 /* value marks the boundary plane that is hit first. */
345
346 /* The axes will be re-labeled using the variable IAXIS. IAXIS(1) */
347 /* will be the index of the primary axis. */
348
349 /* There are a few numeric issues to consider: */
350
351 /* 1) The ratios of the components of the ray's direction vector */
352 /* are computed and stored in the variables S12 and S13. Very */
353 /* small components acting as denominators could cause */
354 /* arithmetic overflow. */
355
356 /* 2) The quantities S12 and S13, while representable as double */
357 /* precision numbers, can be quite large. These quantities */
358 /* may be added repeatedly to the terms AX2ERR and AX3ERR, */
359 /* respectively. These additions could potentially result */
360 /* in arithmetic overflow. */
361
362 /* Both of these problems are addressed by the following observation: */
363
364 /* If a component of the ray direction vector is small enough, and */
365 /* the corresponding component of the ray's vertex is not on a */
366 /* voxel boundary, the ray will exit the grid before reaching a */
367 /* bounding plane orthogonal to that component of the direction */
368 /* vector. */
369
370 /* If the above situation holds, but the ray's vertex is already */
371 /* on a boundary plane orthogonal to the small component, then */
372 /* the ray will exit the grid before hitting a parallel boundary */
373 /* plane. */
374
375 /* So we can safely treat very small direction components as zero. */
376
377
378 /* Check if ray direction vector is a zero vector. */
379
380 if (vzero_(raydir)) {
381 chkin_("XDDA", (ftnlen)4);
382 setmsg_("Ray is the zero vector.", (ftnlen)23);
383 sigerr_("SPICE(ZEROVECTOR)", (ftnlen)17);
384 chkout_("XDDA", (ftnlen)4);
385 return 0;
386 }
387
388 /* Check the voxel grid dimensions. */
389
390 /* Computing MIN */
391 i__1 = min(grdext[0],grdext[1]);
392 if (min(i__1,grdext[2]) < 1) {
393 chkin_("XDDA", (ftnlen)4);
394 setmsg_("Voxel grid dimensions must be strictly positive but are # #"
395 " #.", (ftnlen)62);
396 errint_("#", grdext, (ftnlen)1);
397 errint_("#", &grdext[1], (ftnlen)1);
398 errint_("#", &grdext[2], (ftnlen)1);
399 sigerr_("SPICE(BADDIMENSIONS)", (ftnlen)20);
400 chkout_("XDDA", (ftnlen)4);
401 return 0;
402 }
403
404 /* Make sure the vertex is not too far from the voxel grid. */
405
406 for (i__ = 1; i__ <= 3; ++i__) {
407 if (vertex[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("vertex",
408 i__1, "xdda_", (ftnlen)416)] < grdext[(i__2 = i__ - 1) < 3 &&
409 0 <= i__2 ? i__2 : s_rnge("grdext", i__2, "xdda_", (ftnlen)
410 416)] * -1e-12 || vertex[(i__3 = i__ - 1) < 3 && 0 <= i__3 ?
411 i__3 : s_rnge("vertex", i__3, "xdda_", (ftnlen)416)] > grdext[
412 (i__4 = i__ - 1) < 3 && 0 <= i__4 ? i__4 : s_rnge("grdext",
413 i__4, "xdda_", (ftnlen)416)] * 1.0000000000010001) {
414 chkin_("XDDA", (ftnlen)4);
415 setmsg_("Vertex # # # is outside of voxel grid defined by extent"
416 "s # # #.", (ftnlen)63);
417 errdp_("#", vertex, (ftnlen)1);
418 errdp_("#", &vertex[1], (ftnlen)1);
419 errdp_("#", &vertex[2], (ftnlen)1);
420 errint_("#", grdext, (ftnlen)1);
421 errint_("#", &grdext[1], (ftnlen)1);
422 errint_("#", &grdext[2], (ftnlen)1);
423 sigerr_("SPICE(VERTEXNOTINGRID)", (ftnlen)22);
424 chkout_("XDDA", (ftnlen)4);
425 return 0;
426 }
427 }
428
429 /* The maximum output voxel array size must be positive. */
430
431 if (*maxnvx < 1) {
432 chkin_("XDDA", (ftnlen)4);
433 setmsg_("Maximum voxel list size must be positive but was #.", (
434 ftnlen)51);
435 errint_("#", maxnvx, (ftnlen)1);
436 sigerr_("SPICE(INVALIDSIZE)", (ftnlen)18);
437 chkout_("XDDA", (ftnlen)4);
438 return 0;
439 }
440
441 /* Find the largest component of the direction vector. */
442
443 iaxis[0] = 1;
444 maxcmp = abs(raydir[0]);
445 for (i__ = 2; i__ <= 3; ++i__) {
446 if ((d__1 = raydir[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
447 "raydir", i__1, "xdda_", (ftnlen)459)], abs(d__1)) > maxcmp) {
448 iaxis[0] = i__;
449 maxcmp = (d__1 = raydir[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
450 s_rnge("raydir", i__1, "xdda_", (ftnlen)461)], abs(d__1))
451 ;
452 }
453 }
454
455 /* Set the indices of the orthogonal components of the direction */
456 /* vector. We maintain a right-handed relationship between the axes */
457 /* labeled by IAXIS(1), IAXIS(2), and IAXIS(3): the third axis is */
458 /* the cross product of the first and second. */
459
460 iaxis[1] = next[(i__1 = iaxis[0] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
461 "next", i__1, "xdda_", (ftnlen)472)];
462 iaxis[2] = next[(i__1 = iaxis[1] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
463 "next", i__1, "xdda_", (ftnlen)473)];
464
465 /* Which voxel contains the vertex? Truncate the vertex */
466 /* coordinates to integers. Add 1 to each coord to compensate */
467 /* for 1 based counting. */
468
469 for (i__ = 1; i__ <= 3; ++i__) {
470 intvtx[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("intvtx",
471 i__1, "xdda_", (ftnlen)482)] = (integer) vertex[(i__3 = iaxis[
472 (i__2 = i__ - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("iaxis",
473 i__2, "xdda_", (ftnlen)482)] - 1) < 3 && 0 <= i__3 ? i__3 :
474 s_rnge("vertex", i__3, "xdda_", (ftnlen)482)];
475 icoord[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("icoord",
476 i__1, "xdda_", (ftnlen)484)] = intvtx[(i__2 = i__ - 1) < 3 &&
477 0 <= i__2 ? i__2 : s_rnge("intvtx", i__2, "xdda_", (ftnlen)
478 484)] + 1;
479 icoord[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("icoord",
480 i__1, "xdda_", (ftnlen)486)] = brckti_(&icoord[(i__2 = i__ -
481 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("icoord", i__2, "xdda_", (
482 ftnlen)486)], &c__1, &grdext[(i__4 = iaxis[(i__3 = i__ - 1) <
483 3 && 0 <= i__3 ? i__3 : s_rnge("iaxis", i__3, "xdda_", (
484 ftnlen)486)] - 1) < 3 && 0 <= i__4 ? i__4 : s_rnge("grdext",
485 i__4, "xdda_", (ftnlen)486)]);
486 voxlst[iaxis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("iaxis"
487 , i__1, "xdda_", (ftnlen)488)] - 1] = icoord[(i__2 = i__ - 1)
488 < 3 && 0 <= i__2 ? i__2 : s_rnge("icoord", i__2, "xdda_", (
489 ftnlen)488)];
490 }
491
492 /* Initialize the counter for number of voxels the ray intercepts. */
493 /* The bracketing done above ensures that the coordinates ICOORD of */
494 /* the voxel considered to contain ray's vertex (there is a choice */
495 /* to be made if the vertex lies on a voxel boundary) are within the */
496 /* grid. */
497 *nvx = 1;
498
499 /* Calculate the relative location of vertex within the voxel. The */
500 /* coordinates of a voxel's corners are integer values with each */
501 /* voxel side length 1 (in voxel coords). */
502
503 /* The variable VTXOFF usually has components equal to the */
504 /* fractional parts of the corresponding components of VERTEX( */
505 /* IAXIS(I) ), but the components of VTXOFF may be as large as 1 and */
506 /* are never less than 0. */
507
508 for (i__ = 1; i__ <= 3; ++i__) {
509 d__1 = vertex[(i__3 = iaxis[(i__2 = i__ - 1) < 3 && 0 <= i__2 ? i__2 :
510 s_rnge("iaxis", i__2, "xdda_", (ftnlen)513)] - 1) < 3 && 0 <=
511 i__3 ? i__3 : s_rnge("vertex", i__3, "xdda_", (ftnlen)513)]
512 - (icoord[(i__4 = i__ - 1) < 3 && 0 <= i__4 ? i__4 : s_rnge(
513 "icoord", i__4, "xdda_", (ftnlen)513)] - 1);
514 vtxoff[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("vtxoff",
515 i__1, "xdda_", (ftnlen)513)] = brcktd_(&d__1, &c_b67, &c_b68);
516 }
517
518 /* Compute the lower limit on the magnitudes of RAYDIR( IAXIS(2) ) */
519 /* and of RAYDIR( IAXIS(3) ) for which we'll treat those components */
520 /* of the direction vector as non-zero. */
521
522 limit = 1e-20 / grdext[(i__1 = iaxis[0] - 1) < 3 && 0 <= i__1 ? i__1 :
523 s_rnge("grdext", i__1, "xdda_", (ftnlen)523)] * (d__1 = raydir[(
524 i__2 = iaxis[0] - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("raydir",
525 i__2, "xdda_", (ftnlen)523)], abs(d__1));
526
527 /* If the magnitude of RAYDIR( IAXIS(J) ), J = 2 or 3, is below */
528 /* LIMIT, then the ray can pass through the entire grid in the */
529 /* IAXIS(1) direction without its IAXIS(J) component changing by */
530 /* more than EPS. We'll treat this case as though the IAXIS(J) */
531 /* component of the ray were 0. */
532
533
534 /* Determine the error term initial values and increments. */
535
536
537 ax2err = dpmax_();
538 ax3err = ax2err;
539 s12 = 0.;
540 s13 = 0.;
541
542 /* Compute the initial relative distance measurement AX2ERR */
543 /* for the non-primary axis IAXIS(2). */
544
545 if ((d__1 = raydir[(i__1 = iaxis[1] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
546 "raydir", i__1, "xdda_", (ftnlen)544)], abs(d__1)) > limit) {
547
548 /* For any line segment along the ray, S12 is the ratio of the */
549 /* magnitudes of the projections of the segment in the primary */
550 /* and the IAXIS(2) directions. */
551
552 s12 = (d__1 = raydir[(i__1 = iaxis[0] - 1) < 3 && 0 <= i__1 ? i__1 :
553 s_rnge("raydir", i__1, "xdda_", (ftnlen)550)] / raydir[(i__2 =
554 iaxis[1] - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("raydir",
555 i__2, "xdda_", (ftnlen)550)], abs(d__1));
556 if (raydir[(i__1 = iaxis[0] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
557 "raydir", i__1, "xdda_", (ftnlen)553)] > 0.) {
558
559 /* The primary component of the ray's direction is positive. */
560 /* The distance to the next boundary plane in the primary */
561 /* direction is */
562
563 /* 1.D0 - VTXOFF( IAXIS(1) ) */
564
565 if (raydir[(i__1 = iaxis[1] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
566 "raydir", i__1, "xdda_", (ftnlen)561)] > 0.) {
567
568 /* The IAXIS(2) component of the ray's direction is */
569 /* positive. The distance to the next boundary plane for */
570 /* the that axis is */
571
572 /* 1.D0 - VTXOFF(2) */
573
574 /* The corresponding change along the primary axis is */
575
576 /* S12 * ( 1.D0 - VTXOFF(2) ) */
577
578 /* The "error" term for IAXIS(2) is this value minus the */
579 /* distance from the vertex to the next boundary in the */
580 /* primary direction. */
581
582 ax2err = s12 * (1. - vtxoff[1]) + vtxoff[0] - 1.;
583 } else {
584
585 /* The IAXIS(2) component of the ray's direction is */
586 /* negative. The distance to the next boundary plane for */
587 /* the that axis is */
588
589 /* VTXOFF(2) */
590
591 /* The corresponding change along the primary axis is */
592
593 /* S12 * VTXOFF(2) */
594
595 /* The "error" term for IAXIS(2) is this value minus the */
596 /* distance from the vertex to the next boundary in the */
597 /* primary direction. */
598
599 ax2err = s12 * vtxoff[1] + vtxoff[0] - 1.;
600 }
601 } else {
602
603 /* The primary component of the ray's direction is negative. */
604 /* The distance to the next boundary plane in the primary */
605 /* direction is */
606
607 /* VTXOFF( IAXIS(1) ) */
608 if (raydir[(i__1 = iaxis[1] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
609 "raydir", i__1, "xdda_", (ftnlen)607)] > 0.) {
610
611 /* The IAXIS(2) component of the ray's direction is */
612 /* positive. The distance to the next boundary plane for */
613 /* the that axis is */
614
615 /* 1.D0 - VTXOFF(2) */
616
617 /* The corresponding change along the primary axis is */
618
619 /* S12 * ( 1.D0 - VTXOFF(2) ) */
620
621 /* The "error" term for IAXIS(2) is this value minus the */
622 /* distance from the vertex to the next boundary in the */
623 /* primary direction. */
624 ax2err = s12 * (1. - vtxoff[1]) - vtxoff[0];
625 } else {
626
627 /* The IAXIS(2) component of the ray's direction is */
628 /* negative. The distance to the next boundary plane for */
629 /* the that axis is */
630
631 /* VTXOFF(2) */
632
633 /* The corresponding change along the primary axis is */
634
635 /* S12 * VTXOFF(2) */
636
637 /* The "error" term for IAXIS(2) is this value minus the */
638 /* distance from the vertex to the next boundary in the */
639 /* primary direction. */
640 ax2err = s12 * vtxoff[1] - vtxoff[0];
641 }
642 }
643 }
644
645 /* Computations of AX3ERR are analogous to those of AX2ERR. */
646 /* See the comments above. */
647
648 if ((d__1 = raydir[(i__1 = iaxis[2] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
649 "raydir", i__1, "xdda_", (ftnlen)654)], abs(d__1)) > limit) {
650
651 /* For any line segment along the ray, S13 is the ratio of the */
652 /* magnitudes of the projections of the segment in the primary */
653 /* and the IAXIS(3) directions. */
654
655 s13 = (d__1 = raydir[(i__1 = iaxis[0] - 1) < 3 && 0 <= i__1 ? i__1 :
656 s_rnge("raydir", i__1, "xdda_", (ftnlen)660)] / raydir[(i__2 =
657 iaxis[2] - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge("raydir",
658 i__2, "xdda_", (ftnlen)660)], abs(d__1));
659 if (raydir[(i__1 = iaxis[0] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
660 "raydir", i__1, "xdda_", (ftnlen)662)] > 0.) {
661 if (raydir[(i__1 = iaxis[2] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
662 "raydir", i__1, "xdda_", (ftnlen)664)] > 0.) {
663 ax3err = s13 * (1. - vtxoff[2]) + vtxoff[0] - 1.;
664 } else {
665 ax3err = s13 * vtxoff[2] + vtxoff[0] - 1.;
666 }
667 } else {
668 if (raydir[(i__1 = iaxis[2] - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
669 "raydir", i__1, "xdda_", (ftnlen)672)] > 0.) {
670 ax3err = s13 * (1. - vtxoff[2]) - vtxoff[0];
671 } else {
672 ax3err = s13 * vtxoff[2] - vtxoff[0];
673 }
674 }
675 }
676
677 /* The "steps" set below are the amounts by which any voxel */
678 /* coordinate changes when the "next" voxel is identified. Only one */
679 /* coordinate changes at a time. The magnitude of each coordinate */
680 /* step is always an integer. The signs of the steps are those of */
681 /* the corresponding components of the ray's direction vector. */
682
683 /* We treat direction components smaller than LIMIT as though */
684 /* they were zero. Note that the IAXIS(1) component of the */
685 /* ray will always have magnitude greater than LIMIT. */
686
687 for (i__ = 1; i__ <= 3; ++i__) {
688 if (raydir[(i__2 = iaxis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 :
689 s_rnge("iaxis", i__1, "xdda_", (ftnlen)695)] - 1) < 3 && 0 <=
690 i__2 ? i__2 : s_rnge("raydir", i__2, "xdda_", (ftnlen)695)] >
691 limit) {
692
693 /* Positive component direction, positive step. */
694
695 step[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("step",
696 i__1, "xdda_", (ftnlen)699)] = 1;
697 } else if (raydir[(i__2 = iaxis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ?
698 i__1 : s_rnge("iaxis", i__1, "xdda_", (ftnlen)701)] - 1) < 3
699 && 0 <= i__2 ? i__2 : s_rnge("raydir", i__2, "xdda_", (ftnlen)
700 701)] < -limit) {
701
702 /* Negative component direction, negative step. */
703
704 step[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("step",
705 i__1, "xdda_", (ftnlen)705)] = -1;
706 } else {
707
708 /* No component in this direction, no step. */
709
710 step[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge("step",
711 i__1, "xdda_", (ftnlen)711)] = 0;
712 }
713 }
714
715 /* Follow the ray until it exits the voxel grid. */
716
717 while(zzingrd_(grdext, &voxlst[*nvx * 3 - 3])) {
718 if (ax2err < 0. || ax3err < 0.) {
719
720 /* Ray has crossed over into the next voxel in IAXIS(2) or */
721 /* IAXIS(3) */
722
723 if (ax2err < ax3err) {
724
725 /* The boundary plane orthogonal to axis IAXIS(2) was hit. */
726
727 icoord[1] += step[1];
728 ax2err += s12;
729 ++(*nvx);
730 } else {
731
732 /* The boundary plane orthogonal to axis IAXIS(3) was hit. */
733
734 icoord[2] += step[2];
735 ax3err += s13;
736 ++(*nvx);
737 }
738 } else {
739
740 /* No change in IAXIS(2) or IAXIS(3), step in IAXIS(1). */
741
742 icoord[0] += step[0];
743 ++(*nvx);
744 if (step[1] != 0) {
745 ax2err += -1.;
746 }
747 if (step[2] != 0) {
748 ax3err += -1.;
749 }
750 }
751
752 /* Check we have room in VOXLST. */
753
754 if (*nvx > *maxnvx) {
755 chkin_("XDDA", (ftnlen)4);
756 setmsg_("Index larger than array. Index = #1. Array size = #2.", (
757 ftnlen)53);
758 errint_("#1", nvx, (ftnlen)2);
759 errint_("#2", maxnvx, (ftnlen)2);
760 sigerr_("SPICE(ARRAYTOOSMALL)", (ftnlen)20);
761 chkout_("XDDA", (ftnlen)4);
762 return 0;
763 }
764
765 /* Pack the voxel indices into VOXLST using */
766 /* the values calculated in this loop pass. */
767
768 for (i__ = 1; i__ <= 3; ++i__) {
769 voxlst[iaxis[(i__1 = i__ - 1) < 3 && 0 <= i__1 ? i__1 : s_rnge(
770 "iaxis", i__1, "xdda_", (ftnlen)784)] + *nvx * 3 - 4] =
771 icoord[(i__2 = i__ - 1) < 3 && 0 <= i__2 ? i__2 : s_rnge(
772 "icoord", i__2, "xdda_", (ftnlen)784)];
773 }
774 }
775
776 /* Subtract one off the voxel count since the final voxel */
777 /* exists outside the grid. */
778
779 --(*nvx);
780 return 0;
781 } /* xdda_ */
782
783