BRL-CAD
tor.c
Go to the documentation of this file.
1 /* T O R . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1985-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @addtogroup primitives */
21 /** @{ */
22 /** @file primitives/tor/tor.c
23  *
24  * Intersect a ray with a Torus
25  *
26  * Algorithm:
27  *
28  * Given V, H, A, and B, there is a set of points on this torus
29  *
30  * { (x, y, z) | (x, y, z) is on torus defined by V, H, A, B }
31  *
32  * Through a series of Transformations, this set will be transformed
33  * into a set of points on a unit torus (R1==1) centered at the origin
34  * which lies on the X-Y plane (i.e., H is on the Z axis).
35  *
36  * { (x', y', z') | (x', y', z') is on unit torus at origin }
37  *
38  * The transformation from X to X' is accomplished by:
39  *
40  * X' = S(R(X - V))
41  *
42  * where R(X) = (A/(|A|))
43  * (B/(|B|)) . X
44  * (H/(|H|))
45  *
46  * and S(X) = (1/|A| 0 0)
47  * (0 1/|A| 0) . X
48  * (0 0 1/|A|)
49  * where |A| = R1
50  *
51  * To find the intersection of a line with the torus, consider the
52  * parametric line L:
53  *
54  * L : { P(n) | P + t(n) . D }
55  *
56  * Call W the actual point of intersection between L and the torus.
57  * Let W' be the point of intersection between L' and the unit torus.
58  *
59  * L' : { P'(n) | P' + t(n) . D' }
60  *
61  * W = invR(invS(W')) + V
62  *
63  * Where W' = k D' + P'.
64  *
65  * Given a line and a ratio, alpha, finds the equation of the unit
66  * torus in terms of the variable 't'.
67  *
68  * The equation for the torus is:
69  *
70  * [ X**2 + Y**2 + Z**2 + (1 - alpha**2) ]**2 - 4*(X**2 + Y**2) = 0
71  *
72  * First, find X, Y, and Z in terms of 't' for this line, then
73  * substitute them into the equation above.
74  *
75  * Wx = Dx*t + Px
76  *
77  * Wx**2 = Dx**2 * t**2 + 2 * Dx * Px + Px**2
78  *
79  * The real roots of the equation in 't' are the intersect points
80  * along the parametric line.
81  *
82  * NORMALS. Given the point W on the torus, what is the vector normal
83  * to the tangent plane at that point?
84  *
85  * Map W onto the unit torus, i.e.: W' = S(R(W - V)). In this case,
86  * we find W' by solving the parametric line given k.
87  *
88  * The gradient of the torus at W' is in fact the normal vector.
89  *
90  * Given that the equation for the unit torus is:
91  *
92  * [ X**2 + Y**2 + Z**2 + (1 - alpha**2) ]**2 - 4*(X**2 + Y**2) = 0
93  *
94  * let w = X**2 + Y**2 + Z**2 + (1 - alpha**2), then the equation becomes:
95  *
96  * w**2 - 4*(X**2 + Y**2) = 0
97  *
98  * For f(x, y, z) = 0, the gradient of f() is (df/dx, df/dy, df/dz).
99  *
100  * df/dx = 2 * w * 2 * x - 8 * x = (4 * w - 8) * x
101  * df/dy = 2 * w * 2 * y - 8 * y = (4 * w - 8) * y
102  * df/dz = 2 * w * 2 * z = 4 * w * z
103  *
104  * Note that the normal vector produced above will not have unit
105  * length. Also, to make this useful for the original torus, it will
106  * have to be rotated back to the orientation of the original torus.
107  *
108  */
109 /** @} */
110 
111 #include "common.h"
112 
113 #include <stddef.h>
114 #include <string.h>
115 #include <math.h>
116 #include "bio.h"
117 
118 #include "bu/cv.h"
119 #include "vmath.h"
120 #include "db.h"
121 #include "nmg.h"
122 #include "rtgeom.h"
123 #include "raytrace.h"
124 
125 #include "../../librt_private.h"
126 
127 
128 /*
129  * The TORUS has the following input fields:
130  * V V from origin to center
131  * H Radius Vector, Normal to plane of torus. |H| = R2
132  * A, B perpendicular, to CENTER of torus. |A|==|B|==R1
133  * F5, F6 perpendicular, for inner edge (unused)
134  * F7, F8 perpendicular, for outer edge (unused)
135  *
136  */
137 
138 const struct bu_structparse rt_tor_parse[] = {
139  {"%f", 3, "V", bu_offsetofarray(struct rt_tor_internal, v, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
140  {"%f", 3, "H", bu_offsetofarray(struct rt_tor_internal, h, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
141  {"%f", 1, "r_a", bu_offsetof(struct rt_tor_internal, r_a), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
142  {"%f", 1, "r_h", bu_offsetof(struct rt_tor_internal, r_h), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL},
143  {{'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL}
144 };
145 
146 
147 struct tor_specific {
148  fastf_t tor_alpha; /* 0 < (R2/R1) <= 1 */
149  fastf_t tor_r1; /* for inverse scaling of k values. */
150  fastf_t tor_r2; /* for curvature */
151  vect_t tor_V; /* Vector to center of torus */
152  vect_t tor_N; /* unit normal to plane of torus */
153  mat_t tor_SoR; /* Scale(Rot(vect)) */
154  mat_t tor_invR; /* invRot(vect') */
155 };
156 
157 /**
158  * Compute the bounding RPP for a circular torus.
159  */
160 int
161 rt_tor_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
162  vect_t P, w1; /* for RPP calculation */
163  fastf_t f;
164  struct rt_tor_internal *tip = (struct rt_tor_internal *)ip->idb_ptr;
165  RT_TOR_CK_MAGIC(tip);
166 
167  /* Compute the bounding RPP planes for a circular torus.
168  *
169  * Given a circular torus with vertex V, vector N, and radii r1
170  * and r2. A bounding plane with direction vector P will touch
171  * the surface of the torus at the points:
172  *
173  * V +/- [r2 + r1 * |N x P|] P
174  */
175  /* X */
176  VSET(P, 1.0, 0, 0); /* bounding plane normal */
177  VCROSS(w1, tip->h, P); /* for sin(angle N P) */
178  f = tip->r_h + tip->r_a * MAGNITUDE(w1);
179  VSCALE(w1, P, f);
180  f = fabs(w1[X]);
181  (*min)[X] = tip->v[X] - f;
182  (*max)[X] = tip->v[X] + f;
183 
184  /* Y */
185  VSET(P, 0, 1.0, 0); /* bounding plane normal */
186  VCROSS(w1, tip->h, P); /* for sin(angle N P) */
187  f = tip->r_h + tip->r_a * MAGNITUDE(w1);
188  VSCALE(w1, P, f);
189  f = fabs(w1[Y]);
190  (*min)[Y] = tip->v[Y] - f;
191  (*max)[Y] = tip->v[Y] + f;
192 
193  /* Z */
194  VSET(P, 0, 0, 1.0); /* bounding plane normal */
195  VCROSS(w1, tip->h, P); /* for sin(angle N P) */
196  f = tip->r_h + tip->r_a * MAGNITUDE(w1);
197  VSCALE(w1, P, f);
198  f = fabs(w1[Z]);
199  (*min)[Z] = tip->v[Z] - f;
200  (*max)[Z] = tip->v[Z] + f;
201 
202  return 0;
203 }
204 
205 
206 /**
207  * Given a pointer to a GED database record, and a transformation
208  * matrix, determine if this is a valid torus, and if so, precompute
209  * various terms of the formula.
210  *
211  * Returns -
212  * 0 TOR is OK
213  * !0 Error in description
214  *
215  * Implicit return -
216  * A struct tor_specific is created, and its address is stored in
217  * stp->st_specific for use by rt_tor_shot().
218  */
219 int
220 rt_tor_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
221 {
222  register struct tor_specific *tor;
223  struct rt_tor_internal *tip;
224 
225  mat_t R;
226  fastf_t f;
227 
228  tip = (struct rt_tor_internal *)ip->idb_ptr;
229  RT_TOR_CK_MAGIC(tip);
230 
231  /* Validate that |A| == |B| (for now) */
232  if (!NEAR_EQUAL(tip->r_a, tip->r_b, 0.001)) {
233  bu_log("tor(%s): (|A|=%f) != (|B|=%f) \n",
234  stp->st_name, tip->r_a, tip->r_b);
235  return 1; /* BAD */
236  }
237 
238  /* Validate that A.B == 0, B.H == 0, A.H == 0 */
239  f = VDOT(tip->a, tip->b)/(tip->r_a*tip->r_b);
240 
241  if (!NEAR_ZERO(f, rtip->rti_tol.dist)) {
242  bu_log("tor(%s): A not perpendicular to B, f=%f\n",
243  stp->st_name, f);
244  return 1; /* BAD */
245  }
246  f = VDOT(tip->b, tip->h)/(tip->r_b);
247  if (!NEAR_ZERO(f, rtip->rti_tol.dist)) {
248  bu_log("tor(%s): B not perpendicular to H, f=%f\n",
249  stp->st_name, f);
250  return 1; /* BAD */
251  }
252  f = VDOT(tip->a, tip->h)/(tip->r_a);
253  if (!NEAR_ZERO(f, rtip->rti_tol.dist)) {
254  bu_log("tor(%s): A not perpendicular to H, f=%f\n",
255  stp->st_name, f);
256  return 1; /* BAD */
257  }
258 
259  /* Validate that 0 < r2 <= r1 for alpha computation */
260  if (0.0 >= tip->r_h || tip->r_h > tip->r_a) {
261  bu_log("r1 = %f, r2 = %f\n", tip->r_a, tip->r_h);
262  bu_log("tor(%s): 0 < r2 <= r1 is not true\n", stp->st_name);
263  return 1; /* BAD */
264  }
265 
266  /* Solid is OK, compute constant terms now */
267  BU_GET(tor, struct tor_specific);
268  stp->st_specific = (void *)tor;
269 
270  tor->tor_r1 = tip->r_a;
271  tor->tor_r2 = tip->r_h;
272 
273  VMOVE(tor->tor_V, tip->v);
274  tor->tor_alpha = tip->r_h/tor->tor_r1;
275 
276  /* Compute R and invR matrices */
277  VMOVE(tor->tor_N, tip->h);
278 
279  MAT_IDN(R);
280  VSCALE(&R[0], tip->a, 1.0/tip->r_a);
281  VSCALE(&R[4], tip->b, 1.0/tip->r_b);
282  VMOVE(&R[8], tor->tor_N);
283  bn_mat_inv(tor->tor_invR, R);
284 
285  /* Compute SoR. Here, S = I / r1 */
286  MAT_COPY(tor->tor_SoR, R);
287  tor->tor_SoR[15] *= tor->tor_r1;
288 
289  VMOVE(stp->st_center, tor->tor_V);
290  stp->st_aradius = stp->st_bradius = tor->tor_r1 + tip->r_h;
291 
292  if (stp->st_meth->ft_bbox(ip, &(stp->st_min), &(stp->st_max), &(rtip->rti_tol))) {
293  return 1;
294  }
295 
296  return 0; /* OK */
297 }
298 
299 
300 void
301 rt_tor_print(register const struct soltab *stp)
302 {
303  register const struct tor_specific *tor =
304  (struct tor_specific *)stp->st_specific;
305 
306  bu_log("r2/r1 (alpha) = %f\n", tor->tor_alpha);
307  bu_log("r1 = %f\n", tor->tor_r1);
308  bu_log("r2 = %f\n", tor->tor_r2);
309  VPRINT("V", tor->tor_V);
310  VPRINT("N", tor->tor_N);
311  bn_mat_print("S o R", tor->tor_SoR);
312  bn_mat_print("invR", tor->tor_invR);
313 }
314 
315 
316 /**
317  * Intersect a ray with an torus, where all constant terms have been
318  * precomputed by rt_tor_prep(). If an intersection occurs, one or
319  * two struct seg(s) will be acquired and filled in.
320  *
321  * NOTE: All lines in this function are represented parametrically by
322  * a point, P(x0, y0, z0) and a direction normal, D = ax + by + cz.
323  * Any point on a line can be expressed by one variable 't', where
324  *
325  * X = a*t + x0, e.g., X = Dx*t + Px
326  * Y = b*t + y0,
327  * Z = c*t + z0.
328  *
329  * First, convert the line to the coordinate system of a "standard"
330  * torus. This is a torus which lies in the X-Y plane, circles the
331  * origin, and whose primary radius is one. The secondary radius is
332  * alpha = (R2/R1) of the original torus where (0 < alpha <= 1).
333  *
334  * Then find the equation of that line and the standard torus, which
335  * turns out to be a quartic equation in 't'. Solve the equation
336  * using a general polynomial root finder. Use those values of 't' to
337  * compute the points of intersection in the original coordinate
338  * system.
339  *
340  * Returns -
341  * 0 MISS
342  * >0 HIT
343  */
344 int
345 rt_tor_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
346 {
347  register struct tor_specific *tor =
348  (struct tor_specific *)stp->st_specific;
349  register struct seg *segp;
350  vect_t dprime; /* D' */
351  vect_t pprime; /* P' */
352  vect_t work; /* temporary vector */
353  bn_poly_t C; /* The final equation */
354  bn_complex_t val[4]; /* The complex roots */
355  double k[4]; /* The real roots */
356  register int i;
357  int j;
358  bn_poly_t A, Asqr;
359  bn_poly_t X2_Y2; /* X**2 + Y**2 */
360  vect_t cor_pprime; /* new ray origin */
361  fastf_t cor_proj;
362 
363  /* Convert vector into the space of the unit torus */
364  MAT4X3VEC(dprime, tor->tor_SoR, rp->r_dir);
365  VUNITIZE(dprime);
366 
367  VSUB2(work, rp->r_pt, tor->tor_V);
368  MAT4X3VEC(pprime, tor->tor_SoR, work);
369 
370  /* normalize distance from torus. substitute corrected pprime
371  * which contains a translation along ray direction to closest
372  * approach to vertex of torus. Translating ray origin along
373  * direction of ray to closest pt. to origin of solid's coordinate
374  * system, new ray origin is 'cor_pprime'.
375  */
376  cor_proj = VDOT(pprime, dprime);
377  VSCALE(cor_pprime, dprime, cor_proj);
378  VSUB2(cor_pprime, pprime, cor_pprime);
379 
380  /* Given a line and a ratio, alpha, finds the equation of the unit
381  * torus in terms of the variable 't'.
382  *
383  * The equation for the torus is:
384  *
385  * [ X**2 + Y**2 + Z**2 + (1 - alpha**2) ]**2 - 4*(X**2 + Y**2) = 0
386  *
387  * First, find X, Y, and Z in terms of 't' for this line, then
388  * substitute them into the equation above.
389  *
390  * Wx = Dx*t + Px
391  *
392  * Wx**2 = Dx**2 * t**2 + 2 * Dx * Px + Px**2
393  * [0] [1] [2] dgr=2
394  */
395  X2_Y2.dgr = 2;
396  X2_Y2.cf[0] = dprime[X] * dprime[X] + dprime[Y] * dprime[Y];
397  X2_Y2.cf[1] = 2.0 * (dprime[X] * cor_pprime[X] +
398  dprime[Y] * cor_pprime[Y]);
399  X2_Y2.cf[2] = cor_pprime[X] * cor_pprime[X] +
400  cor_pprime[Y] * cor_pprime[Y];
401 
402  /* A = X2_Y2 + Z2 */
403  A.dgr = 2;
404  A.cf[0] = X2_Y2.cf[0] + dprime[Z] * dprime[Z];
405  A.cf[1] = X2_Y2.cf[1] + 2.0 * dprime[Z] * cor_pprime[Z];
406  A.cf[2] = X2_Y2.cf[2] + cor_pprime[Z] * cor_pprime[Z] +
407  1.0 - tor->tor_alpha * tor->tor_alpha;
408 
409  /* Inline expansion of (void) bn_poly_mul(&Asqr, &A, &A) */
410  /* Both polys have degree two */
411  Asqr.dgr = 4;
412  Asqr.cf[0] = A.cf[0] * A.cf[0];
413  Asqr.cf[1] = A.cf[0] * A.cf[1] + A.cf[1] * A.cf[0];
414  Asqr.cf[2] = A.cf[0] * A.cf[2] + A.cf[1] * A.cf[1] + A.cf[2] * A.cf[0];
415  Asqr.cf[3] = A.cf[1] * A.cf[2] + A.cf[2] * A.cf[1];
416  Asqr.cf[4] = A.cf[2] * A.cf[2];
417 
418  /* Inline expansion of bn_poly_scale(&X2_Y2, 4.0) and
419  * bn_poly_sub(&C, &Asqr, &X2_Y2).
420  */
421  C.dgr = 4;
422  C.cf[0] = Asqr.cf[0];
423  C.cf[1] = Asqr.cf[1];
424  C.cf[2] = Asqr.cf[2] - X2_Y2.cf[0] * 4.0;
425  C.cf[3] = Asqr.cf[3] - X2_Y2.cf[1] * 4.0;
426  C.cf[4] = Asqr.cf[4] - X2_Y2.cf[2] * 4.0;
427 
428  /* It is known that the equation is 4th order. Therefore, if the
429  * root finder returns other than 4 roots, error.
430  */
431  if ((i = rt_poly_roots(&C, val, stp->st_dp->d_namep)) != 4) {
432  if (i > 0) {
433  bu_log("tor: rt_poly_roots() 4!=%d\n", i);
434  bn_pr_roots(stp->st_name, val, i);
435  } else if (i < 0) {
436  static int reported=0;
437  bu_log("The root solver failed to converge on a solution for %s\n", stp->st_dp->d_namep);
438  if (!reported) {
439  VPRINT("while shooting from:\t", rp->r_pt);
440  VPRINT("while shooting at:\t", rp->r_dir);
441  bu_log("Additional torus convergence failure details will be suppressed.\n");
442  reported=1;
443  }
444  }
445  return 0; /* MISS */
446  }
447 
448  /* Only real roots indicate an intersection in real space.
449  *
450  * Look at each root returned; if the imaginary part is zero or
451  * sufficiently close, then use the real part as one value of 't'
452  * for the intersections
453  */
454  for (j=0, i=0; j < 4; j++) {
455  if (NEAR_ZERO(val[j].im, ap->a_rt_i->rti_tol.dist))
456  k[i++] = val[j].re;
457  }
458 
459  /* reverse above translation by adding distance to all 'k' values.
460  */
461  for (j = 0; j < i; ++j)
462  k[j] -= cor_proj;
463 
464  /* Here, 'i' is number of points found */
465  switch (i) {
466  case 0:
467  return 0; /* No hit */
468 
469  default:
470  bu_log("rt_tor_shot: reduced 4 to %d roots\n", i);
471  bn_pr_roots(stp->st_name, val, 4);
472  return 0; /* No hit */
473 
474  case 2:
475  {
476  /* Sort most distant to least distant. */
477  fastf_t u;
478  if ((u=k[0]) < k[1]) {
479  /* bubble larger towards [0] */
480  k[0] = k[1];
481  k[1] = u;
482  }
483  }
484  break;
485  case 4:
486  {
487  register short n;
488  register short lim;
489 
490  /* Inline rt_pt_sort(). Sorts k[] into descending order. */
491  for (lim = i-1; lim > 0; lim--) {
492  for (n = 0; n < lim; n++) {
493  fastf_t u;
494  if ((u=k[n]) < k[n+1]) {
495  /* bubble larger towards [0] */
496  k[n] = k[n+1];
497  k[n+1] = u;
498  }
499  }
500  }
501  }
502  break;
503  }
504 
505  /* Now, t[0] > t[npts-1] */
506  /* k[1] is entry point, and k[0] is farthest exit point */
507  RT_GET_SEG(segp, ap->a_resource);
508  segp->seg_stp = stp;
509  segp->seg_in.hit_dist = k[1]*tor->tor_r1;
510  segp->seg_out.hit_dist = k[0]*tor->tor_r1;
511  segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 0;
512  /* Set aside vector for rt_tor_norm() later */
513  VJOIN1(segp->seg_in.hit_vpriv, pprime, k[1], dprime);
514  VJOIN1(segp->seg_out.hit_vpriv, pprime, k[0], dprime);
515  BU_LIST_INSERT(&(seghead->l), &(segp->l));
516 
517  if (i == 2)
518  return 2; /* HIT */
519 
520  /* 4 points */
521  /* k[3] is entry point, and k[2] is exit point */
522  RT_GET_SEG(segp, ap->a_resource);
523  segp->seg_stp = stp;
524  segp->seg_in.hit_dist = k[3]*tor->tor_r1;
525  segp->seg_out.hit_dist = k[2]*tor->tor_r1;
526  segp->seg_in.hit_surfno = segp->seg_out.hit_surfno = 1;
527  VJOIN1(segp->seg_in.hit_vpriv, pprime, k[3], dprime);
528  VJOIN1(segp->seg_out.hit_vpriv, pprime, k[2], dprime);
529  BU_LIST_INSERT(&(seghead->l), &(segp->l));
530  return 4; /* HIT */
531 }
532 
533 
534 #define RT_TOR_SEG_MISS(SEG) (SEG).seg_stp=(struct soltab *) 0;
535 /**
536  * This is the Becker vector version
537  */
538 void
539 rt_tor_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
540  /* An array of solid pointers */
541  /* An array of ray pointers */
542  /* array of segs (results returned) */
543  /* Number of ray/object pairs */
544 
545 {
546  register int i;
547  register struct tor_specific *tor;
548  vect_t dprime; /* D' */
549  vect_t pprime; /* P' */
550  vect_t work; /* temporary vector */
551  bn_poly_t *C; /* The final equation */
552  bn_complex_t (*val)[4]; /* The complex roots */
553  int num_roots;
554  int num_zero;
555  bn_poly_t A, Asqr;
556  bn_poly_t X2_Y2; /* X**2 + Y**2 */
557  vect_t cor_pprime; /* new ray origin */
558  fastf_t *cor_proj;
559 
560  /* Allocate space for polys and roots */
561  C = (bn_poly_t *)bu_malloc(n * sizeof(bn_poly_t), "tor bn_poly_t");
562  val = (bn_complex_t (*)[4])bu_malloc(n * sizeof(bn_complex_t) * 4,
563  "tor bn_complex_t");
564  cor_proj = (fastf_t *)bu_malloc(n * sizeof(fastf_t), "tor proj");
565 
566  /* Initialize seg_stp to assume hit (zero will then flag miss) */
567  for (i = 0; i < n; i++) segp[i].seg_stp = stp[i];
568 
569  /* for each ray/torus pair */
570  for (i = 0; i < n; i++) {
571  if (segp[i].seg_stp == 0) continue; /* Skip */
572  tor = (struct tor_specific *)stp[i]->st_specific;
573 
574  /* Convert vector into the space of the unit torus */
575  MAT4X3VEC(dprime, tor->tor_SoR, rp[i]->r_dir);
576  VUNITIZE(dprime);
577 
578  /* Use segp[i].seg_in.hit_normal as tmp to hold dprime */
579  VMOVE(segp[i].seg_in.hit_normal, dprime);
580 
581  VSUB2(work, rp[i]->r_pt, tor->tor_V);
582  MAT4X3VEC(pprime, tor->tor_SoR, work);
583 
584  /* Use segp[i].seg_out.hit_normal as tmp to hold pprime */
585  VMOVE(segp[i].seg_out.hit_normal, pprime);
586 
587  /* normalize distance from torus. substitute corrected pprime
588  * which contains a translation along ray direction to closest
589  * approach to vertex of torus. Translating ray origin along
590  * direction of ray to closest pt. to origin of solid's
591  * coordinate system, new ray origin is 'cor_pprime'.
592  */
593  cor_proj[i] = VDOT(pprime, dprime);
594  VSCALE(cor_pprime, dprime, cor_proj[i]);
595  VSUB2(cor_pprime, pprime, cor_pprime);
596 
597  /*
598  * Given a line and a ratio, alpha, finds the equation of the
599  * unit torus in terms of the variable 't'.
600  *
601  * The equation for the torus is:
602  *
603  * [X**2 + Y**2 + Z**2 + (1 - alpha**2)]**2 - 4*(X**2 + Y**2) =0
604  *
605  * First, find X, Y, and Z in terms of 't' for this line, then
606  * substitute them into the equation above.
607  *
608  * Wx = Dx*t + Px
609  *
610  * Wx**2 = Dx**2 * t**2 + 2 * Dx * Px + Px**2
611  * [0] [1] [2] dgr=2
612  */
613  X2_Y2.dgr = 2;
614  X2_Y2.cf[0] = dprime[X] * dprime[X] + dprime[Y] * dprime[Y];
615  X2_Y2.cf[1] = 2.0 * (dprime[X] * cor_pprime[X] +
616  dprime[Y] * cor_pprime[Y]);
617  X2_Y2.cf[2] = cor_pprime[X] * cor_pprime[X] +
618  cor_pprime[Y] * cor_pprime[Y];
619 
620  /* A = X2_Y2 + Z2 */
621  A.dgr = 2;
622  A.cf[0] = X2_Y2.cf[0] + dprime[Z] * dprime[Z];
623  A.cf[1] = X2_Y2.cf[1] + 2.0 * dprime[Z] * cor_pprime[Z];
624  A.cf[2] = X2_Y2.cf[2] + cor_pprime[Z] * cor_pprime[Z] +
625  1.0 - tor->tor_alpha * tor->tor_alpha;
626 
627  /* Inline expansion of (void) bn_poly_mul(&Asqr, &A, &A) */
628  /* Both polys have degree two */
629  Asqr.dgr = 4;
630  Asqr.cf[0] = A.cf[0] * A.cf[0];
631  Asqr.cf[1] = A.cf[0] * A.cf[1] +
632  A.cf[1] * A.cf[0];
633  Asqr.cf[2] = A.cf[0] * A.cf[2] +
634  A.cf[1] * A.cf[1] +
635  A.cf[2] * A.cf[0];
636  Asqr.cf[3] = A.cf[1] * A.cf[2] +
637  A.cf[2] * A.cf[1];
638  Asqr.cf[4] = A.cf[2] * A.cf[2];
639 
640  /* Inline expansion of (void) bn_poly_scale(&X2_Y2, 4.0) */
641  X2_Y2.cf[0] *= 4.0;
642  X2_Y2.cf[1] *= 4.0;
643  X2_Y2.cf[2] *= 4.0;
644 
645  /* Inline expansion of (void) bn_poly_sub(&C, &Asqr, &X2_Y2) */
646  /* offset is known to be 2 */
647  C[i].dgr = 4;
648  C[i].cf[0] = Asqr.cf[0];
649  C[i].cf[1] = Asqr.cf[1];
650  C[i].cf[2] = Asqr.cf[2] - X2_Y2.cf[0];
651  C[i].cf[3] = Asqr.cf[3] - X2_Y2.cf[1];
652  C[i].cf[4] = Asqr.cf[4] - X2_Y2.cf[2];
653  }
654 
655  /* Unfortunately finding the 4th order roots are too ugly to
656  * expand the root solving manually.
657  */
658  for (i = 0; i < n; i++) {
659  if (segp[i].seg_stp == 0) continue; /* Skip */
660 
661  /* It is known that the equation is 4th order. Therefore, if
662  * the root finder returns other than 4 roots, error.
663  */
664  if ((num_roots = rt_poly_roots(&(C[i]), &(val[i][0]), (*stp)->st_dp->d_namep)) != 4) {
665  if (num_roots > 0) {
666  bu_log("tor: rt_poly_roots() 4!=%d\n", num_roots);
667  bn_pr_roots("tor", val[i], num_roots);
668  } else if (num_roots < 0) {
669  static int reported=0;
670  bu_log("The root solver failed to converge on a solution for %s\n", stp[i]->st_dp->d_namep);
671  if (!reported) {
672  VPRINT("while shooting from:\t", rp[i]->r_pt);
673  VPRINT("while shooting at:\t", rp[i]->r_dir);
674  bu_log("Additional torus convergence failure details will be suppressed.\n");
675  reported=1;
676  }
677  }
678  RT_TOR_SEG_MISS(segp[i]); /* MISS */
679  }
680  }
681 
682  /* for each ray/torus pair */
683  for (i = 0; i < n; i++) {
684  if (segp[i].seg_stp == 0) continue; /* Skip */
685 
686  /* Only real roots indicate an intersection in real space.
687  *
688  * Look at each root returned; if the imaginary part is zero
689  * or sufficiently close, then use the real part as one value
690  * of 't' for the intersections
691  */
692  /* Also reverse translation by adding distance to all 'k' values. */
693  /* Reuse C to hold k values */
694  num_zero = 0;
695  if (NEAR_ZERO(val[i][0].im, ap->a_rt_i->rti_tol.dist))
696  C[i].cf[num_zero++] = val[i][0].re - cor_proj[i];
697  if (NEAR_ZERO(val[i][1].im, ap->a_rt_i->rti_tol.dist)) {
698  C[i].cf[num_zero++] = val[i][1].re - cor_proj[i];
699  }
700  if (NEAR_ZERO(val[i][2].im, ap->a_rt_i->rti_tol.dist)) {
701  C[i].cf[num_zero++] = val[i][2].re - cor_proj[i];
702  }
703  if (NEAR_ZERO(val[i][3].im, ap->a_rt_i->rti_tol.dist)) {
704  C[i].cf[num_zero++] = val[i][3].re - cor_proj[i];
705  }
706  C[i].dgr = num_zero;
707 
708  /* Here, 'i' is number of points found */
709  if (num_zero == 0) {
710  RT_TOR_SEG_MISS(segp[i]); /* MISS */
711  } else if (num_zero != 2 && num_zero != 4) {
712  RT_TOR_SEG_MISS(segp[i]); /* MISS */
713  }
714  }
715 
716  /* Process each one segment hit */
717  for (i = 0; i < n; i++) {
718  if (segp[i].seg_stp == 0) continue; /* Skip */
719  if (C[i].dgr != 2) continue; /* Not one segment */
720 
721  tor = (struct tor_specific *)stp[i]->st_specific;
722 
723  /* segp[i].seg_in.hit_normal holds dprime */
724  /* segp[i].seg_out.hit_normal holds pprime */
725  if (C[i].cf[1] < C[i].cf[0]) {
726  /* C[i].cf[1] is entry point */
727  segp[i].seg_in.hit_dist = C[i].cf[1]*tor->tor_r1;
728  segp[i].seg_out.hit_dist = C[i].cf[0]*tor->tor_r1;
729  /* Set aside vector for rt_tor_norm() later */
730  VJOIN1(segp[i].seg_in.hit_vpriv,
731  segp[i].seg_out.hit_normal,
732  C[i].cf[1], segp[i].seg_in.hit_normal);
733  VJOIN1(segp[i].seg_out.hit_vpriv,
734  segp[i].seg_out.hit_normal,
735  C[i].cf[0], segp[i].seg_in.hit_normal);
736  } else {
737  /* C[i].cf[0] is entry point */
738  segp[i].seg_in.hit_dist = C[i].cf[0]*tor->tor_r1;
739  segp[i].seg_out.hit_dist = C[i].cf[1]*tor->tor_r1;
740  /* Set aside vector for rt_tor_norm() later */
741  VJOIN1(segp[i].seg_in.hit_vpriv,
742  segp[i].seg_out.hit_normal,
743  C[i].cf[0], segp[i].seg_in.hit_normal);
744  VJOIN1(segp[i].seg_out.hit_vpriv,
745  segp[i].seg_out.hit_normal,
746  C[i].cf[1], segp[i].seg_in.hit_normal);
747  }
748  }
749 
750  /* Process each two segment hit */
751  for (i = 0; i < n; i++) {
752 
753  if (segp[i].seg_stp == 0) continue; /* Skip */
754  if (C[i].dgr != 4) continue; /* Not two segment */
755 
756  tor = (struct tor_specific *)stp[i]->st_specific;
757 
758  /* Sort most distant to least distant. */
759  rt_pt_sort(C[i].cf, 4);
760  /* Now, t[0] > t[npts-1] */
761 
762  /* segp[i].seg_in.hit_normal holds dprime */
763  VMOVE(dprime, segp[i].seg_in.hit_normal);
764  /* segp[i].seg_out.hit_normal holds pprime */
765  VMOVE(pprime, segp[i].seg_out.hit_normal);
766 
767  /* C[i].cf[1] is entry point */
768  segp[i].seg_in.hit_dist = C[i].cf[1]*tor->tor_r1;
769  segp[i].seg_out.hit_dist = C[i].cf[0]*tor->tor_r1;
770  /* Set aside vector for rt_tor_norm() later */
771  VJOIN1(segp[i].seg_in.hit_vpriv, pprime, C[i].cf[1], dprime);
772  VJOIN1(segp[i].seg_out.hit_vpriv, pprime, C[i].cf[0], dprime);
773  }
774 
775  /* Free tmp space used */
776  bu_free((char *)C, "tor C");
777  bu_free((char *)val, "tor val");
778  bu_free((char *)cor_proj, "tor cor_proj");
779 }
780 
781 
782 /**
783  * Compute the normal to the torus, given a point on the UNIT TORUS
784  * centered at the origin on the X-Y plane. The gradient of the torus
785  * at that point is in fact the normal vector, which will have to be
786  * given unit length. To make this useful for the original torus, it
787  * will have to be rotated back to the orientation of the original
788  * torus.
789  *
790  * Given that the equation for the unit torus is:
791  *
792  * [ X**2 + Y**2 + Z**2 + (1 - alpha**2) ]**2 - 4*(X**2 + Y**2) = 0
793  *
794  * let w = X**2 + Y**2 + Z**2 + (1 - alpha**2), then the equation becomes:
795  *
796  * w**2 - 4*(X**2 + Y**2) = 0
797  *
798  * For f(x, y, z) = 0, the gradient of f() is (df/dx, df/dy, df/dz).
799  *
800  * df/dx = 2 * w * 2 * x - 8 * x = (4 * w - 8) * x
801  * df/dy = 2 * w * 2 * y - 8 * y = (4 * w - 8) * y
802  * df/dz = 2 * w * 2 * z = 4 * w * z
803  *
804  * Since we rescale the gradient (normal) to unity, we divide the
805  * above equations by four here.
806  */
807 void
808 rt_tor_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
809 {
810  register struct tor_specific *tor =
811  (struct tor_specific *)stp->st_specific;
812 
813  fastf_t w;
814  vect_t work;
815 
816  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
817  w = hitp->hit_vpriv[X]*hitp->hit_vpriv[X] +
818  hitp->hit_vpriv[Y]*hitp->hit_vpriv[Y] +
819  hitp->hit_vpriv[Z]*hitp->hit_vpriv[Z] +
820  1.0 - tor->tor_alpha*tor->tor_alpha;
821  VSET(work,
822  (w - 2.0) * hitp->hit_vpriv[X],
823  (w - 2.0) * hitp->hit_vpriv[Y],
824  w * hitp->hit_vpriv[Z]);
825  VUNITIZE(work);
826 
827  MAT3X3VEC(hitp->hit_normal, tor->tor_invR, work);
828 }
829 
830 
831 /**
832  * Return the curvature of the torus.
833  */
834 void
835 rt_tor_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
836 {
837  register struct tor_specific *tor =
838  (struct tor_specific *)stp->st_specific;
839  vect_t w4, w5;
840  fastf_t nx, ny, nz, x_1, y_1, z_1;
841  fastf_t d;
842 
843  nx = tor->tor_N[X];
844  ny = tor->tor_N[Y];
845  nz = tor->tor_N[Z];
846 
847  /* vector from V to hit point */
848  VSUB2(w4, hitp->hit_point, tor->tor_V);
849 
850  if (!NEAR_ZERO(nz, 0.00001)) {
851  z_1 = w4[Z]*nx*nx + w4[Z]*ny*ny - w4[X]*nx*nz - w4[Y]*ny*nz;
852  x_1 = (nx*(z_1-w4[Z])/nz) + w4[X];
853  y_1 = (ny*(z_1-w4[Z])/nz) + w4[Y];
854  } else if (!NEAR_ZERO(ny, 0.00001)) {
855  y_1 = w4[Y]*nx*nx + w4[Y]*nz*nz - w4[X]*nx*ny - w4[Z]*ny*nz;
856  x_1 = (nx*(y_1-w4[Y])/ny) + w4[X];
857  z_1 = (nz*(y_1-w4[Y])/ny) + w4[Z];
858  } else {
859  x_1 = w4[X]*ny*ny + w4[X]*nz*nz - w4[Y]*nx*ny - w4[Z]*nz*nx;
860  y_1 = (ny*(x_1-w4[X])/nx) + w4[Y];
861  z_1 = (nz*(x_1-w4[X])/nx) + w4[Z];
862  }
863  d = sqrt(x_1*x_1 + y_1*y_1 + z_1*z_1);
864 
865  cvp->crv_c1 = (tor->tor_r1 - d) / (d * tor->tor_r2);
866  cvp->crv_c2 = -1.0 / tor->tor_r2;
867 
868  w4[X] = x_1 / d;
869  w4[Y] = y_1 / d;
870  w4[Z] = z_1 / d;
871  VCROSS(w5, tor->tor_N, w4);
872  VCROSS(cvp->crv_pdir, w5, hitp->hit_normal);
873  VUNITIZE(cvp->crv_pdir);
874 }
875 
876 
877 void
878 rt_tor_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
879 {
880  struct tor_specific *tor = (struct tor_specific *) stp->st_specific;
881  vect_t work;
882  vect_t pprime;
883  vect_t pprime2;
884  fastf_t costheta;
885 
886  if (ap) RT_CK_APPLICATION(ap);
887  if (stp) RT_CK_SOLTAB(stp);
888  if (!hitp || !uvp)
889  return;
890  RT_CK_HIT(hitp);
891 
892  VSUB2(work, hitp->hit_point, tor->tor_V);
893  MAT4X3VEC(pprime, tor->tor_SoR, work);
894  /*
895  * -pi/2 <= atan2(x, y) <= pi/2
896  */
897  uvp->uv_u = atan2(pprime[Y], pprime[X]) * M_1_2PI + 0.5;
898 
899  VSET(work, pprime[X], pprime[Y], 0.0);
900  VUNITIZE(work);
901  VSUB2(pprime2, pprime, work);
902  VUNITIZE(pprime2);
903  costheta = VDOT(pprime2, work);
904  uvp->uv_v = atan2(pprime2[Z], costheta) * M_1_2PI + 0.5;
905 }
906 
907 
908 void
909 rt_tor_free(struct soltab *stp)
910 {
911  register struct tor_specific *tor =
912  (struct tor_specific *)stp->st_specific;
913 
914  BU_PUT(tor, struct tor_specific);
915 }
916 
917 
918 /**
919  * Given a circle with a specified radius, determine the minimum
920  * number of straight line segments that the circle can be
921  * approximated with, while still meeting the given maximum
922  * permissible error distance. Form a chord (straight line) by
923  * connecting the start and end points found when sweeping a 'radius'
924  * arc through angle 'theta'.
925  *
926  * The error distance is the distance between where a radius line at
927  * angle theta/2 hits the chord, and where it hits the circle (at
928  * 'radius' distance).
929  *
930  * error_distance = radius * (1 - cos(theta/2))
931  *
932  * or
933  *
934  * theta = 2 * acos(1 - error_distance / radius)
935  *
936  * Returns -
937  * number of segments. Always at least 6.
938  */
939 int
940 rt_num_circular_segments(double maxerr, double radius)
941 {
942  register fastf_t cos_half_theta;
943  register fastf_t half_theta;
944  int n;
945 
946  if (radius <= SMALL_FASTF || maxerr <= SMALL_FASTF || maxerr >= radius) {
947  /* Return a default number of segments */
948  return 6;
949  }
950  cos_half_theta = 1.0 - maxerr / radius;
951 
952  /* There does not seem to be any reasonable way to express the
953  * acos in terms of an atan2(), so extra checking is done. only
954  * proceed if the cosine is between 0 and 1.
955  */
956  if (cos_half_theta <= SMALL_FASTF || cos_half_theta-1.0 >= -SMALL_FASTF) {
957  /* Return a default number of segments */
958  return 6;
959  }
960 
961  half_theta = acos(cos_half_theta);
962 
963  if (half_theta <= SMALL_FASTF) {
964  /* A very large number of segments will be needed. Impose an
965  * upper bound here
966  */
967  return 360*10;
968  }
969  n = (M_PI / half_theta) + 0.99;
970 
971  /* Impose the limits again */
972  if (n <= 6) return 6;
973  if (n >= 360*10) return 360*10;
974  return n;
975 }
976 
977 static int
978 tor_ellipse_points(
979  vect_t ellipse_A,
980  vect_t ellipse_B,
981  const struct rt_view_info *info)
982 {
983  fastf_t avg_radius, circumference;
984 
985  avg_radius = (MAGNITUDE(ellipse_A) + MAGNITUDE(ellipse_B)) / 2.0;
986  circumference = M_2PI * avg_radius;
987 
988  return circumference / info->point_spacing;
989 }
990 
991 int
992 rt_tor_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
993 {
994  vect_t a, b, tor_a, tor_b, tor_h, center;
995  fastf_t mag_a, mag_b, mag_h;
996  struct rt_tor_internal *tor;
997  fastf_t radian, radian_step;
998  int i, num_ellipses, points_per_ellipse;
999 
1000  BU_CK_LIST_HEAD(info->vhead);
1001  RT_CK_DB_INTERNAL(ip);
1002  tor = (struct rt_tor_internal *)ip->idb_ptr;
1003  RT_TOR_CK_MAGIC(tor);
1004 
1005  VMOVE(tor_a, tor->a);
1006  mag_a = tor->r_a;
1007 
1008  VSCALE(tor_h, tor->h, tor->r_h);
1009  mag_h = tor->r_h;
1010 
1011  VCROSS(tor_b, tor_a, tor_h);
1012  VUNITIZE(tor_b);
1013  VSCALE(tor_b, tor_b, mag_a);
1014  mag_b = mag_a;
1015 
1016  /* plot outer circular contour */
1017  VJOIN1(a, tor_a, mag_h / mag_a, tor_a);
1018  VJOIN1(b, tor_b, mag_h / mag_b, tor_b);
1019 
1020  points_per_ellipse = tor_ellipse_points(a, b, info);
1021  if (points_per_ellipse < 6) {
1022  points_per_ellipse = 6;
1023  }
1024 
1025  plot_ellipse(info->vhead, tor->v, a, b, points_per_ellipse);
1026 
1027  /* plot inner circular contour */
1028  VJOIN1(a, tor_a, -1.0 * mag_h / mag_a, tor_a);
1029  VJOIN1(b, tor_b, -1.0 * mag_h / mag_b, tor_b);
1030 
1031  points_per_ellipse = tor_ellipse_points(a, b, info);
1032  if (points_per_ellipse < 6) {
1033  points_per_ellipse = 6;
1034  }
1035 
1036  plot_ellipse(info->vhead, tor->v, a, b, points_per_ellipse);
1037 
1038  /* Draw parallel circles to show the primitive's most extreme points along
1039  * +h/-h.
1040  */
1041  points_per_ellipse = tor_ellipse_points(tor_a, tor_b, info);
1042  if (points_per_ellipse < 6) {
1043  points_per_ellipse = 6;
1044  }
1045 
1046  VADD2(center, tor->v, tor_h);
1047  plot_ellipse(info->vhead, center, tor_a, tor_b, points_per_ellipse);
1048 
1049  VJOIN1(center, tor->v, -1.0, tor_h);
1050  plot_ellipse(info->vhead, center, tor_a, tor_b, points_per_ellipse);
1051 
1052  /* draw circular radial cross sections */
1053  VMOVE(b, tor_h);
1054 
1055  num_ellipses = primitive_curve_count(ip, info);
1056  if (num_ellipses < 3) {
1057  num_ellipses = 3;
1058  }
1059 
1060  radian_step = M_2PI / num_ellipses;
1061  radian = 0;
1062  for (i = 0; i < num_ellipses; ++i) {
1063  ellipse_point_at_radian(center, tor->v, tor_a, tor_b, radian);
1064 
1065  VJOIN1(a, center, -1.0, tor->v);
1066  VUNITIZE(a);
1067  VSCALE(a, a, mag_h);
1068 
1069  plot_ellipse(info->vhead, center, a, b, points_per_ellipse);
1070 
1071  radian += radian_step;
1072  }
1073 
1074  return 0;
1075 }
1076 
1077 /**
1078  * The TORUS has the following input fields:
1079  * ti.v V from origin to center
1080  * ti.h Radius Vector, Normal to plane of torus
1081  * ti.a, ti.b perpendicular, to CENTER of torus (for top, bottom)
1082  */
1083 int
1084 rt_tor_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *UNUSED(tol), const struct rt_view_info *UNUSED(info))
1085 {
1086  fastf_t alpha;
1087  fastf_t beta;
1088  fastf_t cos_alpha, sin_alpha;
1089  fastf_t cos_beta, sin_beta;
1090  fastf_t dist_to_rim;
1091  struct rt_tor_internal *tip;
1092  int w;
1093  int nw = 8;
1094  int len;
1095  int nlen = 16;
1096  fastf_t *pts;
1097  vect_t G;
1098  vect_t radius;
1099  vect_t edge;
1100  fastf_t rel;
1101 
1102  BU_CK_LIST_HEAD(vhead);
1103  RT_CK_DB_INTERNAL(ip);
1104  tip = (struct rt_tor_internal *)ip->idb_ptr;
1105  RT_TOR_CK_MAGIC(tip);
1106 
1107  rel = primitive_get_absolute_tolerance(ttol, 2.0 * (tip->r_a + tip->r_h));
1108 
1109  if (ttol->rel <= 0.0 || ttol->rel >= 1.0) {
1110  rel = 0.0; /* none */
1111  } else {
1112  /* Convert relative tolerance to absolute tolerance
1113  * by scaling w.r.t. the torus diameter.
1114  */
1115  rel = ttol->rel * 2 * (tip->r_a+tip->r_h);
1116  }
1117  /* Take tighter of two (absolute) tolerances */
1118  if (ttol->abs <= 0.0) {
1119  /* No absolute tolerance given */
1120  if (rel <= 0.0) {
1121  /* User has no tolerance for this kind of drink */
1122  nw = 8;
1123  nlen = 16;
1124  } else {
1125  /* Use the absolute-ized relative tolerance */
1126  nlen = rt_num_circular_segments(rel, tip->r_a);
1127  nw = rt_num_circular_segments(rel, tip->r_h);
1128  }
1129  } else {
1130  /* Absolute tolerance was given */
1131  if (rel <= 0.0 || rel > ttol->abs)
1132  rel = ttol->abs;
1133  nlen = rt_num_circular_segments(rel, tip->r_a);
1134  nw = rt_num_circular_segments(rel, tip->r_h);
1135  }
1136 
1137  /* Implement surface-normal tolerance, if given:
1138  *
1139  * nseg = (2 * pi) / (2 * tol)
1140  *
1141  * For a facet which subtends angle theta, surface normal is exact
1142  * in the center, and off by theta/2 at the edges. Note: 1 degree
1143  * tolerance requires 180*180 tessellation!
1144  */
1145  if (ttol->norm > 0.0) {
1146  register int nseg;
1147  nseg = (M_PI / ttol->norm) + 0.99;
1148  if (nseg > nlen) nlen = nseg;
1149  if (nseg > nw) nw = nseg;
1150  }
1151 
1152  /* Compute the points on the surface of the torus */
1153  dist_to_rim = tip->r_h/tip->r_a;
1154  pts = (fastf_t *)bu_malloc(nw * nlen * sizeof(point_t),
1155  "rt_tor_plot pts[]");
1156 
1157 #define TOR_PT(www, lll) ((((www)%nw)*nlen)+((lll)%nlen))
1158 #define TOR_PTA(ww, ll) (&pts[TOR_PT(ww, ll)*3])
1159 #define TOR_NORM_A(ww, ll) (&norms[TOR_PT(ww, ll)*3])
1160 
1161  for (len = 0; len < nlen; len++) {
1162  beta = M_2PI * len / nlen;
1163  cos_beta = cos(beta);
1164  sin_beta = sin(beta);
1165  /* G always points out to rim, along radius vector */
1166  VCOMB2(radius, cos_beta, tip->a, sin_beta, tip->b);
1167  /* We assume that |radius| = |A|. Circular */
1168  VSCALE(G, radius, dist_to_rim);
1169  for (w = 0; w < nw; w++) {
1170  alpha = M_2PI * w / nw;
1171  cos_alpha = cos(alpha);
1172  sin_alpha = sin(alpha);
1173  VCOMB2(edge, cos_alpha, G, sin_alpha*tip->r_h, tip->h);
1174  VADD3(TOR_PTA(w, len), tip->v, edge, radius);
1175  }
1176  }
1177 
1178  /* Draw lengthwise (around outside rim) */
1179  for (w = 0; w < nw; w++) {
1180  len = nlen-1;
1181  RT_ADD_VLIST(vhead, TOR_PTA(w, len), BN_VLIST_LINE_MOVE);
1182  for (len = 0; len < nlen; len++) {
1183  RT_ADD_VLIST(vhead, TOR_PTA(w, len), BN_VLIST_LINE_DRAW);
1184  }
1185  }
1186  /* Draw around the "width" (1 cross section) */
1187  for (len = 0; len < nlen; len++) {
1188  w = nw-1;
1189  RT_ADD_VLIST(vhead, TOR_PTA(w, len), BN_VLIST_LINE_MOVE);
1190  for (w = 0; w < nw; w++) {
1191  RT_ADD_VLIST(vhead, TOR_PTA(w, len), BN_VLIST_LINE_DRAW);
1192  }
1193  }
1194 
1195  bu_free((char *)pts, "rt_tor_plot pts[]");
1196  return 0;
1197 }
1198 
1199 
1200 int
1201 rt_tor_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
1202 {
1203  fastf_t alpha;
1204  fastf_t beta;
1205  fastf_t cos_alpha, sin_alpha;
1206  fastf_t cos_beta, sin_beta;
1207  fastf_t dist_to_rim;
1208  struct rt_tor_internal *tip;
1209  int w;
1210  int nw = 6;
1211  int len;
1212  int nlen = 6;
1213  fastf_t *pts;
1214  vect_t G;
1215  vect_t radius;
1216  vect_t edge;
1217  struct shell *s;
1218  struct vertex **verts;
1219  struct faceuse **faces;
1220  fastf_t *norms;
1221  struct vertex **vertp[4];
1222  int nfaces;
1223  int i;
1224  fastf_t rel;
1225 
1226  RT_CK_DB_INTERNAL(ip);
1227  tip = (struct rt_tor_internal *)ip->idb_ptr;
1228  RT_TOR_CK_MAGIC(tip);
1229 
1230  if (ttol->rel <= 0.0 || ttol->rel >= 1.0) {
1231  rel = 0.0; /* none */
1232  } else {
1233  /* Convert relative tolerance to absolute tolerance by scaling
1234  * w.r.t. the torus diameter.
1235  */
1236  rel = ttol->rel * 2 * (tip->r_a+tip->r_h);
1237  }
1238  /* Take tighter of two (absolute) tolerances */
1239  if (ttol->abs <= 0.0) {
1240  /* No absolute tolerance given */
1241  if (rel <= 0.0) {
1242  /* User has no tolerance for this kind of drink */
1243  nw = 8;
1244  nlen = 16;
1245  } else {
1246  /* Use the absolute-ized relative tolerance */
1247  nlen = rt_num_circular_segments(rel, tip->r_a);
1248  nw = rt_num_circular_segments(rel, tip->r_h);
1249  }
1250  } else {
1251  /* Absolute tolerance was given */
1252  if (rel <= 0.0 || rel > ttol->abs)
1253  rel = ttol->abs;
1254  nlen = rt_num_circular_segments(rel, tip->r_a);
1255  nw = rt_num_circular_segments(rel, tip->r_h);
1256  }
1257 
1258  /* Implement surface-normal tolerance, if given
1259  *
1260  * nseg = (2 * pi) / (2 * tol)
1261  *
1262  * For a facet which subtends angle theta, surface normal is exact
1263  * in the center, and off by theta/2 at the edges. Note: 1 degree
1264  * tolerance requires 180*180 tessellation!
1265  */
1266  if (ttol->norm > 0.0) {
1267  register int nseg;
1268  nseg = (M_PI / ttol->norm) + 0.99;
1269  if (nseg > nlen) nlen = nseg;
1270  if (nseg > nw) nw = nseg;
1271  }
1272 
1273  /* Compute the points on the surface of the torus */
1274  dist_to_rim = tip->r_h/tip->r_a;
1275  pts = (fastf_t *)bu_malloc(nw * nlen * sizeof(point_t),
1276  "rt_tor_tess pts[]");
1277  norms = (fastf_t *)bu_malloc(nw * nlen * sizeof(vect_t), "rt_tor_tess: norms[]");
1278 
1279  for (len = 0; len < nlen; len++) {
1280  beta = M_2PI * len / nlen;
1281  cos_beta = cos(beta);
1282  sin_beta = sin(beta);
1283  /* G always points out to rim, along radius vector */
1284  VCOMB2(radius, cos_beta, tip->a, sin_beta, tip->b);
1285  /* We assume that |radius| = |A|. Circular */
1286  VSCALE(G, radius, dist_to_rim);
1287  for (w = 0; w < nw; w++) {
1288  alpha = M_2PI * w / nw;
1289  cos_alpha = cos(alpha);
1290  sin_alpha = sin(alpha);
1291  VCOMB2(edge, cos_alpha, G, sin_alpha*tip->r_h, tip->h);
1292  VADD3(TOR_PTA(w, len), tip->v, edge, radius);
1293 
1294  VMOVE(TOR_NORM_A(w, len), edge);
1295  VUNITIZE(TOR_NORM_A(w, len));
1296  }
1297  }
1298 
1299  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
1300  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
1301 
1302  verts = (struct vertex **)bu_calloc(nw*nlen, sizeof(struct vertex *),
1303  "rt_tor_tess *verts[]");
1304  faces = (struct faceuse **)bu_calloc(nw*nlen, sizeof(struct faceuse *),
1305  "rt_tor_tess *faces[]");
1306 
1307  /* Build the topology of the torus */
1308  /* Note that increasing 'w' goes to the left (alpha goes ccw) */
1309  nfaces = 0;
1310  for (w = 0; w < nw; w++) {
1311  for (len = 0; len < nlen; len++) {
1312  vertp[0] = &verts[ TOR_PT(w+0, len+0) ];
1313  vertp[1] = &verts[ TOR_PT(w+0, len+1) ];
1314  vertp[2] = &verts[ TOR_PT(w+1, len+1) ];
1315  vertp[3] = &verts[ TOR_PT(w+1, len+0) ];
1316  if ((faces[nfaces++] = nmg_cmface(s, vertp, 4)) == (struct faceuse *)0) {
1317  bu_log("rt_tor_tess() nmg_cmface failed, w=%d/%d, len=%d/%d\n",
1318  w, nw, len, nlen);
1319  nfaces--;
1320  }
1321  }
1322  }
1323 
1324  /* Associate vertex geometry */
1325  for (w = 0; w < nw; w++) {
1326  for (len = 0; len < nlen; len++) {
1327  nmg_vertex_gv(verts[TOR_PT(w, len)], TOR_PTA(w, len));
1328  }
1329  }
1330 
1331  /* Associate face geometry */
1332  for (i=0; i < nfaces; i++) {
1333  if (nmg_fu_planeeqn(faces[i], tol) < 0)
1334  return -1; /* FAIL */
1335  }
1336 
1337  /* Associate vertexuse normals */
1338  for (w=0; w<nw; w++) {
1339  for (len=0; len<nlen; len++) {
1340  struct vertexuse *vu;
1341  vect_t rev_norm;
1342 
1343  VREVERSE(rev_norm, TOR_NORM_A(w, len));
1344 
1345  for (BU_LIST_FOR(vu, vertexuse, &verts[TOR_PT(w, len)]->vu_hd)) {
1346  struct faceuse *fu;
1347 
1348  NMG_CK_VERTEXUSE(vu);
1349 
1350  fu = nmg_find_fu_of_vu(vu);
1351  NMG_CK_FACEUSE(fu);
1352 
1353  if (fu->orientation == OT_SAME)
1354  nmg_vertexuse_nv(vu, TOR_NORM_A(w, len));
1355  else if (fu->orientation == OT_OPPOSITE)
1356  nmg_vertexuse_nv(vu, rev_norm);
1357  }
1358  }
1359  }
1360 
1361  /* kill zero length edgeuse */
1362  (void)nmg_keu_zl(s, tol);
1363 
1364  /* Compute "geometry" for region and shell */
1365  nmg_region_a(*r, tol);
1366 
1367  bu_free((char *)pts, "rt_tor_tess pts[]");
1368  bu_free((char *)verts, "rt_tor_tess *verts[]");
1369  bu_free((char *)faces, "rt_tor_tess *faces[]");
1370  bu_free((char *)norms, "rt_tor_tess norms[]");
1371  return 0;
1372 }
1373 
1374 
1375 /**
1376  * Import a torus from the database format to the internal format.
1377  * Apply modeling transformations at the same time.
1378  */
1379 int
1380 rt_tor_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
1381 {
1382  struct rt_tor_internal *tip;
1383  union record *rp;
1384  fastf_t vec[3*4];
1385  vect_t axb;
1386  register fastf_t f;
1387 
1388  if (dbip) RT_CK_DBI(dbip);
1389 
1390  BU_CK_EXTERNAL(ep);
1391  rp = (union record *)ep->ext_buf;
1392  /* Check record type */
1393  if (rp->u_id != ID_SOLID) {
1394  bu_log("rt_tor_import4: defective record\n");
1395  return -1;
1396  }
1397 
1398  RT_CK_DB_INTERNAL(ip);
1399  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1400  ip->idb_type = ID_TOR;
1401  ip->idb_meth = &OBJ[ID_TOR];
1402  BU_ALLOC(ip->idb_ptr, struct rt_tor_internal);
1403 
1404  tip = (struct rt_tor_internal *)ip->idb_ptr;
1405  tip->magic = RT_TOR_INTERNAL_MAGIC;
1406 
1407  /* Convert from database to internal format */
1408  flip_fastf_float(vec, rp->s.s_values, 4, dbip->dbi_version < 0 ? 1 : 0);
1409 
1410  /* Apply modeling transformations */
1411  if (mat == NULL) mat = bn_mat_identity;
1412  MAT4X3PNT(tip->v, mat, &vec[0*3]);
1413  MAT4X3VEC(tip->h, mat, &vec[1*3]);
1414  MAT4X3VEC(tip->a, mat, &vec[2*3]);
1415  MAT4X3VEC(tip->b, mat, &vec[3*3]);
1416 
1417  /* Make the vectors unit length */
1418  tip->r_a = MAGNITUDE(tip->a);
1419  tip->r_b = MAGNITUDE(tip->b);
1420  tip->r_h = MAGNITUDE(tip->h);
1421  if (tip->r_a <= SMALL_FASTF || tip->r_b <= SMALL_FASTF || tip->r_h <= SMALL_FASTF) {
1422  bu_log("rt_tor_import4: zero length A, B, or H vector\n");
1423  return -1;
1424  }
1425  /* In memory, the H vector is unit length */
1426  f = 1.0/tip->r_h;
1427  VSCALE(tip->h, tip->h, f);
1428 
1429  /* If H does not point in the direction of A cross B, reverse H. */
1430  /* Somehow, database records have been written with this problem. */
1431  VCROSS(axb, tip->a, tip->b);
1432  if (VDOT(axb, tip->h) < 0) {
1433  VREVERSE(tip->h, tip->h);
1434  }
1435 
1436  return 0; /* OK */
1437 }
1438 
1439 
1440 int
1441 rt_tor_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1442 {
1443  struct rt_tor_internal *tip;
1444 
1445  /* must be double for export */
1446  double vec[2*3+2];
1447 
1448  if (dbip) RT_CK_DBI(dbip);
1449 
1450  RT_CK_DB_INTERNAL(ip);
1451  if (ip->idb_type != ID_TOR) return -1;
1452  tip = (struct rt_tor_internal *)ip->idb_ptr;
1453  RT_TOR_CK_MAGIC(tip);
1454 
1455  BU_CK_EXTERNAL(ep);
1456  ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * (2*3+2);
1457  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "tor external");
1458 
1459  /* scale values into local buffer */
1460  VSCALE(&vec[0*3], tip->v, local2mm);
1461  VMOVE(&vec[1*3], tip->h); /* UNIT vector, not scaled */
1462  vec[2*3+0] = tip->r_a*local2mm; /* r1 */
1463  vec[2*3+1] = tip->r_h*local2mm; /* r2 */
1464 
1465  /* convert from internal (host) to database (network) format */
1466  bu_cv_htond(ep->ext_buf, (unsigned char *)vec, 2*3+2);
1467 
1468  return 0;
1469 }
1470 
1471 
1472 /**
1473  * The name will be added by the caller.
1474  */
1475 int
1476 rt_tor_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1477 {
1478  struct rt_tor_internal *tip;
1479  union record *rec;
1480  vect_t norm;
1481  vect_t cross1, cross2;
1482  fastf_t r1, r2;
1483  fastf_t r3, r4;
1484  double m2;
1485 
1486  if (dbip) RT_CK_DBI(dbip);
1487 
1488  RT_CK_DB_INTERNAL(ip);
1489  if (ip->idb_type != ID_TOR) return -1;
1490  tip = (struct rt_tor_internal *)ip->idb_ptr;
1491  RT_TOR_CK_MAGIC(tip);
1492 
1493  BU_CK_EXTERNAL(ep);
1494  ep->ext_nbytes = sizeof(union record);
1495  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "tor external");
1496  rec = (union record *)ep->ext_buf;
1497 
1498  rec->s.s_id = ID_SOLID;
1499  rec->s.s_type = TOR;
1500 
1501  r1 = tip->r_a;
1502  r2 = tip->r_h;
1503 
1504  /* Validate that 0 < r2 <= r1 */
1505  if (r2 <= 0.0) {
1506  bu_log("rt_tor_export4: illegal r2=%.12e <= 0\n", r2);
1507  return -1;
1508  }
1509  if (r2 > r1) {
1510  bu_log("rt_tor_export4: illegal r2=%.12e > r1=%.12e\n",
1511  r2, r1);
1512  return -1;
1513  }
1514 
1515  r1 *= local2mm;
1516  r2 *= local2mm;
1517  VSCALE(&rec->s.s_values[0*3], tip->v, local2mm);
1518 
1519  VMOVE(norm, tip->h);
1520  m2 = MAGNITUDE(norm); /* F2 is NORMAL to torus */
1521  if (m2 <= SQRT_SMALL_FASTF) {
1522  bu_log("rt_tor_export4: normal magnitude is zero!\n");
1523  return -1; /* failure */
1524  }
1525  m2 = 1.0/m2;
1526  VSCALE(norm, norm, m2); /* Give normal unit length */
1527  VSCALE(&rec->s.s_values[1*3], norm, r2); /* F2: normal radius len */
1528 
1529  /* Create two mutually perpendicular vectors, perpendicular to Norm */
1530  /* Ensure that AxB points in direction of N */
1531  bn_vec_ortho(cross1, norm);
1532  VCROSS(cross2, norm, cross1);
1533  VUNITIZE(cross2);
1534 
1535  /* F3, F4 are perpendicular, goto center of solid part */
1536  VSCALE(&rec->s.s_values[2*3], cross1, r1);
1537  VSCALE(&rec->s.s_values[3*3], cross2, r1);
1538 
1539  /*
1540  * The rest of these provide no real extra information,
1541  * and exist for compatibility with old versions of MGED.
1542  */
1543  r3=r1-r2; /* Radius to inner circular edge */
1544  r4=r1+r2; /* Radius to outer circular edge */
1545 
1546  /* F5, F6 are perpendicular, goto inner edge of ellipse */
1547  VSCALE(&rec->s.s_values[4*3], cross1, r3);
1548  VSCALE(&rec->s.s_values[5*3], cross2, r3);
1549 
1550  /* F7, F8 are perpendicular, goto outer edge of ellipse */
1551  VSCALE(&rec->s.s_values[6*3], cross1, r4);
1552  VSCALE(&rec->s.s_values[7*3], cross2, r4);
1553 
1554  return 0;
1555 }
1556 
1557 
1558 /**
1559  * Taken from the database record:
1560  * v vertex (point) of center of torus.
1561  * h unit vector in the normal direction of the torus
1562  * major radius of ring from 'v' to center of ring
1563  * minor radius of the ring
1564  *
1565  * Calculate:
1566  * 2nd radius of ring (==1st radius)
1567  * ring unit vector 1
1568  * ring unit vector 2
1569  */
1570 int
1571 rt_tor_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
1572 {
1573  struct rt_tor_internal *tip;
1574 
1575  /* must be double for import and export */
1576  struct rec {
1577  double v[3];
1578  double h[3];
1579  double ra; /* r1 */
1580  double rh; /* r2 */
1581  } rec;
1582 
1583  if (dbip) RT_CK_DBI(dbip);
1584 
1585  BU_CK_EXTERNAL(ep);
1586  BU_ASSERT_LONG(ep->ext_nbytes, ==, SIZEOF_NETWORK_DOUBLE * (2*3+2));
1587 
1588  RT_CK_DB_INTERNAL(ip);
1589 
1590  if (mat == NULL) mat = bn_mat_identity;
1591  if (bn_mat_is_non_unif (mat)) {
1592  bu_log("------------------ WARNING ----------------\nNon-uniform matrix transform on torus. Ignored\n");
1593  }
1594 
1595  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1596  ip->idb_type = ID_TOR;
1597  ip->idb_meth = &OBJ[ID_TOR];
1598  BU_ALLOC(ip->idb_ptr, struct rt_tor_internal);
1599 
1600  tip = (struct rt_tor_internal *)ip->idb_ptr;
1601  tip->magic = RT_TOR_INTERNAL_MAGIC;
1602 
1603  bu_cv_ntohd((unsigned char *)&rec, ep->ext_buf, 2*3+2);
1604 
1605  /* Apply modeling transformations */
1606  MAT4X3PNT(tip->v, mat, rec.v);
1607  MAT4X3VEC(tip->h, mat, rec.h);
1608  VUNITIZE(tip->h); /* just to be sure */
1609 
1610  tip->r_a = rec.ra / mat[15];
1611  tip->r_h = rec.rh / mat[15];
1612 
1613  /* Prepare the extra information */
1614  tip->r_b = tip->r_a;
1615 
1616  /* Calculate two mutually perpendicular vectors, perpendicular to N */
1617  bn_vec_ortho(tip->a, tip->h); /* a has unit length */
1618  VCROSS(tip->b, tip->h, tip->a); /* |A| = |H| = 1, so |B|=1 */
1619 
1620  VSCALE(tip->a, tip->a, tip->r_a);
1621  VSCALE(tip->b, tip->b, tip->r_b);
1622  return 0;
1623 }
1624 
1625 
1626 /**
1627  * Make human-readable formatted presentation of this solid. First
1628  * line describes type of solid. Additional lines are indented one
1629  * tab, and give parameter values.
1630  */
1631 int
1632 rt_tor_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1633 {
1634  register struct rt_tor_internal *tip =
1635  (struct rt_tor_internal *)ip->idb_ptr;
1636  double r3, r4;
1637 
1638  RT_TOR_CK_MAGIC(tip);
1639  bu_vls_strcat(str, "torus (TOR)\n");
1640 
1641  bu_vls_printf(str, "\tV (%g, %g, %g), r1=%g (A), r2=%g (H)\n",
1642  INTCLAMP(tip->v[X] * mm2local),
1643  INTCLAMP(tip->v[Y] * mm2local),
1644  INTCLAMP(tip->v[Z] * mm2local),
1645  INTCLAMP(tip->r_a * mm2local),
1646  INTCLAMP(tip->r_h * mm2local));
1647 
1648  bu_vls_printf(str, "\tN=(%g, %g, %g)\n",
1649  INTCLAMP(tip->h[X] * mm2local),
1650  INTCLAMP(tip->h[Y] * mm2local),
1651  INTCLAMP(tip->h[Z] * mm2local));
1652 
1653  if (!verbose) return 0;
1654 
1655  bu_vls_printf(str, "\tA=(%g, %g, %g)\n",
1656  INTCLAMP(tip->a[X] * mm2local / tip->r_a),
1657  INTCLAMP(tip->a[Y] * mm2local / tip->r_a),
1658  INTCLAMP(tip->a[Z] * mm2local / tip->r_a));
1659 
1660  bu_vls_printf(str, "\tB=(%g, %g, %g)\n",
1661  INTCLAMP(tip->b[X] * mm2local / tip->r_b),
1662  INTCLAMP(tip->b[Y] * mm2local / tip->r_b),
1663  INTCLAMP(tip->b[Z] * mm2local / tip->r_b));
1664 
1665  r3 = tip->r_a - tip->r_h;
1666  bu_vls_printf(str, "\tvector to inner edge = (%g, %g, %g)\n",
1667  INTCLAMP(tip->a[X] * mm2local / tip->r_a * r3),
1668  INTCLAMP(tip->a[Y] * mm2local / tip->r_a * r3),
1669  INTCLAMP(tip->a[Z] * mm2local / tip->r_a * r3));
1670 
1671  r4 = tip->r_a + tip->r_h;
1672  bu_vls_printf(str, "\tvector to outer edge = (%g, %g, %g)\n",
1673  INTCLAMP(tip->a[X] * mm2local / tip->r_a * r4),
1674  INTCLAMP(tip->a[Y] * mm2local / tip->r_a * r4),
1675  INTCLAMP(tip->a[Z] * mm2local / tip->r_a * r4));
1676 
1677  return 0;
1678 }
1679 
1680 
1681 /**
1682  * Free the storage associated with the rt_db_internal version of this
1683  * solid.
1684  */
1685 void
1687 {
1688  register struct rt_tor_internal *tip;
1689 
1690  RT_CK_DB_INTERNAL(ip);
1691 
1692  tip = (struct rt_tor_internal *)ip->idb_ptr;
1693  RT_TOR_CK_MAGIC(tip);
1694 
1695  bu_free((char *)tip, "rt_tor_internal");
1696  ip->idb_ptr = ((void *)0); /* sanity */
1697 }
1698 
1699 
1700 int
1701 rt_tor_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
1702 {
1703  if (ip) RT_CK_DB_INTERNAL(ip);
1704 
1705  return 0; /* OK */
1706 }
1707 
1708 
1709 void
1710 rt_tor_surf_area(fastf_t *area, const struct rt_db_internal *ip)
1711 {
1712  struct rt_tor_internal *tip = (struct rt_tor_internal *)ip->idb_ptr;
1713  RT_TOR_CK_MAGIC(tip);
1714  /* r_h: radius of torus tube
1715  * r_a: radius from axis of rotation to center of tube
1716  */
1717  *area = 4.0 * M_PI * M_PI * tip->r_h * tip->r_a;
1718 }
1719 
1720 
1721 void
1722 rt_tor_volume(fastf_t *vol, const struct rt_db_internal *ip)
1723 {
1724  struct rt_tor_internal *tip = (struct rt_tor_internal *)ip->idb_ptr;
1725  RT_TOR_CK_MAGIC(tip);
1726  *vol = 2.0 * M_PI * M_PI * (tip->r_h * tip->r_h) * tip->r_a;
1727 }
1728 
1729 
1730 void
1731 rt_tor_centroid(point_t *cent, const struct rt_db_internal *ip)
1732 {
1733  struct rt_tor_internal *tip = (struct rt_tor_internal *)ip->idb_ptr;
1734  RT_TOR_CK_MAGIC(tip);
1735  VMOVE(*cent,tip->v);
1736 }
1737 
1738 /*
1739  * Local Variables:
1740  * mode: C
1741  * tab-width: 8
1742  * indent-tabs-mode: t
1743  * c-file-style: "stroustrup"
1744  * End:
1745  * ex: shiftwidth=4 tabstop=8
1746  */
Definition: db_flip.c:35
char * d_namep
pointer to name string
Definition: raytrace.h:859
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
Definition: raytrace.h:800
int rt_tor_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: tor.c:1701
fastf_t cf[BN_MAX_POLY_DEGREE+1]
Definition: poly.h:50
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
#define SIZEOF_NETWORK_DOUBLE
Definition: cv.h:48
mat_t tor_invR
Definition: tor.c:154
struct hit seg_in
IN information.
Definition: raytrace.h:370
int(* ft_bbox)(struct rt_db_internal *, point_t *, point_t *, const struct bn_tol *)
Definition: raytrace.h:2196
fastf_t C[2 *MAX_CNT+1][2 *MAX_CNT+1]
Definition: dsp_brep.cpp:38
const struct bu_structparse rt_tor_parse[]
Definition: tor.c:138
Definition: list.h:118
int nmg_fu_planeeqn(struct faceuse *fu, const struct bn_tol *tol)
Definition: nmg_mod.c:1311
const struct directory * st_dp
Directory entry of solid.
Definition: raytrace.h:436
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
double dist
>= 0
Definition: tol.h:73
#define RT_TOR_SEG_MISS(SEG)
Definition: tor.c:534
vect_t crv_pdir
Principle direction.
Definition: raytrace.h:307
const mat_t bn_mat_identity
Matrix and vector functionality.
Definition: mat.c:46
fastf_t uv_u
Range 0..1.
Definition: raytrace.h:341
#define TOR_NORM_A(ww, ll)
if lu s
Definition: nmg_mod.c:3860
void rt_tor_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
Definition: tor.c:835
void bu_vls_strcat(struct bu_vls *vp, const char *s)
Definition: vls.c:368
#define VSET(a, b, c, d)
Definition: color.c:53
int nmg_keu_zl(struct shell *s, const struct bn_tol *tol)
Definition: nmg_mk.c:3089
Definition: raytrace.h:215
void bn_pr_roots(const char *title, const struct bn_complex roots[], int n)
fastf_t tor_alpha
Definition: tor.c:148
#define M_PI
Definition: fft.h:35
Definition: pc.h:108
Definition: raytrace.h:368
#define BU_ASSERT_LONG(_lhs, _relation, _rhs)
Definition: defines.h:240
mat_t tor_SoR
Definition: tor.c:153
Definition: raytrace.h:248
void nmg_vertex_gv(struct vertex *v, const fastf_t *pt)
Definition: nmg_mk.c:1668
struct bn_complex bn_complex_t
Complex numbers.
#define SMALL_FASTF
Definition: defines.h:342
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
int rt_tor_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
Definition: tor.c:345
Header file for the BRL-CAD common definitions.
void nmg_vertexuse_nv(struct vertexuse *vu, const fastf_t *norm)
Definition: nmg_mk.c:1719
Definition: poly.h:47
#define ID_TOR
Toroid.
Definition: raytrace.h:459
void flip_fastf_float(fastf_t *ff, const dbfloat_t *fp, int n, int flip)
Definition: db_flip.c:74
void rt_tor_print(register const struct soltab *stp)
Definition: tor.c:301
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
void bu_cv_htond(unsigned char *out, const unsigned char *in, size_t count)
void rt_pt_sort(fastf_t t[], int npts)
Definition: tgc.c:1394
struct bu_list l
Definition: raytrace.h:369
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
void rt_tor_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: tor.c:1710
double rel
rel dist tol
Definition: raytrace.h:181
int rt_tor_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: tor.c:1380
if(share_geom)
Definition: nmg_mod.c:3829
int idb_major_type
Definition: raytrace.h:192
size_t dgr
Definition: poly.h:49
fastf_t tor_r2
Definition: tor.c:150
Definition: color.c:49
struct rt_i * a_rt_i
this librt instance
Definition: raytrace.h:1588
int bn_mat_is_non_unif(const mat_t m)
Definition: mat.c:1146
vect_t hit_vpriv
PRIVATE vector for xxx_*()
Definition: raytrace.h:253
struct bu_list * vhead
Definition: raytrace.h:1926
#define RT_ADD_VLIST(hd, pnt, draw)
Definition: raytrace.h:1865
void bn_mat_print(const char *title, const mat_t m)
Definition: mat.c:81
#define TOR_PTA(ww, ll)
int rt_tor_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: tor.c:1441
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
#define RT_CK_HIT(_p)
Definition: raytrace.h:259
#define RT_TOR_INTERNAL_MAGIC
Definition: magic.h:112
int rt_tor_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol, const struct rt_view_info *info)
Definition: tor.c:1084
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
void rt_tor_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
Definition: tor.c:878
void bn_mat_inv(mat_t output, const mat_t input)
fastf_t primitive_get_absolute_tolerance(const struct rt_tess_tol *ttol, fastf_t rel_to_abs)
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
uint8_t * ext_buf
Definition: parse.h:216
#define BU_GET(_ptr, _type)
Definition: malloc.h:201
point_t hit_point
DEPRECATED: Intersection point, use VJOIN1 hit_dist.
Definition: raytrace.h:251
point_t st_max
max X, Y, Z of bounding RPP
Definition: raytrace.h:438
#define SQRT_SMALL_FASTF
Definition: defines.h:346
int rt_tor_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: tor.c:220
struct hit seg_out
OUT information.
Definition: raytrace.h:371
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
#define UNUSED(parameter)
Definition: common.h:239
void rt_tor_ifree(struct rt_db_internal *ip)
Definition: tor.c:1686
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
struct faceuse * nmg_find_fu_of_vu(const struct vertexuse *vu)
Definition: nmg_info.c:304
Support for uniform tolerances.
Definition: tol.h:71
int rt_tor_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: tor.c:1201
int rt_tor_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: tor.c:161
void rt_tor_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
Definition: tor.c:808
#define bu_offsetofarray(_t, _a, _d, _i)
Definition: parse.h:65
struct nmgregion * nmg_mrsv(struct model *m)
Definition: nmg_mk.c:306
#define BU_STRUCTPARSE_FUNC_NULL
Definition: parse.h:153
ustring alpha
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
fastf_t primitive_curve_count(struct rt_db_internal *ip, const struct rt_view_info *info)
#define bu_offsetof(_t, _m)
Definition: parse.h:64
#define RT_CK_DBI(_p)
Definition: raytrace.h:829
int rt_tor_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: tor.c:1632
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
fastf_t point_spacing
Definition: raytrace.h:1934
void * idb_ptr
Definition: raytrace.h:195
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
vect_t tor_N
Definition: tor.c:152
point_t st_min
min X, Y, Z of bounding RPP
Definition: raytrace.h:437
void bu_cv_ntohd(unsigned char *out, const unsigned char *in, size_t count)
const struct rt_functab OBJ[]
Definition: table.c:159
int rt_poly_roots(bn_poly_t *eqn, bn_complex_t roots[], const char *name)
void rt_tor_volume(fastf_t *vol, const struct rt_db_internal *ip)
Definition: tor.c:1722
void bn_vec_ortho(vect_t out, const vect_t in)
#define R
Definition: msr.c:55
int rt_tor_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: tor.c:1571
double abs
absolute dist tol
Definition: raytrace.h:180
#define RT_CK_SOLTAB(_p)
Definition: raytrace.h:453
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
#define A
Definition: msr.c:51
Definition: color.c:51
int dbi_version
PRIVATE: use db_version()
Definition: raytrace.h:824
void plot_ellipse(struct bu_list *vhead, const vect_t t, const vect_t a, const vect_t b, int num_points)
double norm
normal tol
Definition: raytrace.h:182
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
vect_t hit_normal
DEPRECATED: Surface Normal at hit_point, use RT_HIT_NORMAL.
Definition: raytrace.h:252
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
int rt_num_circular_segments(double maxerr, double radius)
Definition: tor.c:940
struct faceuse * nmg_cmface(struct shell *s, struct vertex ***verts, int n)
Definition: nmg_mod.c:979
Complex numbers.
Definition: complex.h:39
const struct rt_functab * st_meth
pointer to per-solid methods
Definition: raytrace.h:428
int rt_tor_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
Definition: tor.c:992
size_t ext_nbytes
Definition: parse.h:210
void ellipse_point_at_radian(point_t result, const vect_t center, const vect_t axis_a, const vect_t axis_b, fastf_t radian)
fastf_t tor_r1
Definition: tor.c:149
void rt_tor_centroid(point_t *cent, const struct rt_db_internal *ip)
Definition: tor.c:1731
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
HIDDEN void verbose(struct human_data_t *dude)
Definition: human.c:2008
Definition: vls.h:56
int rt_tor_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: tor.c:1476
double fastf_t
Definition: defines.h:300
#define TOR_PT(www, lll)
#define VPRINT(a, b)
Definition: raytrace.h:1881
void rt_tor_free(struct soltab *stp)
Definition: tor.c:909
vect_t tor_V
Definition: tor.c:151
Definition: color.c:50
#define BU_LIST_FIRST(structure, hp)
Definition: list.h:312
fastf_t uv_v
Range 0..1.
Definition: raytrace.h:342
void rt_tor_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
Definition: tor.c:539
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
void nmg_region_a(struct nmgregion *r, const struct bn_tol *tol)
Definition: nmg_mk.c:2557