BRL-CAD
part.c
Go to the documentation of this file.
1 /* P A R T . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1990-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/part/part.c
23  *
24  * Intersect a ray with a "particle" solid, which can have three main
25  * forms: sphere, hemisphere-tipped cylinder (lozenge), and
26  * hemisphere-tipped cone. This code draws on the examples of g_rec
27  * (Davisson) & g_sph (Dykstra).
28  *
29  * Algorithm for the hemisphere-tipped cylinder and cone cases -
30  *
31  * Given V, H, vrad, and hrad, there is a set of points on this
32  * cylinder
33  *
34  * { (x, y, z) | (x, y, z) is on cylinder }
35  *
36  * Through a series of Affine Transformations, this set of points will
37  * be transformed into a set of points on a unit cylinder (or cone)
38  * with the transformed base (V') located at the origin with a
39  * transformed radius of 1 (vrad'). The height of the cylinder (or
40  * cone) along the +Z axis is +1 (i.e., H' = (0, 0, 1)), with a
41  * transformed radius of hrad/vrad.
42  *
43  *
44  * { (x', y', z') | (x', y', z') is on cylinder at origin }
45  *
46  * The transformation from X to X' is accomplished by:
47  *
48  * finding two unit vectors A and B mutually perpendicular, and
49  * perp. to H.
50  *
51  * X' = S(R(X - V))
52  *
53  * where R(X) rotates H to the +Z axis, and S(X) scales vrad' to 1 and
54  * |H'| to 1.
55  *
56  * where R(X) = (A/(|A|))
57  * (B/(|B|)) . X
58  * (H/(|H|))
59  *
60  * and S(X) = (1/|A| 0 0)
61  * (0 1/|B| 0) . X
62  * (0 0 1/|H|)
63  *
64  * To find the intersection of a line with the surface of the cylinder,
65  * consider the parametric line L:
66  *
67  * L : { P(n) | P + t(n) . D }
68  *
69  * Call W the actual point of intersection between L and the cylinder.
70  * Let W' be the point of intersection between L' and the unit cylinder.
71  *
72  * L' : { P'(n) | P' + t(n) . D' }
73  *
74  * W = invR(invS(W')) + V
75  *
76  * Where W' = k D' + P'.
77  *
78  * If Dx' and Dy' are both 0, then there is no hit on the cylinder;
79  * but the end spheres need checking.
80  *
81  * The equation for the unit cylinder ranging along Z is
82  *
83  * x**2 + y**2 - r**2 = 0
84  *
85  * and the equation for a unit cone ranging along Z is
86  *
87  * x**2 + y**2 - f(z)**2 = 0
88  *
89  * where in this case f(z) linearly interpolates the radius of the
90  * cylinder from vrad (r1) to hrad (r2) as z ranges from 0 to 1, i.e.:
91  *
92  * f(z) = (r2-r1)/1 * z + r1
93  *
94  * let m = (r2-r1)/1, and substitute:
95  *
96  * x**2 + y**2 - (m*z+r1)**2 = 0 .
97  *
98  * For the cylinder case, r1 == r2, so m == 0, and everything
99  * simplifies.
100  *
101  * The parametric formulation for line L' is P' + t * D', or
102  *
103  * x = Px' + t * Dx'
104  * y = Py' + t * Dy'
105  * z = Pz' + t * Dz' .
106  *
107  * Substituting these definitions into the formula for the unit cone gives
108  *
109  * (Px'+t*Dx')**2 + (Py'+t*Dy')**2 + (m*(Pz'+t*Dz')+r1)**2 = 0
110  *
111  * Expanding and regrouping terms gives a quadratic in "t" which has
112  * the form
113  *
114  * a * t**2 + b * t + c = 0
115  *
116  * where
117  *
118  * a = Dx'**2 + Dy'**2 - m**2 * Dz'**2
119  * b = 2 * (Px'*Dx' + Py'*Dy' - m**2 * Pz'*Dz' - m*r1*Dz')
120  * c = Px'**2 + Py'**2 - m**2 * Pz'**2 - 2*m*r1*Pz' - r1**2
121  *
122  * Line L' hits the infinitely tall unit cone at point(s) W' which
123  * correspond to the roots of the quadratic. The quadratic formula
124  * yields values for "t"
125  *
126  * t = [ -b +/- sqrt(b** - 4 * a * c) ] / (2 * a)
127  *
128  * This parameter "t" can be substituted into the formulas for either
129  * L' or L, because affine transformations preserve distances along
130  * lines.
131  *
132  * Now, D' = S(R(D))
133  * and P' = S(R(P - V))
134  *
135  * Substituting,
136  *
137  * W = V + invR(invS[ k *(S(R(D))) + S(R(P - V)) ])
138  * = V + invR((k * R(D)) + R(P - V))
139  * = V + k * D + P - V
140  * = k * D + P
141  *
142  * Note that ``t'' is constant, and is the same in the formulations
143  * for both W and W'.
144  *
145  * The hit at ``t'' is a hit on the height=1 unit cylinder IFF
146  * 0 <= Wz' <= 1.
147  *
148  * NORMALS. Given the point W on the surface of the cylinder, what is
149  * the vector normal to the tangent plane at that point?
150  *
151  * Map W onto the unit cylinder, i.e.: W' = S(R(W - V)).
152  *
153  * Plane on unit cylinder at W' has a normal vector N' of the same
154  * value as W' in x and y, with z set to zero, i.e., (Wx', Wy', 0)
155  *
156  * The plane transforms back to the tangent plane at W, and this new
157  * plane (on the original cylinder) has a normal vector of N, viz:
158  *
159  * N = inverse[ transpose(invR o invS) ] (N')
160  * = inverse[ transpose(invS) o transpose(invR) ] (N')
161  * = inverse[ inverse(S) o R ] (N')
162  * = invR o S (N')
163  *
164  * Note that the normal vector produced above will not have unit length.
165  *
166  * THE HEMISPHERES.
167  *
168  * THE "EQUIVALENT CONE":
169  *
170  * In order to have exact matching of the surface normals at the join
171  * between the conical body of the particle and the hemispherical end,
172  * it is necessary to alter the cone to form an "equivalent cone",
173  * where the end caps of the cone are both shifted away from the large
174  * hemisphere and towards the smaller one. This makes the cone end
175  * where it is tangent to the hemisphere. The calculation for theta
176  * come from a diagram drawn by PJT on 18-Nov-99.
177  */
178 /** @} */
179 
180 #include "common.h"
181 
182 #include <stddef.h>
183 #include <string.h>
184 #include <math.h>
185 #include "bio.h"
186 
187 #include "bu/cv.h"
188 #include "vmath.h"
189 #include "db.h"
190 #include "rtgeom.h"
191 #include "raytrace.h"
192 #include "nmg.h"
193 #include "../../librt_private.h"
194 
195 
197  struct rt_part_internal part_int;
198  mat_t part_SoR; /* Scale(Rot(vect)) */
199  mat_t part_invRoS; /* invRot(Scale(vect)) */
202  /* For the "equivalent cone" */
203  fastf_t part_v_hdist; /* dist of base plate on unit cone */
205  fastf_t part_v_erad; /* radius of equiv. particle */
207 };
208 
209 
210 /* hit_surfno flags for which end was hit */
211 #define RT_PARTICLE_SURF_VSPHERE 1
212 #define RT_PARTICLE_SURF_BODY 2
213 #define RT_PARTICLE_SURF_HSPHERE 3
214 
215 const struct bu_structparse rt_part_parse[] = {
216  { "%f", 3, "V", bu_offsetofarray(struct rt_part_internal, part_V, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
217  { "%f", 3, "H", bu_offsetofarray(struct rt_part_internal, part_H, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
218  { "%f", 1, "r_v", bu_offsetof(struct rt_part_internal, part_vrad), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
219  { "%f", 1, "r_h", bu_offsetof(struct rt_part_internal, part_hrad), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
220  { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
221 };
222 
223 /**
224  * Compute the bounding RPP for a particle
225  */
226 int
227 rt_part_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
228  struct rt_part_internal *pip;
229  vect_t tip_pt, tmp_min, tmp_max;
230  RT_CK_DB_INTERNAL(ip);
231  pip = (struct rt_part_internal *)ip->idb_ptr;
232  RT_PART_CK_MAGIC(pip);
233 
234  (*min)[X] = pip->part_V[X] - pip->part_vrad;
235  (*max)[X] = pip->part_V[X] + pip->part_vrad;
236  (*min)[Y] = pip->part_V[Y] - pip->part_vrad;
237  (*max)[Y] = pip->part_V[Y] + pip->part_vrad;
238  (*min)[Z] = pip->part_V[Z] - pip->part_vrad;
239  (*max)[Z] = pip->part_V[Z] + pip->part_vrad;
240 
241  /* If we've got a sphere, we're done */
242  if (pip->part_type == RT_PARTICLE_TYPE_SPHERE) {
243  return 0; /* OK */
244  }
245 
246  VADD2(tip_pt, pip->part_V, pip->part_H);
247  tmp_min[X] = tip_pt[X] - pip->part_hrad;
248  tmp_max[X] = tip_pt[X] + pip->part_hrad;
249  tmp_min[Y] = tip_pt[Y] - pip->part_hrad;
250  tmp_max[Y] = tip_pt[Y] + pip->part_hrad;
251  tmp_min[Z] = tip_pt[Z] - pip->part_hrad;
252  tmp_max[Z] = tip_pt[Z] + pip->part_hrad;
253  VMINMAX((*min), (*max), tmp_min);
254  VMINMAX((*min), (*max), tmp_max);
255 
256  return 0;
257 }
258 
259 
260 /**
261  * Given a pointer to a GED database record, and a transformation matrix,
262  * determine if this is a valid particle, and if so, precompute various
263  * terms of the formula.
264  *
265  * Returns -
266  * 0 particle is OK
267  * !0 Error in description
268  *
269  * Implicit return -
270  * A struct part_specific is created, and its address is stored in
271  * stp->st_specific for use by part_shot().
272  */
273 int
274 rt_part_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
275 {
276  register struct part_specific *part;
277  struct rt_part_internal *pip;
278  vect_t Hunit;
279  vect_t a, b;
280  mat_t R, Rinv;
281  mat_t S;
282  fastf_t inv_hlen;
283  fastf_t hlen;
284  fastf_t hlen_sq;
285  fastf_t r_diff;
286  fastf_t r_diff_sq;
287  fastf_t sin_theta;
288  fastf_t cos_theta;
289 
290  if (rtip) RT_CK_RTI(rtip);
291  RT_CK_DB_INTERNAL(ip);
292  pip = (struct rt_part_internal *)ip->idb_ptr;
293  RT_PART_CK_MAGIC(pip);
294 
295  BU_GET(part, struct part_specific);
296  stp->st_specific = (void *)part;
297  part->part_int = *pip; /* struct copy */
298  pip = &part->part_int;
299 
300  if (rt_part_bbox(ip, &(stp->st_min), &(stp->st_max), &rtip->rti_tol)) return 1;
301 
302  if (pip->part_type == RT_PARTICLE_TYPE_SPHERE) {
303  /* Compute bounding sphere*/
304  VMOVE(stp->st_center, pip->part_V);
305  stp->st_aradius = stp->st_bradius = pip->part_vrad;
306  return 0; /* OK */
307  }
308 
309  /* Compute some essential terms */
310  hlen_sq = MAGSQ(pip->part_H);
311  if (hlen_sq <= SMALL) {
312  bu_log("part(%s): 0-length H vector\n", stp->st_dp->d_namep);
313  return 1; /* BAD */
314  }
315  hlen = sqrt(hlen_sq);
316  inv_hlen = 1/hlen;
317  VSCALE(Hunit, pip->part_H, inv_hlen);
318  bn_vec_ortho(a, Hunit);
319  VCROSS(b, Hunit, a);
320 
321  /*
322  * Compute parameters for the "equivalent cone"
323  */
324 
325  /* Calculate changes in terms of the "unit particle" */
326  if (pip->part_vrad >= pip->part_hrad) {
327  /* V end is larger, H end is smaller */
328  r_diff = (pip->part_vrad - pip->part_hrad) * inv_hlen;
329  r_diff_sq = r_diff * r_diff;
330  if (r_diff_sq > 1) {
331  /* No "equivalent cone" is possible, theta=90deg */
332  sin_theta = 1;
333  cos_theta = 0;
334  } else {
335  sin_theta = sqrt(1 - r_diff_sq);
336  cos_theta = fabs(r_diff);
337  }
338 
339  part->part_v_erad = pip->part_vrad / sin_theta;
340  part->part_h_erad = pip->part_hrad / sin_theta;
341 
342  /* Move both plates towards H hemisphere */
343  part->part_v_hdist = cos_theta * pip->part_vrad * inv_hlen;
344  part->part_h_hdist = 1 + cos_theta * pip->part_hrad * inv_hlen;
345  } else {
346  /* H end is larger, V end is smaller */
347  r_diff = (pip->part_hrad - pip->part_vrad) * inv_hlen;
348  r_diff_sq = r_diff * r_diff;
349  if (r_diff_sq > 1) {
350  /* No "equivalent cone" is possible, theta=90deg */
351  sin_theta = 1;
352  cos_theta = 0;
353  } else {
354  sin_theta = sqrt(1 - r_diff_sq);
355  cos_theta = fabs(r_diff);
356  }
357 
358  part->part_v_erad = pip->part_vrad / sin_theta;
359  part->part_h_erad = pip->part_hrad / sin_theta;
360 
361  /* Move both plates towards V hemisphere */
362  part->part_v_hdist = -cos_theta * pip->part_vrad * inv_hlen;
363  part->part_h_hdist = 1 - cos_theta * pip->part_hrad * inv_hlen;
364  }
365  /* Thanks to matrix S, vrad_prime is always 1 */
366 /*#define VRAD_PRIME 1 */
367 /*#define HRAD_PRIME (part->part_int.part_hrad / part->part_int.part_vrad) */
368  part->part_vrad_prime = 1;
369  part->part_hrad_prime = part->part_h_erad / part->part_v_erad;
370 
371  /* Compute R and Rinv */
372  MAT_IDN(R);
373  VMOVE(&R[0], a); /* has unit length */
374  VMOVE(&R[4], b); /* has unit length */
375  VMOVE(&R[8], Hunit);
376  bn_mat_trn(Rinv, R);
377 
378  /* Compute scale matrix S, where |A| = |B| = equiv_Vradius */
379  MAT_IDN(S);
380  S[ 0] = 1.0 / part->part_v_erad;
381  S[ 5] = S[0];
382  S[10] = inv_hlen;
383 
384  bn_mat_mul(part->part_SoR, S, R);
385  bn_mat_mul(part->part_invRoS, Rinv, S);
386 
387  /* RPP and bounding sphere */
388  VJOIN1(stp->st_center, pip->part_V, 0.5, pip->part_H);
389 
390  stp->st_aradius = stp->st_bradius = pip->part_vrad;
391 
392  /* Determine bounding sphere from the RPP */
393  {
394  register fastf_t f;
395  vect_t work;
396  VSUB2SCALE(work, stp->st_max, stp->st_min, 0.5); /* radius */
397  f = work[X];
398  if (work[Y] > f) f = work[Y];
399  if (work[Z] > f) f = work[Z];
400  stp->st_aradius = f;
401  stp->st_bradius = MAGNITUDE(work);
402  }
403  return 0; /* OK */
404 }
405 
406 
407 void
408 rt_part_print(register const struct soltab *stp)
409 {
410  register const struct part_specific *part =
411  (struct part_specific *)stp->st_specific;
412 
413  VPRINT("part_V", part->part_int.part_V);
414  VPRINT("part_H", part->part_int.part_H);
415  bu_log("part_vrad=%g\n", part->part_int.part_vrad);
416  bu_log("part_hrad=%g\n", part->part_int.part_hrad);
417 
418  switch (part->part_int.part_type) {
419  case RT_PARTICLE_TYPE_SPHERE:
420  bu_log("part_type = SPHERE\n");
421  break;
422  case RT_PARTICLE_TYPE_CYLINDER:
423  bu_log("part_type = CYLINDER\n");
424  bn_mat_print("part_SoR", part->part_SoR);
425  bn_mat_print("part_invRoS", part->part_invRoS);
426  break;
427  case RT_PARTICLE_TYPE_CONE:
428  bu_log("part_type = CONE\n");
429  bn_mat_print("part_SoR", part->part_SoR);
430  bn_mat_print("part_invRoS", part->part_invRoS);
431  break;
432  default:
433  bu_log("part_type = %d ???\n", part->part_int.part_type);
434  break;
435  }
436 }
437 
438 
439 /**
440  * Intersect a ray with a part.
441  * If an intersection occurs, a struct seg will be acquired
442  * and filled in.
443  *
444  * Returns -
445  * 0 MISS
446  * >0 HIT
447  */
448 int
449 rt_part_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
450 {
451  register struct part_specific *part =
452  (struct part_specific *)stp->st_specific;
453  struct seg *segp;
454  vect_t dprime; /* D' */
455  point_t pprime; /* P' */
456  point_t xlated; /* translated ray start point */
457  fastf_t t1, t2; /* distance constants of solution */
458  fastf_t f;
459  struct hit hits[4]; /* 4 potential hit points */
460  register struct hit *hitp = &hits[0];
461  int check_v, check_h;
462 
463  if (part->part_int.part_type == RT_PARTICLE_TYPE_SPHERE) {
464  vect_t ov; /* ray origin to center (V - P) */
465  fastf_t vrad_sq;
466  fastf_t magsq_ov; /* length squared of ov */
467  fastf_t b; /* second term of quadratic eqn */
468  fastf_t root; /* root of radical */
469 
470  VSUB2(ov, part->part_int.part_V, rp->r_pt);
471  b = VDOT(rp->r_dir, ov);
472  magsq_ov = MAGSQ(ov);
473 
474  if (magsq_ov >= (vrad_sq = part->part_int.part_vrad *
475  part->part_int.part_vrad)) {
476  /* ray origin is outside of sphere */
477  if (b < 0) {
478  /* ray direction is away from sphere */
479  return 0; /* No hit */
480  }
481  root = b*b - magsq_ov + vrad_sq;
482  if (root <= 0) {
483  /* no real roots */
484  return 0; /* No hit */
485  }
486  } else {
487  root = b*b - magsq_ov + vrad_sq;
488  }
489  root = sqrt(root);
490 
491  RT_GET_SEG(segp, ap->a_resource);
492  segp->seg_stp = stp;
493 
494  /* we know root is positive, so we know the smaller t */
495  segp->seg_in.hit_magic = RT_HIT_MAGIC;
496  segp->seg_in.hit_dist = b - root;
497  segp->seg_in.hit_surfno = RT_PARTICLE_SURF_VSPHERE;
498  segp->seg_out.hit_magic = RT_HIT_MAGIC;
499  segp->seg_out.hit_dist = b + root;
500  segp->seg_out.hit_surfno = RT_PARTICLE_SURF_VSPHERE;
501  BU_LIST_INSERT(&(seghead->l), &(segp->l));
502  return 2; /* HIT */
503  }
504 
505  /* Transform ray to coordinate system of unit cone at origin */
506  MAT4X3VEC(dprime, part->part_SoR, rp->r_dir);
507  VSUB2(xlated, rp->r_pt, part->part_int.part_V);
508  MAT4X3VEC(pprime, part->part_SoR, xlated);
509 
510  if (ZERO(dprime[X]) && ZERO(dprime[Y])) {
511  check_v = check_h = 1;
512  goto check_hemispheres;
513  }
514  check_v = check_h = 0;
515 
516  /* Find roots of the equation, using formula for quadratic */
517  /* Note that vrad' = 1 and hrad' = hrad/vrad */
518  if (part->part_int.part_type == RT_PARTICLE_TYPE_CYLINDER) {
519  /* Cylinder case, hrad == vrad, m = 0 */
520  fastf_t a, b, c;
521  fastf_t root; /* root of radical */
522 
523  a = dprime[X]*dprime[X] + dprime[Y]*dprime[Y];
524  b = dprime[X]*pprime[X] + dprime[Y]*pprime[Y];
525  c = pprime[X]*pprime[X] + pprime[Y]*pprime[Y] - 1;
526  if ((root = b*b - a * c) <= 0)
527  goto check_hemispheres;
528  root = sqrt(root);
529  t1 = (root-b) / a;
530  t2 = -(root+b) / a;
531  } else {
532  /* Cone case */
533  fastf_t a, b, c;
534  fastf_t root; /* root of radical */
535  fastf_t m, msq;
536 
537  m = part->part_hrad_prime - part->part_vrad_prime;
538 
539  /* This quadratic has had a factor of 2 divided out of "b"
540  * throughout. More efficient, but the same answers.
541  */
542  a = dprime[X]*dprime[X] + dprime[Y]*dprime[Y] -
543  (msq = m*m) * dprime[Z]*dprime[Z];
544  b = dprime[X]*pprime[X] + dprime[Y]*pprime[Y] -
545  msq * dprime[Z]*pprime[Z] -
546  m * dprime[Z]; /* * part->part_vrad_prime */
547  c = pprime[X]*pprime[X] + pprime[Y]*pprime[Y] -
548  msq * pprime[Z]*pprime[Z] -
549  2 * m * pprime[Z] - 1;
550  /* was: ... -2m * vrad' * Pz' - vrad'**2 */
551 
552  if ((root = b*b - a * c) <= 0)
553  goto check_hemispheres;
554  root = sqrt(root);
555 
556  t1 = (root-b) / a;
557  t2 = -(root+b) / a;
558  }
559 
560  /*
561  * t1 and t2 are potential solutions to intersection with side.
562  * Find hit' point, see if Z values fall in range.
563  */
564  if ((f = pprime[Z] + t1 * dprime[Z]) >= part->part_v_hdist) {
565  check_h = 1; /* may also hit off end */
566  if (f <= part->part_h_hdist) {
567  hitp->hit_magic = RT_HIT_MAGIC;
568  /** VJOIN1(hitp->hit_vpriv, pprime, t1, dprime); **/
569  hitp->hit_vpriv[X] = pprime[X] + t1 * dprime[X];
570  hitp->hit_vpriv[Y] = pprime[Y] + t1 * dprime[Y];
571  hitp->hit_vpriv[Z] = f;
572  hitp->hit_dist = t1;
573  hitp->hit_surfno = RT_PARTICLE_SURF_BODY;
574  hitp++;
575  }
576  } else {
577  check_v = 1;
578  }
579 
580  if ((f = pprime[Z] + t2 * dprime[Z]) >= part->part_v_hdist) {
581  check_h = 1; /* may also hit off end */
582  if (f <= part->part_h_hdist) {
583  hitp->hit_magic = RT_HIT_MAGIC;
584  /** VJOIN1(hitp->hit_vpriv, pprime, t2, dprime); **/
585  hitp->hit_vpriv[X] = pprime[X] + t2 * dprime[X];
586  hitp->hit_vpriv[Y] = pprime[Y] + t2 * dprime[Y];
587  hitp->hit_vpriv[Z] = f;
588  hitp->hit_dist = t2;
589  hitp->hit_surfno = RT_PARTICLE_SURF_BODY;
590  hitp++;
591  }
592  } else {
593  check_v = 1;
594  }
595 
596  /*
597  * Check for hitting the end hemispheres.
598  */
599  check_hemispheres:
600  if (check_v) {
601  vect_t ov; /* ray origin to center (V - P) */
602  fastf_t rad_sq;
603  fastf_t magsq_ov; /* length squared of ov */
604  fastf_t b;
605  fastf_t root; /* root of radical */
606 
607  /*
608  * First, consider a hit on V hemisphere.
609  */
610  VSUB2(ov, part->part_int.part_V, rp->r_pt);
611  b = VDOT(rp->r_dir, ov);
612  magsq_ov = MAGSQ(ov);
613  if (magsq_ov >= (rad_sq = part->part_int.part_vrad *
614  part->part_int.part_vrad)) {
615  /* ray origin is outside of sphere */
616  if (b < 0) {
617  /* ray direction is away from sphere */
618  goto do_check_h;
619  }
620  root = b*b - magsq_ov + rad_sq;
621  if (root <= 0) {
622  /* no real roots */
623  goto do_check_h;
624  }
625  } else {
626  root = b*b - magsq_ov + rad_sq;
627  }
628  root = sqrt(root);
629  t1 = b - root;
630  /* see if hit'[Z] is below V end of cylinder */
631  if (pprime[Z] + t1 * dprime[Z] <= part->part_v_hdist) {
632  hitp->hit_magic = RT_HIT_MAGIC;
633  hitp->hit_dist = t1;
634  hitp->hit_surfno = RT_PARTICLE_SURF_VSPHERE;
635  hitp++;
636  }
637  t2 = b + root;
638  if (pprime[Z] + t2 * dprime[Z] <= part->part_v_hdist) {
639  hitp->hit_magic = RT_HIT_MAGIC;
640  hitp->hit_dist = t2;
641  hitp->hit_surfno = RT_PARTICLE_SURF_VSPHERE;
642  hitp++;
643  }
644  }
645 
646  do_check_h:
647  if (check_h) {
648  vect_t ov; /* ray origin to center (V - P) */
649  fastf_t rad_sq;
650  fastf_t magsq_ov; /* length squared of ov */
651  fastf_t b; /* second term of quadratic eqn */
652  fastf_t root; /* root of radical */
653 
654  /*
655  * Next, consider a hit on H hemisphere
656  */
657  VADD2(ov, part->part_int.part_V, part->part_int.part_H);
658  VSUB2(ov, ov, rp->r_pt);
659  b = VDOT(rp->r_dir, ov);
660  magsq_ov = MAGSQ(ov);
661  if (magsq_ov >= (rad_sq = part->part_int.part_hrad *
662  part->part_int.part_hrad)) {
663  /* ray origin is outside of sphere */
664  if (b < 0) {
665  /* ray direction is away from sphere */
666  goto out;
667  }
668  root = b*b - magsq_ov + rad_sq;
669  if (root <= 0) {
670  /* no real roots */
671  goto out;
672  }
673  } else {
674  root = b*b - magsq_ov + rad_sq;
675  }
676  root = sqrt(root);
677  t1 = b - root;
678  /* see if hit'[Z] is above H end of cylinder */
679  if (pprime[Z] + t1 * dprime[Z] >= part->part_h_hdist) {
680  hitp->hit_magic = RT_HIT_MAGIC;
681  hitp->hit_dist = t1;
682  hitp->hit_surfno = RT_PARTICLE_SURF_HSPHERE;
683  hitp++;
684  }
685  t2 = b + root;
686  if (pprime[Z] + t2 * dprime[Z] >= part->part_h_hdist) {
687  hitp->hit_magic = RT_HIT_MAGIC;
688  hitp->hit_dist = t2;
689  hitp->hit_surfno = RT_PARTICLE_SURF_HSPHERE;
690  hitp++;
691  }
692  }
693  out:
694  if (hitp == &hits[0])
695  return 0; /* MISS */
696  if (hitp == &hits[1]) {
697  /* Only one hit, make it a 0-thickness segment */
698  hits[1] = hits[0]; /* struct copy */
699  hitp++;
700  } else if (hitp > &hits[2]) {
701  /*
702  * More than two intersections found.
703  * This can happen when a ray grazes down along a tangent
704  * line; the intersection interval from the hemisphere
705  * may not quite join up with the interval from the cone.
706  * Since particles are convex, all we need to do is to
707  * return the maximum extent of the ray.
708  * Do this by sorting the intersections,
709  * and using the minimum and maximum values.
710  */
711  rt_hitsort(hits, hitp - &hits[0]);
712 
713  /* [0] is minimum, make [1] be maximum (hitp is +1 off end) */
714  hits[1] = hitp[-1]; /* struct copy */
715  }
716 
717  if (hits[0].hit_dist < hits[1].hit_dist) {
718  /* entry is [0], exit is [1] */
719  RT_GET_SEG(segp, ap->a_resource);
720  segp->seg_stp = stp;
721  segp->seg_in = hits[0]; /* struct copy */
722  segp->seg_out = hits[1]; /* struct copy */
723  BU_LIST_INSERT(&(seghead->l), &(segp->l));
724  } else {
725  /* entry is [1], exit is [0] */
726  RT_GET_SEG(segp, ap->a_resource);
727  segp->seg_stp = stp;
728  segp->seg_in = hits[1]; /* struct copy */
729  segp->seg_out = hits[0]; /* struct copy */
730  BU_LIST_INSERT(&(seghead->l), &(segp->l));
731  }
732  return 2; /* HIT */
733 }
734 
735 
736 /**
737  * Given ONE ray distance, return the normal and entry/exit point.
738  */
739 void
740 rt_part_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
741 {
742  register struct part_specific *part =
743  (struct part_specific *)stp->st_specific;
744 
745  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
746  switch (hitp->hit_surfno) {
748  VSUB2(hitp->hit_normal, hitp->hit_point, part->part_int.part_V);
749  VUNITIZE(hitp->hit_normal);
750  break;
752  VSUB3(hitp->hit_normal, hitp->hit_point,
753  part->part_int.part_V, part->part_int.part_H);
754  VUNITIZE(hitp->hit_normal);
755  break;
757  /* compute it */
758  if (part->part_int.part_type == RT_PARTICLE_TYPE_CYLINDER) {
759  /* The X' and Y' components of hit' are N' */
760  hitp->hit_vpriv[Z] = 0;
761  MAT4X3VEC(hitp->hit_normal, part->part_invRoS,
762  hitp->hit_vpriv);
763  VUNITIZE(hitp->hit_normal);
764  } else {
765  /* The cone case */
766  fastf_t s, m;
767  vect_t unorm;
768  /* vpriv[Z] ranges from 0 to 1 (roughly) */
769  /* Rescale X' and Y' into unit circle */
770  m = part->part_hrad_prime - part->part_vrad_prime;
771  s = 1/(part->part_vrad_prime + m * hitp->hit_vpriv[Z]);
772  unorm[X] = hitp->hit_vpriv[X] * s;
773  unorm[Y] = hitp->hit_vpriv[Y] * s;
774  /* Z' is constant, from slope of cylinder wall*/
775  unorm[Z] = -m / sqrt(m*m+1);
776  MAT4X3VEC(hitp->hit_normal, part->part_invRoS, unorm);
777  VUNITIZE(hitp->hit_normal);
778  }
779  break;
780  }
781 }
782 
783 
784 /**
785  * Return the curvature of the particle.
786  * There are two cases: hitting a hemisphere, and hitting the cylinder.
787  */
788 void
789 rt_part_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
790 {
791  register struct part_specific *part =
792  (struct part_specific *)stp->st_specific;
793  point_t hit_local; /* hit_point, with V as origin */
794  point_t hit_unit; /* hit_point in unit coords, +Z along H */
795 
796  switch (hitp->hit_surfno) {
798  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
799  cvp->crv_c1 = cvp->crv_c2 = -part->part_int.part_vrad;
800  break;
802  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
803  cvp->crv_c1 = cvp->crv_c2 = -part->part_int.part_hrad;
804  break;
806  /* Curvature in only one direction, around H */
807  VCROSS(cvp->crv_pdir, hitp->hit_normal, part->part_int.part_H);
808  VUNITIZE(cvp->crv_pdir);
809  /* Interpolate radius between vrad and hrad */
810  VSUB2(hit_local, hitp->hit_point, part->part_int.part_V);
811  MAT4X3VEC(hit_unit, part->part_SoR, hit_local);
812  /* hit_unit[Z] ranges from 0 at V to 1 at H, interpolate */
813  cvp->crv_c1 = -(
814  part->part_v_erad * hit_unit[Z] +
815  part->part_h_erad * (1 - hit_unit[Z]));
816  cvp->crv_c2 = 0;
817  break;
818  }
819 }
820 
821 
822 /**
823  * For a hit on the surface of a particle, return the (u, v) coordinates
824  * of the hit point, 0 <= u, v <= 1.
825  * u = azimuth
826  * v = elevation along H
827  *
828  * The 'u' coordinate wraps around the particle, once.
829  * The 'v' coordinate covers the 'height' of the particle,
830  * from V-r1 to (V+H)+r2.
831  *
832  * hit_point has already been computed.
833  */
834 void
835 rt_part_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
836 {
837  register const struct part_specific *part =
838  (struct part_specific *)stp->st_specific;
839  point_t hit_local; /* hit_point, with V as origin */
840  point_t hit_unit; /* hit_poit in unit coords, +Z along H */
841  fastf_t hsize;
842  fastf_t hmag_inv;
843  fastf_t vrad_unit;
844  fastf_t r;
845  fastf_t minrad;
846 
847  RT_PART_CK_MAGIC(&part->part_int.part_magic);
848 
849  hmag_inv = 1.0/MAGNITUDE(part->part_int.part_H);
850  hsize = 1 + (vrad_unit = part->part_v_erad*hmag_inv) +
851  part->part_h_erad*hmag_inv;
852 
853  /* Transform hit point into unit particle coords */
854  VSUB2(hit_local, hitp->hit_point, part->part_int.part_V);
855  MAT4X3VEC(hit_unit, part->part_SoR, hit_local);
856  /* normalize 0..1 */
857  uvp->uv_v = (hit_unit[Z] + vrad_unit) / hsize;
858 
859  /* U is azimuth, atan() range: -pi to +pi */
860  uvp->uv_u = bn_atan2(hit_unit[Y], hit_unit[X]) * M_1_2PI;
861  if (uvp->uv_u < 0)
862  uvp->uv_u += 1.0;
863 
864  /* approximation: beam_r / (solid_circumference = 2 * pi * radius) */
865  minrad = part->part_v_erad;
866  V_MIN(minrad, part->part_h_erad);
867  r = ap->a_rbeam + ap->a_diverge * hitp->hit_dist;
868  uvp->uv_du = uvp->uv_dv =
869  M_1_2PI * r / minrad;
870 }
871 
872 
873 void
874 rt_part_free(register struct soltab *stp)
875 {
876  register struct part_specific *part =
877  (struct part_specific *)stp->st_specific;
878 
879  BU_PUT(part, struct part_specific);
880  stp->st_specific = ((void *)0);
881 }
882 
883 
884 /**
885  * Produce a crude approximation to a hemisphere,
886  * 8 points around the rim [0]..[7],
887  * 4 points around a midway latitude [8]..[11], and
888  * 1 point at the pole [12].
889  *
890  * For the dome, connect up:
891  * 0 8 12 10 4
892  * 2 9 12 11 6
893  */
894 HIDDEN void
895 rt_part_hemisphere(register point_t (*ov), register fastf_t *v, fastf_t *a, fastf_t *b, fastf_t *h)
896 {
897  /* M_SQRT1_2 is cos45 */
898 
899  /* This is the top of the dome */
900  VADD2(ov[12], v, h);
901 
902  VADD2(ov[0], v, a);
903  VJOIN2(ov[1], v, M_SQRT1_2, a, M_SQRT1_2, b);
904  VADD2(ov[2], v, b);
905  VJOIN2(ov[3], v, -M_SQRT1_2, a, M_SQRT1_2, b);
906  VSUB2(ov[4], v, a);
907  VJOIN2(ov[5], v, -M_SQRT1_2, a, -M_SQRT1_2, b);
908  VSUB2(ov[6], v, b);
909  VJOIN2(ov[7], v, M_SQRT1_2, a, -M_SQRT1_2, b);
910 
911  VJOIN2(ov[8], v, M_SQRT1_2, a, M_SQRT1_2, h);
912  VJOIN2(ov[10], v, -M_SQRT1_2, a, M_SQRT1_2, h);
913 
914  VJOIN2(ov[9], v, M_SQRT1_2, b, M_SQRT1_2, h);
915  VJOIN2(ov[11], v, -M_SQRT1_2, b, M_SQRT1_2, h);
916  /* Obviously, this could be optimized quite a lot more */
917 }
918 
919 
920 int
921 rt_part_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *UNUSED(ttol), const struct bn_tol *UNUSED(tol), const struct rt_view_info *UNUSED(info))
922 {
923  struct rt_part_internal *pip;
924  point_t tail;
925  point_t sphere_rim[16];
926  point_t vhemi[13];
927  point_t hhemi[13];
928  vect_t a, b, c; /* defines coord sys of part */
929  vect_t as, bs, hs; /* scaled by radius */
930  vect_t Hunit;
931  register int i;
932 
933  BU_CK_LIST_HEAD(vhead);
934  RT_CK_DB_INTERNAL(ip);
935  pip = (struct rt_part_internal *)ip->idb_ptr;
936  RT_PART_CK_MAGIC(pip);
937 
938  if (pip->part_type == RT_PARTICLE_TYPE_SPHERE) {
939  /* For the sphere, 3 rings of 16 points */
940  VSET(a, pip->part_vrad, 0, 0);
941  VSET(b, 0, pip->part_vrad, 0);
942  VSET(c, 0, 0, pip->part_vrad);
943 
944  rt_ell_16pts(&sphere_rim[0][X], pip->part_V, a, b);
945  RT_ADD_VLIST(vhead, sphere_rim[15], BN_VLIST_LINE_MOVE);
946  for (i=0; i<16; i++) {
947  RT_ADD_VLIST(vhead, sphere_rim[i], BN_VLIST_LINE_DRAW);
948  }
949 
950  rt_ell_16pts(&sphere_rim[0][X], pip->part_V, b, c);
951  RT_ADD_VLIST(vhead, sphere_rim[15], BN_VLIST_LINE_MOVE);
952  for (i=0; i<16; i++) {
953  RT_ADD_VLIST(vhead, sphere_rim[i], BN_VLIST_LINE_DRAW);
954  }
955 
956  rt_ell_16pts(&sphere_rim[0][X], pip->part_V, a, c);
957  RT_ADD_VLIST(vhead, sphere_rim[15], BN_VLIST_LINE_MOVE);
958  for (i=0; i<16; i++) {
959  RT_ADD_VLIST(vhead, sphere_rim[i], BN_VLIST_LINE_DRAW);
960  }
961  return 0; /* OK */
962  }
963 
964  VMOVE(Hunit, pip->part_H);
965  VUNITIZE(Hunit);
966  bn_vec_perp(a, Hunit);
967  VUNITIZE(a);
968  VCROSS(b, Hunit, a);
969  VUNITIZE(b);
970 
971  VSCALE(as, a, pip->part_vrad);
972  VSCALE(bs, b, pip->part_vrad);
973  VSCALE(hs, Hunit, -pip->part_vrad);
974  rt_part_hemisphere(vhemi, pip->part_V, as, bs, hs);
975 
976  VADD2(tail, pip->part_V, pip->part_H);
977  VSCALE(as, a, pip->part_hrad);
978  VSCALE(bs, b, pip->part_hrad);
979  VSCALE(hs, Hunit, pip->part_hrad);
980  rt_part_hemisphere(hhemi, tail, as, bs, hs);
981 
982  /* Draw V end hemisphere */
983  RT_ADD_VLIST(vhead, vhemi[0], BN_VLIST_LINE_MOVE);
984  for (i=7; i >= 0; i--) {
985  RT_ADD_VLIST(vhead, vhemi[i], BN_VLIST_LINE_DRAW);
986  }
987  RT_ADD_VLIST(vhead, vhemi[8], BN_VLIST_LINE_DRAW);
988  RT_ADD_VLIST(vhead, vhemi[12], BN_VLIST_LINE_DRAW);
989  RT_ADD_VLIST(vhead, vhemi[10], BN_VLIST_LINE_DRAW);
990  RT_ADD_VLIST(vhead, vhemi[4], BN_VLIST_LINE_DRAW);
991  RT_ADD_VLIST(vhead, vhemi[2], BN_VLIST_LINE_MOVE);
992  RT_ADD_VLIST(vhead, vhemi[9], BN_VLIST_LINE_DRAW);
993  RT_ADD_VLIST(vhead, vhemi[12], BN_VLIST_LINE_DRAW);
994  RT_ADD_VLIST(vhead, vhemi[11], BN_VLIST_LINE_DRAW);
995  RT_ADD_VLIST(vhead, vhemi[6], BN_VLIST_LINE_DRAW);
996 
997  /* Draw H end hemisphere */
998  RT_ADD_VLIST(vhead, hhemi[0], BN_VLIST_LINE_MOVE);
999  for (i=7; i >= 0; i--) {
1000  RT_ADD_VLIST(vhead, hhemi[i], BN_VLIST_LINE_DRAW);
1001  }
1002  RT_ADD_VLIST(vhead, hhemi[8], BN_VLIST_LINE_DRAW);
1003  RT_ADD_VLIST(vhead, hhemi[12], BN_VLIST_LINE_DRAW);
1004  RT_ADD_VLIST(vhead, hhemi[10], BN_VLIST_LINE_DRAW);
1005  RT_ADD_VLIST(vhead, hhemi[4], BN_VLIST_LINE_DRAW);
1006  RT_ADD_VLIST(vhead, hhemi[2], BN_VLIST_LINE_MOVE);
1007  RT_ADD_VLIST(vhead, hhemi[9], BN_VLIST_LINE_DRAW);
1008  RT_ADD_VLIST(vhead, hhemi[12], BN_VLIST_LINE_DRAW);
1009  RT_ADD_VLIST(vhead, hhemi[11], BN_VLIST_LINE_DRAW);
1010  RT_ADD_VLIST(vhead, hhemi[6], BN_VLIST_LINE_DRAW);
1011 
1012  /* Draw 4 connecting lines */
1013  RT_ADD_VLIST(vhead, vhemi[0], BN_VLIST_LINE_MOVE);
1014  RT_ADD_VLIST(vhead, hhemi[0], BN_VLIST_LINE_DRAW);
1015  RT_ADD_VLIST(vhead, vhemi[2], BN_VLIST_LINE_MOVE);
1016  RT_ADD_VLIST(vhead, hhemi[2], BN_VLIST_LINE_DRAW);
1017  RT_ADD_VLIST(vhead, vhemi[4], BN_VLIST_LINE_MOVE);
1018  RT_ADD_VLIST(vhead, hhemi[4], BN_VLIST_LINE_DRAW);
1019  RT_ADD_VLIST(vhead, vhemi[6], BN_VLIST_LINE_MOVE);
1020  RT_ADD_VLIST(vhead, hhemi[6], BN_VLIST_LINE_DRAW);
1021 
1022  return 0;
1023 }
1024 
1025 
1026 struct part_state {
1027  struct shell *s;
1033 };
1034 
1035 
1038  int nverts;
1039  struct vertex **vp;
1040  vect_t *norms;
1041  int nfaces;
1042  struct faceuse **fu;
1043 };
1044 
1045 
1046 /**
1047  * Based upon the tesselator for the ellipsoid.
1048  *
1049  * Break the particle into three parts:
1050  * Upper hemisphere 0..nsegs H North
1051  * middle cylinder nsegs..nsegs+1
1052  * lower hemisphere nsegs+1..nstrips-1 V South
1053  */
1054 int
1055 rt_part_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
1056 {
1057  struct rt_part_internal *pip;
1058  mat_t R;
1059  mat_t S;
1060  mat_t invR;
1061  mat_t invS;
1062  vect_t zz;
1063  vect_t hcenter;
1064  struct part_state state;
1065  register int i;
1066  fastf_t radius;
1067  int nsegs;
1068  int nstrips;
1069  struct part_vert_strip *strips;
1070  int j;
1071  struct vertex **vertp[5];
1072  int faceno;
1073  int stripno;
1074  int boff; /* base offset */
1075  int toff; /* top offset */
1076  int blim; /* base subscript limit */
1077  int tlim; /* top subscript limit */
1078  fastf_t dtol; /* Absolutized relative tolerance */
1079 
1080  RT_CK_DB_INTERNAL(ip);
1081  pip = (struct rt_part_internal *)ip->idb_ptr;
1082  RT_PART_CK_MAGIC(pip);
1083 
1084  if (pip->part_type == RT_PARTICLE_TYPE_SPHERE)
1085  return -1;
1086  /* For now, concentrate on the most important kind. */
1087 
1088  VADD2(hcenter, pip->part_V, pip->part_H);
1089 
1090  /* Compute R and Rinv matrices */
1091  /* R is the same for both upper and lower hemispheres */
1092  /* R is rotation from model coords to unit sphere */
1093  /* For rotation of the particle, +H (model) becomes +Z (unit sph) */
1094  VSET(zz, 0, 0, 1);
1095  bn_mat_fromto(R, pip->part_H, zz, tol);
1096  bn_mat_trn(invR, R); /* inv of rot mat is trn */
1097 
1098  /*** Upper (H) ***/
1099 
1100  /* Compute S and invS matrices */
1101  /* invS is just 1/diagonal elements */
1102  MAT_IDN(S);
1103  S[ 0] = S[ 5] = S[10] = 1/pip->part_hrad;
1104  bn_mat_inv(invS, S);
1105 
1106  /* invRinvS, for converting points from unit sphere to model */
1107  bn_mat_mul(state.upper_invRinvS, invR, invS);
1108 
1109  /* invRoS, for converting normals from unit sphere to model */
1110  bn_mat_mul(state.upper_invRoS, invR, S);
1111 
1112  /*** Lower (V) ***/
1113 
1114  /* Compute S and invS matrices */
1115  /* invS is just 1/diagonal elements */
1116  MAT_IDN(S);
1117  S[ 0] = S[ 5] = S[10] = 1/pip->part_vrad;
1118  bn_mat_inv(invS, S);
1119 
1120  /* invRinvS, for converting points from unit sphere to model */
1121  bn_mat_mul(state.lower_invRinvS, invR, invS);
1122 
1123  /* invRoS, for converting normals from unit sphere to model */
1124  bn_mat_mul(state.lower_invRoS, invR, S);
1125 
1126  /* Find the larger of two hemispheres */
1127  radius = pip->part_vrad;
1128  if (pip->part_hrad > radius)
1129  radius = pip->part_hrad;
1130 
1131  dtol = primitive_get_absolute_tolerance(ttol, radius);
1132 
1133  if (dtol > radius) {
1134  dtol = radius;
1135  }
1136 
1137  /*
1138  * Convert distance tolerance into a maximum permissible
1139  * angle tolerance. 'radius' is largest radius.
1140  */
1141  state.theta_tol = 2.0 * acos(1.0 - dtol / radius);
1142 
1143  /* To ensure normal tolerance, remain below this angle */
1144  if (ttol->norm > 0.0 && ttol->norm < state.theta_tol) {
1145  state.theta_tol = ttol->norm;
1146  }
1147 
1148  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
1149  state.s = BU_LIST_FIRST(shell, &(*r)->s_hd);
1150 
1151  /* Find the number of segments to divide 90 degrees worth into */
1152  nsegs = M_PI_2 / state.theta_tol + 0.999;
1153  if (nsegs < 2) nsegs = 2;
1154 
1155  /* Find total number of strips of vertices that will be needed.
1156  * nsegs for each hemisphere, plus one equator each.
1157  * The two equators will be stitched together to make the cylinder.
1158  * Note that faces are listed in the stripe ABOVE, i.e., toward
1159  * the poles. Thus, strips[0] will have 4 faces.
1160  */
1161  nstrips = 2 * nsegs + 2;
1162  strips = (struct part_vert_strip *)bu_calloc(nstrips,
1163  sizeof(struct part_vert_strip), "strips[]");
1164 
1165  /* North pole (Upper hemisphere, H end) */
1166  strips[0].nverts = 1;
1167  strips[0].nverts_per_strip = 0;
1168  strips[0].nfaces = 4;
1169  /* South pole (Lower hemisphere, V end) */
1170  strips[nstrips-1].nverts = 1;
1171  strips[nstrips-1].nverts_per_strip = 0;
1172  strips[nstrips-1].nfaces = 4;
1173  /* upper equator (has faces) */
1174  strips[nsegs].nverts = nsegs * 4;
1175  strips[nsegs].nverts_per_strip = nsegs;
1176  strips[nsegs].nfaces = nsegs * 4;
1177  /* lower equator (no faces) */
1178  strips[nsegs+1].nverts = nsegs * 4;
1179  strips[nsegs+1].nverts_per_strip = nsegs;
1180  strips[nsegs+1].nfaces = 0;
1181 
1182  for (i=1; i<nsegs; i++) {
1183  strips[i].nverts_per_strip =
1184  strips[nstrips-1-i].nverts_per_strip = i;
1185  strips[i].nverts =
1186  strips[nstrips-1-i].nverts = i * 4;
1187  strips[i].nfaces =
1188  strips[nstrips-1-i].nfaces = (2 * i + 1)*4;
1189  }
1190  /* All strips have vertices and normals */
1191  for (i=0; i<nstrips; i++) {
1192  strips[i].vp = (struct vertex **)bu_calloc(strips[i].nverts,
1193  sizeof(struct vertex *), "strip vertex[]");
1194  strips[i].norms = (vect_t *)bu_calloc(strips[i].nverts,
1195  sizeof(vect_t), "strip normals[]");
1196  }
1197  /* All strips have faces, except for the one (marked) equator */
1198  for (i=0; i < nstrips; i++) {
1199  if (strips[i].nfaces <= 0) {
1200  strips[i].fu = (struct faceuse **)NULL;
1201  continue;
1202  }
1203  strips[i].fu = (struct faceuse **)bu_calloc(strips[i].nfaces,
1204  sizeof(struct faceuse *), "strip faceuse[]");
1205  }
1206 
1207  /* First, build the triangular mesh topology */
1208  /* Do the top. "toff" in i-1 is UP, towards +B */
1209  for (i = 1; i <= nsegs; i++) {
1210  faceno = 0;
1211  tlim = strips[i-1].nverts;
1212  blim = strips[i].nverts;
1213  for (stripno=0; stripno<4; stripno++) {
1214  toff = stripno * strips[i-1].nverts_per_strip;
1215  boff = stripno * strips[i].nverts_per_strip;
1216 
1217  /* Connect this quarter strip */
1218  for (j = 0; j < strips[i].nverts_per_strip; j++) {
1219 
1220  /* "Right-side-up" triangle */
1221  vertp[0] = &(strips[i].vp[j+boff]);
1222  vertp[1] = &(strips[i].vp[(j+1+boff)%blim]);
1223  vertp[2] = &(strips[i-1].vp[(j+toff)%tlim]);
1224  if ((strips[i-1].fu[faceno++] = nmg_cmface(state.s, vertp, 3)) == 0) {
1225  bu_log("rt_part_tess() nmg_cmface failure\n");
1226  goto fail;
1227  }
1228  if (j+1 >= strips[i].nverts_per_strip) break;
1229 
1230  /* Follow with interior "Up-side-down" triangle */
1231  vertp[0] = &(strips[i].vp[(j+1+boff)%blim]);
1232  vertp[1] = &(strips[i-1].vp[(j+1+toff)%tlim]);
1233  vertp[2] = &(strips[i-1].vp[(j+toff)%tlim]);
1234  if ((strips[i-1].fu[faceno++] = nmg_cmface(state.s, vertp, 3)) == 0) {
1235  bu_log("rt_part_tess() nmg_cmface failure\n");
1236  goto fail;
1237  }
1238  }
1239  }
1240  }
1241 
1242  /* Connect the two equators with rectangular (4-point) faces */
1243  i = nsegs+1;
1244  {
1245  faceno = 0;
1246  tlim = strips[i-1].nverts;
1247  blim = strips[i].nverts;
1248  {
1249  /* Connect this whole strip */
1250  for (j = 0; j < strips[i].nverts; j++) {
1251 
1252  /* "Right-side-up" rectangle */
1253  vertp[3] = &(strips[i].vp[(j)%blim]);
1254  vertp[2] = &(strips[i-1].vp[(j)%tlim]);
1255  vertp[1] = &(strips[i-1].vp[(j+1)%tlim]);
1256  vertp[0] = &(strips[i].vp[(j+1)%blim]);
1257  if ((strips[i-1].fu[faceno++] = nmg_cmface(state.s, vertp, 4)) == 0) {
1258  bu_log("rt_part_tess() nmg_cmface failure\n");
1259  goto fail;
1260  }
1261  }
1262  }
1263  }
1264 
1265  /* Do the bottom. Everything is upside down. "toff" in i+1 is DOWN */
1266  for (i = nsegs+1; i < nstrips; i++) {
1267  faceno = 0;
1268  tlim = strips[i+1].nverts;
1269  blim = strips[i].nverts;
1270  for (stripno=0; stripno<4; stripno++) {
1271  toff = stripno * strips[i+1].nverts_per_strip;
1272  boff = stripno * strips[i].nverts_per_strip;
1273 
1274  /* Connect this quarter strip */
1275  for (j = 0; j < strips[i].nverts_per_strip; j++) {
1276 
1277  /* "Right-side-up" triangle */
1278  vertp[0] = &(strips[i].vp[j+boff]);
1279  vertp[1] = &(strips[i+1].vp[(j+toff)%tlim]);
1280  vertp[2] = &(strips[i].vp[(j+1+boff)%blim]);
1281  if ((strips[i+1].fu[faceno++] = nmg_cmface(state.s, vertp, 3)) == 0) {
1282  bu_log("rt_part_tess() nmg_cmface failure\n");
1283  goto fail;
1284  }
1285  if (j+1 >= strips[i].nverts_per_strip) break;
1286 
1287  /* Follow with interior "Up-side-down" triangle */
1288  vertp[0] = &(strips[i].vp[(j+1+boff)%blim]);
1289  vertp[1] = &(strips[i+1].vp[(j+toff)%tlim]);
1290  vertp[2] = &(strips[i+1].vp[(j+1+toff)%tlim]);
1291  if ((strips[i+1].fu[faceno++] = nmg_cmface(state.s, vertp, 3)) == 0) {
1292  bu_log("rt_part_tess() nmg_cmface failure\n");
1293  goto fail;
1294  }
1295  }
1296  }
1297  }
1298 
1299  /* Compute the geometry of each vertex.
1300  * Start with the location in the unit sphere, and project back.
1301  * i=0 is "straight up" along +H.
1302  */
1303  for (i=0; i < nstrips; i++) {
1304  double alpha; /* decline down from B to A */
1305  double beta; /* angle around equator (azimuth) */
1306  fastf_t cos_alpha, sin_alpha;
1307  fastf_t cos_beta, sin_beta;
1308  point_t sphere_pt;
1309  point_t model_pt;
1310 
1311  /* Compensate for extra equator */
1312  if (i <= nsegs)
1313  alpha = (((double)i) / (nstrips-1-1));
1314  else
1315  alpha = (((double)i-1) / (nstrips-1-1));
1316  cos_alpha = cos(alpha*M_PI);
1317  sin_alpha = sin(alpha*M_PI);
1318  for (j=0; j < strips[i].nverts; j++) {
1319 
1320  beta = ((double)j) / strips[i].nverts;
1321  cos_beta = cos(beta*M_2PI);
1322  sin_beta = sin(beta*M_2PI);
1323  VSET(sphere_pt,
1324  cos_beta * sin_alpha,
1325  sin_beta * sin_alpha,
1326  cos_alpha);
1327  /* Convert from ideal sphere coordinates */
1328  if (i <= nsegs) {
1329  MAT4X3PNT(model_pt, state.upper_invRinvS, sphere_pt);
1330  VADD2(model_pt, model_pt, hcenter);
1331  /* Convert sphere normal to ellipsoid normal */
1332  MAT4X3VEC(strips[i].norms[j], state.upper_invRoS, sphere_pt);
1333  } else {
1334  MAT4X3PNT(model_pt, state.lower_invRinvS, sphere_pt);
1335  VADD2(model_pt, model_pt, pip->part_V);
1336  /* Convert sphere normal to ellipsoid normal */
1337  MAT4X3VEC(strips[i].norms[j], state.lower_invRoS, sphere_pt);
1338  }
1339 
1340  /* May not be unit length anymore */
1341  VUNITIZE(strips[i].norms[j]);
1342  /* Associate vertex geometry */
1343  nmg_vertex_gv(strips[i].vp[j], model_pt);
1344  }
1345  }
1346 
1347  /* Associate face geometry. lower Equator has no faces */
1348  for (i=0; i < nstrips; i++) {
1349  for (j=0; j < strips[i].nfaces; j++) {
1350  if (nmg_fu_planeeqn(strips[i].fu[j], tol) < 0)
1351  goto fail;
1352  }
1353  }
1354 
1355  /* Associate normals with vertexuses */
1356  for (i=0; i < nstrips; i++) {
1357  for (j=0; j < strips[i].nverts; j++) {
1358  struct faceuse *fu;
1359  struct vertexuse *vu;
1360  vect_t norm_opp;
1361 
1362  NMG_CK_VERTEX(strips[i].vp[j]);
1363  VREVERSE(norm_opp, strips[i].norms[j]);
1364 
1365  for (BU_LIST_FOR(vu, vertexuse, &strips[i].vp[j]->vu_hd)) {
1366  fu = nmg_find_fu_of_vu(vu);
1367  NMG_CK_FACEUSE(fu);
1368  /* get correct direction of normals depending on
1369  * faceuse orientation
1370  */
1371  if (fu->orientation == OT_SAME)
1372  nmg_vertexuse_nv(vu, strips[i].norms[j]);
1373  else if (fu->orientation == OT_OPPOSITE)
1374  nmg_vertexuse_nv(vu, norm_opp);
1375  }
1376  }
1377  }
1378 
1379  /* Compute "geometry" for region and shell */
1380  nmg_region_a(*r, tol);
1381 
1382  /* Release memory */
1383  /* All strips have vertices and normals */
1384  for (i=0; i<nstrips; i++) {
1385  bu_free((char *)strips[i].vp, "strip vertex[]");
1386  bu_free((char *)strips[i].norms, "strip norms[]");
1387  }
1388  /* All strips have faces, except for equator */
1389  for (i=0; i < nstrips; i++) {
1390  if (strips[i].fu == (struct faceuse **)0) continue;
1391  bu_free((char *)strips[i].fu, "strip faceuse[]");
1392  }
1393  bu_free((char *)strips, "strips[]");
1394  return 0;
1395  fail:
1396  /* Release memory */
1397  /* All strips have vertices and normals */
1398  for (i=0; i<nstrips; i++) {
1399  bu_free((char *)strips[i].vp, "strip vertex[]");
1400  bu_free((char *)strips[i].norms, "strip norms[]");
1401  }
1402  /* All strips have faces, except for equator */
1403  for (i=0; i < nstrips; i++) {
1404  if (strips[i].fu == (struct faceuse **)0) continue;
1405  bu_free((char *)strips[i].fu, "strip faceuse[]");
1406  }
1407  bu_free((char *)strips, "strips[]");
1408  return -1;
1409 }
1410 
1411 
1412 int
1413 rt_part_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
1414 {
1415  fastf_t maxrad, minrad;
1416  union record *rp;
1417  struct rt_part_internal *part;
1418 
1419  /* must be double for import and export */
1420  double v[ELEMENTS_PER_POINT];
1421  double h[ELEMENTS_PER_VECT];
1422  double vrad;
1423  double hrad;
1424 
1425  if (dbip) RT_CK_DBI(dbip);
1426 
1427  BU_CK_EXTERNAL(ep);
1428  rp = (union record *)ep->ext_buf;
1429  /* Check record type */
1430  if (rp->u_id != DBID_PARTICLE) {
1431  bu_log("rt_part_import4: defective record\n");
1432  return -1;
1433  }
1434 
1435  /* Convert from database to internal format */
1436  bu_cv_ntohd((unsigned char *)v, rp->part.p_v, ELEMENTS_PER_POINT);
1437  bu_cv_ntohd((unsigned char *)h, rp->part.p_h, ELEMENTS_PER_VECT);
1438  bu_cv_ntohd((unsigned char *)&vrad, rp->part.p_vrad, 1);
1439  bu_cv_ntohd((unsigned char *)&hrad, rp->part.p_hrad, 1);
1440 
1441  RT_CK_DB_INTERNAL(ip);
1442  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1443  ip->idb_type = ID_PARTICLE;
1444  ip->idb_meth = &OBJ[ID_PARTICLE];
1445  BU_ALLOC(ip->idb_ptr, struct rt_part_internal);
1446 
1447  part = (struct rt_part_internal *)ip->idb_ptr;
1448  part->part_magic = RT_PART_INTERNAL_MAGIC;
1449 
1450  /* Apply modeling transformations */
1451  if (mat == NULL) mat = bn_mat_identity;
1452  MAT4X3PNT(part->part_V, mat, v);
1453  MAT4X3VEC(part->part_H, mat, h);
1454  if ((part->part_vrad = vrad / mat[15]) < 0) {
1455  bu_log("unable to import particle, negative v radius\n");
1456  bu_free(ip->idb_ptr, "rt_part_internal");
1457  ip->idb_ptr=NULL;
1458  return -2;
1459  }
1460  if ((part->part_hrad = hrad / mat[15]) < 0) {
1461  bu_log("unable to import particle, negative h radius\n");
1462  bu_free(ip->idb_ptr, "rt_part_internal");
1463  ip->idb_ptr=NULL;
1464  return -3;
1465  }
1466 
1467  if (part->part_vrad > part->part_hrad) {
1468  maxrad = part->part_vrad;
1469  minrad = part->part_hrad;
1470  } else {
1471  maxrad = part->part_hrad;
1472  minrad = part->part_vrad;
1473  }
1474  if (maxrad <= 0) {
1475  bu_log("unable to import particle, negative radius\n");
1476  bu_free(ip->idb_ptr, "rt_part_internal");
1477  ip->idb_ptr=NULL;
1478  return -4;
1479  }
1480 
1481  if (MAGSQ(part->part_H) * 1000000 < maxrad * maxrad) {
1482  /* Height vector is insignificant, particle is a sphere */
1483  part->part_vrad = part->part_hrad = maxrad;
1484  VSETALL(part->part_H, 0); /* sanity */
1485  part->part_type = RT_PARTICLE_TYPE_SPHERE;
1486  return 0; /* OK */
1487  }
1488 
1489  if ((maxrad - minrad) / maxrad < 0.001) {
1490  /* radii are nearly equal, particle is a cylinder (lozenge) */
1491  part->part_vrad = part->part_hrad = maxrad;
1492  part->part_type = RT_PARTICLE_TYPE_CYLINDER;
1493  return 0; /* OK */
1494  }
1495 
1496  part->part_type = RT_PARTICLE_TYPE_CONE;
1497  return 0; /* OK */
1498 }
1499 
1500 
1501 int
1502 rt_part_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1503 {
1504  struct rt_part_internal *pip;
1505  union record *rec;
1506 
1507  /* must be double for import and export */
1508  double vert[ELEMENTS_PER_POINT];
1509  double hi[ELEMENTS_PER_VECT];
1510  double vrad;
1511  double hrad;
1512 
1513  if (dbip) RT_CK_DBI(dbip);
1514 
1515  RT_CK_DB_INTERNAL(ip);
1516  if (ip->idb_type != ID_PARTICLE) return -1;
1517  pip = (struct rt_part_internal *)ip->idb_ptr;
1518  RT_PART_CK_MAGIC(pip);
1519 
1520  BU_CK_EXTERNAL(ep);
1521  ep->ext_nbytes = sizeof(union record);
1522  ep->ext_buf = (uint8_t*)bu_calloc(1, ep->ext_nbytes, "part external");
1523  rec = (union record *)ep->ext_buf;
1524 
1525  /* Convert from user units to mm */
1526  VSCALE(vert, pip->part_V, local2mm);
1527  VSCALE(hi, pip->part_H, local2mm);
1528  vrad = pip->part_vrad * local2mm;
1529  hrad = pip->part_hrad * local2mm;
1530  /* pip->part_type is not converted -- internal only */
1531 
1532  rec->part.p_id = DBID_PARTICLE;
1533  bu_cv_htond(rec->part.p_v, (unsigned char *)vert, ELEMENTS_PER_POINT);
1534  bu_cv_htond(rec->part.p_h, (unsigned char *)hi, ELEMENTS_PER_VECT);
1535  bu_cv_htond(rec->part.p_vrad, (unsigned char *)&vrad, 1);
1536  bu_cv_htond(rec->part.p_hrad, (unsigned char *)&hrad, 1);
1537 
1538  return 0;
1539 }
1540 
1541 
1542 int
1543 rt_part_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
1544 {
1545  fastf_t maxrad, minrad;
1546  struct rt_part_internal *part;
1547 
1548  /* must be double for import and export */
1549  double vec[8];
1550 
1551  if (dbip) RT_CK_DBI(dbip);
1552 
1553  BU_CK_EXTERNAL(ep);
1554 
1556 
1557  RT_CK_DB_INTERNAL(ip);
1558  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1559  ip->idb_type = ID_PARTICLE;
1560  ip->idb_meth = &OBJ[ID_PARTICLE];
1561  BU_ALLOC(ip->idb_ptr, struct rt_part_internal);
1562 
1563  part = (struct rt_part_internal *)ip->idb_ptr;
1564  part->part_magic = RT_PART_INTERNAL_MAGIC;
1565 
1566  /* Convert from database (network) to internal (host) format */
1567  bu_cv_ntohd((unsigned char *)vec, ep->ext_buf, 8);
1568 
1569  /* Apply modeling transformations */
1570  if (mat == NULL) mat = bn_mat_identity;
1571  MAT4X3PNT(part->part_V, mat, &vec[0*3]);
1572  MAT4X3VEC(part->part_H, mat, &vec[1*3]);
1573  if ((part->part_vrad = vec[2*3] / mat[15]) < 0) {
1574  bu_free(ip->idb_ptr, "rt_part_internal");
1575  ip->idb_ptr=NULL;
1576  bu_log("unable to import particle, negative v radius\n");
1577  return -2;
1578  }
1579  if ((part->part_hrad = vec[2*3+1] / mat[15]) < 0) {
1580  bu_free(ip->idb_ptr, "rt_part_internal");
1581  ip->idb_ptr=NULL;
1582  bu_log("unable to import particle, negative h radius\n");
1583  return -3;
1584  }
1585 
1586  if (part->part_vrad > part->part_hrad) {
1587  maxrad = part->part_vrad;
1588  minrad = part->part_hrad;
1589  } else {
1590  maxrad = part->part_hrad;
1591  minrad = part->part_vrad;
1592  }
1593  if (maxrad <= 0) {
1594  bu_free(ip->idb_ptr, "rt_part_internal");
1595  ip->idb_ptr=NULL;
1596  bu_log("unable to import particle, negative radius\n");
1597  return -4;
1598  }
1599 
1600  if (MAGSQ(part->part_H) * 1000000 < maxrad * maxrad) {
1601  /* Height vector is insignificant, particle is a sphere */
1602  part->part_vrad = part->part_hrad = maxrad;
1603  VSETALL(part->part_H, 0); /* sanity */
1604  part->part_type = RT_PARTICLE_TYPE_SPHERE;
1605  return 0; /* OK */
1606  }
1607 
1608  if ((maxrad - minrad) / maxrad < 0.001) {
1609  /* radii are nearly equal, particle is a cylinder (lozenge) */
1610  part->part_vrad = part->part_hrad = maxrad;
1611  part->part_type = RT_PARTICLE_TYPE_CYLINDER;
1612  return 0; /* OK */
1613  }
1614 
1615  part->part_type = RT_PARTICLE_TYPE_CONE;
1616  return 0; /* OK */
1617 }
1618 
1619 
1620 int
1621 rt_part_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1622 {
1623  struct rt_part_internal *pip;
1624 
1625  /* must be double for import and export */
1626  double vec[8];
1627 
1628  if (dbip) RT_CK_DBI(dbip);
1629 
1630  RT_CK_DB_INTERNAL(ip);
1631  if (ip->idb_type != ID_PARTICLE) return -1;
1632  pip = (struct rt_part_internal *)ip->idb_ptr;
1633  RT_PART_CK_MAGIC(pip);
1634 
1635  BU_CK_EXTERNAL(ep);
1637  ep->ext_buf = (uint8_t*)bu_malloc(ep->ext_nbytes, "part external");
1638 
1639  /* scale 'em into local buffer */
1640  VSCALE(&vec[0*3], pip->part_V, local2mm);
1641  VSCALE(&vec[1*3], pip->part_H, local2mm);
1642  vec[2*3] = pip->part_vrad * local2mm;
1643  vec[2*3+1] = pip->part_hrad * local2mm;
1644 
1645  /* !!! should make sure the values are proper (no negative radius) */
1646 
1647  /* Convert from internal (host) to database (network) format */
1648  bu_cv_htond(ep->ext_buf, (unsigned char *)vec, 8);
1649 
1650  return 0;
1651 }
1652 
1653 
1654 /**
1655  * Make human-readable formatted presentation of this solid.
1656  * First line describes type of solid.
1657  * Additional lines are indented one tab, and give parameter values.
1658  */
1659 int
1660 rt_part_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1661 {
1662  register struct rt_part_internal *pip =
1663  (struct rt_part_internal *)ip->idb_ptr;
1664  char buf[256];
1665 
1666  RT_PART_CK_MAGIC(pip);
1667  switch (pip->part_type) {
1668  case RT_PARTICLE_TYPE_SPHERE:
1669  bu_vls_strcat(str, "spherical particle\n");
1670  if (!verbose)
1671  return 0;
1672  sprintf(buf, "\tV (%g, %g, %g)\n",
1673  INTCLAMP(pip->part_V[X] * mm2local),
1674  INTCLAMP(pip->part_V[Y] * mm2local),
1675  INTCLAMP(pip->part_V[Z] * mm2local));
1676  bu_vls_strcat(str, buf);
1677  sprintf(buf, "\tradius = %g\n",
1678  INTCLAMP(pip->part_vrad * mm2local));
1679  bu_vls_strcat(str, buf);
1680 
1681  break;
1682  case RT_PARTICLE_TYPE_CYLINDER:
1683  bu_vls_strcat(str, "cylindrical particle (lozenge)\n");
1684  if (!verbose)
1685  return 0;
1686  sprintf(buf, "\tV (%g, %g, %g)\n",
1687  INTCLAMP(pip->part_V[X] * mm2local),
1688  INTCLAMP(pip->part_V[Y] * mm2local),
1689  INTCLAMP(pip->part_V[Z] * mm2local));
1690  bu_vls_strcat(str, buf);
1691  sprintf(buf, "\tH (%g, %g, %g)\n",
1692  INTCLAMP(pip->part_H[X] * mm2local),
1693  INTCLAMP(pip->part_H[Y] * mm2local),
1694  INTCLAMP(pip->part_H[Z] * mm2local));
1695  bu_vls_strcat(str, buf);
1696  sprintf(buf, "\tradius = %g\n",
1697  INTCLAMP(pip->part_vrad * mm2local));
1698  bu_vls_strcat(str, buf);
1699 
1700  break;
1701  case RT_PARTICLE_TYPE_CONE:
1702  bu_vls_strcat(str, "conical particle\n");
1703  if (!verbose)
1704  return 0;
1705  sprintf(buf, "\tV (%g, %g, %g)\n",
1706  INTCLAMP(pip->part_V[X] * mm2local),
1707  INTCLAMP(pip->part_V[Y] * mm2local),
1708  INTCLAMP(pip->part_V[Z] * mm2local));
1709  bu_vls_strcat(str, buf);
1710  sprintf(buf, "\tH (%g, %g, %g)\n",
1711  INTCLAMP(pip->part_H[X] * mm2local),
1712  INTCLAMP(pip->part_H[Y] * mm2local),
1713  INTCLAMP(pip->part_H[Z] * mm2local));
1714  bu_vls_strcat(str, buf);
1715  sprintf(buf, "\tv end radius = %g\n",
1716  INTCLAMP(pip->part_vrad * mm2local));
1717  bu_vls_strcat(str, buf);
1718  sprintf(buf, "\th end radius = %g\n",
1719  INTCLAMP(pip->part_hrad * mm2local));
1720  bu_vls_strcat(str, buf);
1721  break;
1722  default:
1723  bu_vls_strcat(str, "Unknown particle type\n");
1724  return -1;
1725  }
1726  return 0;
1727 }
1728 
1729 
1730 /**
1731  * Free the storage associated with the rt_db_internal version of this solid.
1732  */
1733 void
1735 {
1736  RT_CK_DB_INTERNAL(ip);
1737 
1738  bu_free(ip->idb_ptr, "particle ifree");
1739  ip->idb_ptr = ((void *)0);
1740 }
1741 
1742 
1743 int
1744 rt_part_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
1745 {
1746  if (ip) RT_CK_DB_INTERNAL(ip);
1747 
1748  return 0; /* OK */
1749 }
1750 
1751 
1752 void
1753 rt_part_volume(fastf_t *vol, const struct rt_db_internal *ip)
1754 {
1755  fastf_t vrad, hrad, mag_h;
1756  struct rt_part_internal *pip = (struct rt_part_internal *)ip->idb_ptr;
1757  RT_PART_CK_MAGIC(pip);
1758 
1759  vrad = pip->part_vrad;
1760  hrad = pip->part_hrad;
1761  mag_h = MAGNITUDE(pip->part_H);
1762 
1763  if (EQUAL(vrad, hrad)) {
1764  *vol = M_PI * vrad * vrad * (4.0/3.0 * vrad + mag_h);
1765  } else {
1766  fastf_t vrad3, hrad3, mid_section;
1767  vrad3 = vrad * vrad * vrad;
1768  hrad3 = hrad * hrad * hrad;
1769  mid_section = M_PI * mag_h * (hrad * hrad + vrad * vrad + hrad * vrad) / 3.0;
1770  *vol = 2.0/3.0 * M_PI * (vrad3 + hrad3) + mid_section;
1771  }
1772 }
1773 
1774 
1775 void
1776 rt_part_surf_area(fastf_t *area, const struct rt_db_internal *ip)
1777 {
1778  fastf_t vrad, hrad, mag_h;
1779  struct rt_part_internal *pip = (struct rt_part_internal *)ip->idb_ptr;
1780  RT_PART_CK_MAGIC(pip);
1781 
1782  vrad = pip->part_vrad;
1783  hrad = pip->part_hrad;
1784  mag_h = MAGNITUDE(pip->part_H);
1785 
1786  if (EQUAL(vrad, hrad)) {
1787  *area = M_2PI * vrad * (2.0 * vrad + mag_h);
1788  } else {
1789  fastf_t mid_section;
1790  mid_section = M_PI * ((vrad + hrad) * sqrt((vrad - hrad) * (vrad - hrad) + mag_h * mag_h));
1791  *area = M_2PI * (vrad * vrad + hrad * hrad) + mid_section;
1792  }
1793 }
1794 
1795 
1796 /*
1797  * Local Variables:
1798  * mode: C
1799  * tab-width: 8
1800  * indent-tabs-mode: t
1801  * c-file-style: "stroustrup"
1802  * End:
1803  * ex: shiftwidth=4 tabstop=8
1804  */
mat_t lower_invRoS
Definition: part.c:1031
int rt_part_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: part.c:921
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
void rt_part_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
Definition: part.c:835
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
fastf_t part_v_hdist
Definition: part.c:203
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
mat_t part_invRoS
Definition: part.c:199
#define SMALL
Definition: defines.h:351
#define RT_CK_RTI(_p)
Definition: raytrace.h:1833
mat_t upper_invRinvS
Definition: part.c:1028
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
int rt_part_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
Definition: part.c:449
if lu s
Definition: nmg_mod.c:3860
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
#define VSETALL(a, s)
Definition: color.c:54
Definition: raytrace.h:215
#define M_PI
Definition: fft.h:35
Definition: pc.h:108
void rt_hitsort(struct hit h[], int nh)
Definition: raytrace.h:368
fastf_t part_v_erad
Definition: part.c:205
#define BU_ASSERT_LONG(_lhs, _relation, _rhs)
Definition: defines.h:240
Definition: raytrace.h:248
void nmg_vertex_gv(struct vertex *v, const fastf_t *pt)
Definition: nmg_mk.c:1668
#define ID_PARTICLE
Particle system solid.
Definition: raytrace.h:474
int rt_part_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: part.c:1744
struct faceuse ** fu
Definition: part.c:1042
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
int rt_part_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: part.c:227
Header file for the BRL-CAD common definitions.
void nmg_vertexuse_nv(struct vertexuse *vu, const fastf_t *norm)
Definition: nmg_mk.c:1719
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)
#define HIDDEN
Definition: common.h:86
struct bu_list l
Definition: raytrace.h:369
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
#define RT_PARTICLE_SURF_VSPHERE
Definition: part.c:211
int rt_part_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: part.c:1660
if(share_geom)
Definition: nmg_mod.c:3829
#define RT_PARTICLE_SURF_BODY
Definition: part.c:212
int idb_major_type
Definition: raytrace.h:192
int rt_part_import4(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: part.c:1413
Definition: color.c:49
vect_t hit_vpriv
PRIVATE vector for xxx_*()
Definition: raytrace.h:253
#define RT_ADD_VLIST(hd, pnt, draw)
Definition: raytrace.h:1865
void rt_part_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: part.c:1776
void bn_mat_print(const char *title, const mat_t m)
Definition: mat.c:81
#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
mat_t lower_invRinvS
Definition: part.c:1030
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
void rt_part_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
Definition: part.c:740
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
fastf_t a_diverge
slope of beam divergence/mm
Definition: raytrace.h:1600
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
struct rt_part_internal part_int
Definition: part.c:197
fastf_t uv_dv
delta in v
Definition: raytrace.h:344
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)
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
mat_t part_SoR
Definition: part.c:198
void bn_mat_fromto(mat_t m, const fastf_t *from, const fastf_t *to, const struct bn_tol *tol)
Definition: mat.c:639
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
void rt_part_volume(fastf_t *vol, const struct rt_db_internal *ip)
Definition: part.c:1753
const struct bu_structparse rt_part_parse[]
Definition: part.c:215
int rt_part_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: part.c:1621
void bn_mat_trn(mat_t om, const mat_t im)
#define UNUSED(parameter)
Definition: common.h:239
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
mat_t upper_invRoS
Definition: part.c:1029
goto out
Definition: nmg_mod.c:3846
void bn_mat_mul(mat_t o, const mat_t a, const mat_t b)
struct faceuse * nmg_find_fu_of_vu(const struct vertexuse *vu)
Definition: nmg_info.c:304
Support for uniform tolerances.
Definition: tol.h:71
double bn_atan2(double x, double y)
Definition: mat.c:104
HIDDEN void rt_part_hemisphere(register point_t(*ov), register fastf_t *v, fastf_t *a, fastf_t *b, fastf_t *h)
Definition: part.c:895
#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
void rt_part_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
Definition: part.c:789
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
struct shell * s
Definition: part.c:1027
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
#define bu_offsetof(_t, _m)
Definition: parse.h:64
int rt_part_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
Definition: part.c:1543
#define RT_CK_DBI(_p)
Definition: raytrace.h:829
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
int rt_part_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: part.c:1055
#define ZERO(val)
Definition: units.c:38
void * idb_ptr
Definition: raytrace.h:195
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
point_t st_min
min X, Y, Z of bounding RPP
Definition: raytrace.h:437
int rt_part_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: part.c:1502
void bu_cv_ntohd(unsigned char *out, const unsigned char *in, size_t count)
const struct rt_functab OBJ[]
Definition: table.c:159
void rt_part_ifree(struct rt_db_internal *ip)
Definition: part.c:1734
fastf_t part_hrad_prime
Definition: part.c:201
void bn_vec_ortho(vect_t out, const vect_t in)
#define R
Definition: msr.c:55
#define RT_PARTICLE_SURF_HSPHERE
Definition: part.c:213
fastf_t part_vrad_prime
Definition: part.c:200
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t uv_du
delta in u
Definition: raytrace.h:343
void rt_part_free(register struct soltab *stp)
Definition: part.c:874
fastf_t part_h_hdist
Definition: part.c:204
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
fastf_t theta_tol
Definition: part.c:1032
Definition: color.c:51
double norm
normal tol
Definition: raytrace.h:182
int hit_surfno
solid-specific surface indicator
Definition: raytrace.h:255
void bn_vec_perp(vect_t new_vec, const vect_t old_vec)
Definition: mat.c:616
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
fastf_t a_rbeam
initial beam radius (mm)
Definition: raytrace.h:1599
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
vect_t * norms
Definition: part.c:1040
#define BU_CK_EXTERNAL(_p)
Definition: parse.h:224
struct faceuse * nmg_cmface(struct shell *s, struct vertex ***verts, int n)
Definition: nmg_mod.c:979
void rt_ell_16pts(fastf_t *ov, fastf_t *V, fastf_t *A, fastf_t *B)
Definition: ell.c:606
size_t ext_nbytes
Definition: parse.h:210
fastf_t part_h_erad
Definition: part.c:206
#define RT_PART_INTERNAL_MAGIC
Definition: magic.h:102
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
struct vertex ** vp
Definition: part.c:1039
Definition: vls.h:56
#define M_SQRT1_2
Definition: fft.h:38
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
int rt_part_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: part.c:274
int nverts_per_strip
Definition: part.c:1037
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
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
#define RT_HIT_MAGIC
Definition: magic.h:161
void rt_part_print(register const struct soltab *stp)
Definition: part.c:408
void nmg_region_a(struct nmgregion *r, const struct bn_tol *tol)
Definition: nmg_mk.c:2557