BRL-CAD
rhc.c
Go to the documentation of this file.
1 /* R H C . 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/rhc/rhc.c
23  *
24  * Intersect a ray with a Right Hyperbolic Cylinder.
25  *
26  * Algorithm -
27  *
28  * Given V, H, R, and B, there is a set of points on this rhc
29  *
30  * { (x, y, z) | (x, y, z) is on rhc }
31  *
32  * Through a series of Affine Transformations, this set of points will
33  * be transformed into a set of points on an rhc located at the origin
34  * with a rectangular halfwidth R of 1 along the Y axis, a height H of
35  * +1 along the -X axis, a distance B of 1 along the -Z axis between
36  * the vertex V and the tip of the hyperbola, and a distance c between
37  * the tip of the hyperbola and the vertex of the asymptotic cone.
38  *
39  *
40  * { (x', y', z') | (x', y', z') is on rhc at origin }
41  *
42  * The transformation from X to X' is accomplished by:
43  *
44  * X' = S(R(X - V))
45  *
46  * where R(X) = (H/(-|H|))
47  * (R/(|R|)) . X
48  * (B/(-|B|))
49  *
50  * and S(X) = (1/|H| 0 0)
51  * (0 1/|R| 0) . X
52  * (0 0 1/|B|)
53  *
54  * To find the intersection of a line with the surface of the rhc,
55  * consider the parametric line L:
56  *
57  * L : { P(n) | P + t(n) . D }
58  *
59  * Call W the actual point of intersection between L and the rhc.
60  * Let W' be the point of intersection between L' and the unit rhc.
61  *
62  * L' : { P'(n) | P' + t(n) . D' }
63  *
64  * W = invR(invS(W')) + V
65  *
66  * Where W' = k D' + P'.
67  *
68  * If Dy' and Dz' are both 0, then there is no hit on the rhc; but the
69  * end plates need checking. If there is now only 1 hit point, the
70  * top plate needs to be checked as well.
71  *
72  * Line L' hits the infinitely long canonical rhc at W' when
73  *
74  * A * k**2 + B * k + C = 0
75  *
76  * where
77  *
78  * A = Dz'**2 - Dy'**2 * (1 + 2*c')
79  * B = 2 * ((1 + c' + Pz') * Dz' - (1 + 2*c') * Dy' * Py'
80  * C = (Pz' + c' + 1)**2 - (1 + 2*c') * Py'**2 - c'**2
81  * b = |Breadth| = 1.0
82  * h = |Height| = 1.0
83  * r = 1.0
84  * c' = c / |Breadth|
85  *
86  * The quadratic formula yields k (which is constant):
87  *
88  * k = [ -B +/- sqrt(B**2 - 4 * A * C)] / (2 * A)
89  *
90  * Now, D' = S(R(D))
91  * and P' = S(R(P - V))
92  *
93  * Substituting,
94  *
95  * W = V + invR(invS[ k *(S(R(D))) + S(R(P - V)) ])
96  * = V + invR((k * R(D)) + R(P - V))
97  * = V + k * D + P - V
98  * = k * D + P
99  *
100  * Note that ``k'' is constant, and is the same in the formulations
101  * for both W and W'.
102  *
103  * The hit at ``k'' is a hit on the canonical rhc IFF
104  * -1 <= Wx' <= 0 and -1 <= Wz' <= 0.
105  *
106  * NORMALS. Given the point W on the surface of the rhc, what is the
107  * vector normal to the tangent plane at that point?
108  *
109  * Map W onto the unit rhc, i.e.: W' = S(R(W - V)).
110  *
111  * Plane on unit rhc at W' has a normal vector N' where
112  *
113  * N' = <0, Wy'*(1 + 2*c), -z-c-1>.
114  *
115  * The plane transforms back to the tangent plane at W, and this new
116  * plane (on the original rhc) has a normal vector of N, viz:
117  *
118  * N = inverse[ transpose(inverse[ S o R ]) ] (N')
119  *
120  * because if H is perpendicular to plane Q, and matrix M maps from
121  * Q to Q', then inverse[ transpose(M) ] (H) is perpendicular to Q'.
122  * Here, H and Q are in "prime space" with the unit sphere.
123  * [Somehow, the notation here is backwards].
124  * So, the mapping matrix M = inverse(S o R), because
125  * S o R maps from normal space to the unit sphere.
126  *
127  * N = inverse[ transpose(inverse[ S o R ]) ] (N')
128  * = inverse[ transpose(invR o invS) ] (N')
129  * = inverse[ transpose(invS) o transpose(invR) ] (N')
130  * = inverse[ inverse(S) o R ] (N')
131  * = invR o S (N')
132  *
133  * because inverse(R) = transpose(R), so R = transpose(invR),
134  * and S = transpose(S).
135  *
136  * Note that the normal vector produced above will not have unit
137  * length.
138  *
139  * THE TOP AND END PLATES.
140  *
141  * If Dz' == 0, line L' is parallel to the top plate, so there is no
142  * hit on the top plate. Otherwise, rays intersect the top plate
143  * with k = (0 - Pz')/Dz'. The solution is within the top plate
144  * IFF -1 <= Wx' <= 0 and -1 <= Wy' <= 1.
145  *
146  * If Dx' == 0, line L' is parallel to the end plates, so there is no
147  * hit on the end plates. Otherwise, rays intersect the front plate
148  * line L' hits the front plate with k = (0 - Px') / Dx', and and hits
149  * the back plate with k = (-1 - Px') / Dx'.
150  *
151  * The solution W' is within an end plate IFF
152  *
153  * (Wz' + c + 1)**2 - (Wy'**2 * (2*c + 1) >= c**2 and Wz' <= 1.0
154  *
155  * The normal for a hit on the top plate is -Bunit.
156  * The normal for a hit on the front plate is -Hunit, and
157  * the normal for a hit on the back plate is +Hunit.
158  *
159  */
160 /** @} */
161 
162 #include "common.h"
163 
164 #include <stddef.h>
165 #include <string.h>
166 #include <math.h>
167 #include "bio.h"
168 
169 #include "bu/cv.h"
170 #include "vmath.h"
171 #include "db.h"
172 #include "nmg.h"
173 #include "rtgeom.h"
174 #include "raytrace.h"
175 
176 #include "../../librt_private.h"
177 
178 static int
179 rhc_is_valid(struct rt_rhc_internal *rhc);
180 
181 struct rhc_specific {
182  point_t rhc_V; /* vector to rhc origin */
183  vect_t rhc_Bunit; /* unit B vector */
184  vect_t rhc_Hunit; /* unit H vector */
185  vect_t rhc_Runit; /* unit vector, B x H */
186  mat_t rhc_SoR; /* Scale(Rot(vect)) */
187  mat_t rhc_invRoS; /* invRot(Scale(vect)) */
188  fastf_t rhc_b; /* |B| */
189  fastf_t rhc_c; /* c */
190  fastf_t rhc_cprime; /* c / |B| */
191  fastf_t rhc_rsq; /* r * r */
192 };
193 
194 
195 const struct bu_structparse rt_rhc_parse[] = {
196  { "%f", 3, "V", bu_offsetofarray(struct rt_rhc_internal, rhc_V, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
197  { "%f", 3, "H", bu_offsetofarray(struct rt_rhc_internal, rhc_H, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
198  { "%f", 3, "B", bu_offsetofarray(struct rt_rhc_internal, rhc_B, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
199  { "%f", 1, "r", bu_offsetof(struct rt_rhc_internal, rhc_r), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
200  { "%f", 1, "c", bu_offsetof(struct rt_rhc_internal, rhc_c), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
201  { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
202 };
203 
204 
205 /**
206  * Calculate the bounding RPP for an RHC
207  */
208 int
209 rt_rhc_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol))
210 {
211 
212  struct rt_rhc_internal *xip;
213  vect_t rinv, rvect, rv2, working;
214 
215  RT_CK_DB_INTERNAL(ip);
216  xip = (struct rt_rhc_internal *)ip->idb_ptr;
217  RT_RHC_CK_MAGIC(xip);
218 
219  VSETALL((*min), INFINITY);
220  VSETALL((*max), -INFINITY);
221 
222  VCROSS(rvect, xip->rhc_H, xip->rhc_B);
223  VREVERSE(rinv, rvect);
224  VUNITIZE(rvect);
225  VUNITIZE(rinv);
226  VSCALE(rvect, rvect, xip->rhc_r);
227  VSCALE(rinv, rinv, xip->rhc_r);
228 
229  VADD2(working, xip->rhc_V, rvect);
230  VMINMAX((*min), (*max), working);
231 
232  VADD2(working, xip->rhc_V, rinv);
233  VMINMAX((*min), (*max), working);
234 
235  VADD3(working, xip->rhc_V, rvect, xip->rhc_H);
236  VMINMAX((*min), (*max), working);
237 
238  VADD3(working, xip->rhc_V, rinv, xip->rhc_H);
239  VMINMAX((*min), (*max), working);
240 
241  VADD2(rv2, xip->rhc_V, xip->rhc_B);
242 
243  VADD2(working, rv2, rvect);
244  VMINMAX((*min), (*max), working);
245 
246  VADD2(working, rv2, rinv);
247  VMINMAX((*min), (*max), working);
248 
249  VADD3(working, rv2, rvect, xip->rhc_H);
250  VMINMAX((*min), (*max), working);
251 
252  VADD3(working, rv2, rinv, xip->rhc_H);
253  VMINMAX((*min), (*max), working);
254 
255  return 0;
256 }
257 
258 
259 /**
260  * Given a pointer to a GED database record, and a transformation matrix,
261  * determine if this is a valid RHC, and if so, precompute various
262  * terms of the formula.
263  *
264  * Returns -
265  * 0 RHC is OK
266  * !0 Error in description
267  *
268  * Implicit return -
269  * A struct rhc_specific is created, and its address is stored in
270  * stp->st_specific for use by rhc_shot().
271  */
272 int
273 rt_rhc_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
274 {
275  struct rt_rhc_internal *xip;
276  struct rhc_specific *rhc;
277  fastf_t magsq_b, magsq_h, magsq_r;
278  fastf_t mag_b, mag_h, mag_r;
279  mat_t R;
280  mat_t Rinv;
281  mat_t S;
282  vect_t invsq; /* [ 1/(|H|**2), 1/(|R|**2), 1/(|B|**2) ] */
283 
284  RT_CK_DB_INTERNAL(ip);
285  RT_CK_RTI(rtip);
286 
287  xip = (struct rt_rhc_internal *)ip->idb_ptr;
288  if (!rhc_is_valid(xip)) {
289  return 1;
290  }
291 
292  /* compute |B| |H| */
293  mag_b = sqrt(magsq_b = MAGSQ(xip->rhc_B));
294  mag_h = sqrt(magsq_h = MAGSQ(xip->rhc_H));
295  mag_r = xip->rhc_r;
296  magsq_r = mag_r * mag_r;
297 
298  stp->st_id = ID_RHC; /* set soltab ID */
299  stp->st_meth = &OBJ[ID_RHC];
300 
301  BU_GET(rhc, struct rhc_specific);
302  stp->st_specific = (void *)rhc;
303  rhc->rhc_b = mag_b;
304  rhc->rhc_rsq = magsq_r;
305  rhc->rhc_c = xip->rhc_c;
306 
307  /* make unit vectors in B, H, and BxH directions */
308  VMOVE(rhc->rhc_Hunit, xip->rhc_H);
309  VUNITIZE(rhc->rhc_Hunit);
310  VMOVE(rhc->rhc_Bunit, xip->rhc_B);
311  VUNITIZE(rhc->rhc_Bunit);
312  VCROSS(rhc->rhc_Runit, rhc->rhc_Bunit, rhc->rhc_Hunit);
313 
314  VMOVE(rhc->rhc_V, xip->rhc_V);
315  rhc->rhc_cprime = xip->rhc_c / mag_b;
316 
317  /* Compute R and Rinv matrices */
318  MAT_IDN(R);
319  VREVERSE(&R[0], rhc->rhc_Hunit);
320  VMOVE(&R[4], rhc->rhc_Runit);
321  VREVERSE(&R[8], rhc->rhc_Bunit);
322  bn_mat_trn(Rinv, R); /* inv of rot mat is trn */
323 
324  /* Compute S */
325  VSET(invsq, 1.0 / magsq_h, 1.0 / magsq_r, 1.0 / magsq_b);
326  MAT_IDN(S);
327  S[ 0] = sqrt(invsq[0]);
328  S[ 5] = sqrt(invsq[1]);
329  S[10] = sqrt(invsq[2]);
330 
331  /* Compute SoR and invRoS */
332  bn_mat_mul(rhc->rhc_SoR, S, R);
333  bn_mat_mul(rhc->rhc_invRoS, Rinv, S);
334 
335  /* Compute bounding sphere and RPP */
336  /* bounding sphere center */
337  VJOIN2(stp->st_center, rhc->rhc_V,
338  mag_h / 2.0, rhc->rhc_Hunit,
339  mag_b / 2.0, rhc->rhc_Bunit);
340  /* bounding radius */
341  stp->st_bradius = 0.5 * sqrt(magsq_h + 4.0 * magsq_r + magsq_b);
342  /* approximate bounding radius */
343  stp->st_aradius = stp->st_bradius;
344 
345  if (rt_rhc_bbox(ip, &(stp->st_min), &(stp->st_max), &rtip->rti_tol)) {
346  return 1;
347  }
348 
349  return 0; /* OK */
350 }
351 
352 
353 void
354 rt_rhc_print(const struct soltab *stp)
355 {
356  const struct rhc_specific *rhc =
357  (struct rhc_specific *)stp->st_specific;
358 
359  VPRINT("V", rhc->rhc_V);
360  VPRINT("Bunit", rhc->rhc_Bunit);
361  VPRINT("Hunit", rhc->rhc_Hunit);
362  VPRINT("Runit", rhc->rhc_Runit);
363  bn_mat_print("S o R", rhc->rhc_SoR);
364  bn_mat_print("invR o S", rhc->rhc_invRoS);
365 }
366 
367 
368 /* hit_surfno is set to one of these */
369 #define RHC_NORM_BODY (1) /* compute normal */
370 #define RHC_NORM_TOP (2) /* copy rhc_N */
371 #define RHC_NORM_FRT (3) /* copy reverse rhc_N */
372 #define RHC_NORM_BACK (4)
373 
374 /**
375  * Intersect a ray with a rhc.
376  * If an intersection occurs, a struct seg will be acquired
377  * and filled in.
378  *
379  * Returns -
380  * 0 MISS
381  * >0 HIT
382  */
383 int
384 rt_rhc_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
385 {
386  struct rhc_specific *rhc =
387  (struct rhc_specific *)stp->st_specific;
388  vect_t dprime; /* D' */
389  vect_t pprime; /* P' */
390  fastf_t k1, k2; /* distance constants of solution */
391  fastf_t x;
392  vect_t xlated; /* translated vector */
393  struct hit hits[3]; /* 2 potential hit points */
394  struct hit *hitp; /* pointer to hit point */
395 
396  hitp = &hits[0];
397 
398  /* out, Mat, vect */
399  MAT4X3VEC(dprime, rhc->rhc_SoR, rp->r_dir);
400  VSUB2(xlated, rp->r_pt, rhc->rhc_V);
401  MAT4X3VEC(pprime, rhc->rhc_SoR, xlated);
402 
403  x = rhc->rhc_cprime;
404 
405  if (ZERO(dprime[Y]) && ZERO(dprime[Z])) {
406  goto check_plates;
407  }
408 
409  /* Find roots of the equation, using formula for quadratic */
410  {
411  fastf_t a, b, c; /* coeffs of polynomial */
412  fastf_t disc; /* disc of radical */
413 
414  a = dprime[Z] * dprime[Z] - dprime[Y] * dprime[Y] * (1 + 2 * x);
415  b = 2 * ((pprime[Z] + x + 1) * dprime[Z]
416  - (2 * x + 1) * dprime[Y] * pprime[Y]);
417  c = (pprime[Z] + x + 1) * (pprime[Z] + x + 1)
418  - (2 * x + 1) * pprime[Y] * pprime[Y] - x * x;
419 
420  if (!NEAR_ZERO(a, RT_PCOEF_TOL)) {
421  disc = b * b - 4 * a * c;
422 
423  if (disc <= 0) {
424  goto check_plates;
425  }
426 
427  disc = sqrt(disc);
428 
429  k1 = (-b + disc) / (2.0 * a);
430  k2 = (-b - disc) / (2.0 * a);
431 
432  /*
433  * k1 and k2 are potential solutions to intersection with side.
434  * See if they fall in range.
435  */
436  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
437 
438  if (hitp->hit_vpriv[X] >= -1.0
439  && hitp->hit_vpriv[X] <= 0.0
440  && hitp->hit_vpriv[Z] >= -1.0
441  && hitp->hit_vpriv[Z] <= 0.0) {
442  hitp->hit_magic = RT_HIT_MAGIC;
443  hitp->hit_dist = k1;
444  hitp->hit_surfno = RHC_NORM_BODY; /* compute N */
445  hitp++;
446  }
447 
448  VJOIN1(hitp->hit_vpriv, pprime, k2, dprime); /* hit' */
449 
450  if (hitp->hit_vpriv[X] >= -1.0
451  && hitp->hit_vpriv[X] <= 0.0
452  && hitp->hit_vpriv[Z] >= -1.0
453  && hitp->hit_vpriv[Z] <= 0.0) {
454  hitp->hit_magic = RT_HIT_MAGIC;
455  hitp->hit_dist = k2;
456  hitp->hit_surfno = RHC_NORM_BODY; /* compute N */
457  hitp++;
458  }
459  } else if (!NEAR_ZERO(b, RT_PCOEF_TOL)) {
460  k1 = -c / b;
461  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
462 
463  if (hitp->hit_vpriv[X] >= -1.0
464  && hitp->hit_vpriv[X] <= 0.0
465  && hitp->hit_vpriv[Z] >= -1.0
466  && hitp->hit_vpriv[Z] <= 0.0) {
467  hitp->hit_magic = RT_HIT_MAGIC;
468  hitp->hit_dist = k1;
469  hitp->hit_surfno = RHC_NORM_BODY; /* compute N */
470  hitp++;
471  }
472  }
473  }
474 
475 
476  /*
477  * Check for hitting the top and end plates.
478  */
479 check_plates:
480 
481  /* check front and back plates */
482  if (hitp < &hits[2] && !ZERO(dprime[X])) {
483  /* 0 or 1 hits so far, this is worthwhile */
484  k1 = -pprime[X] / dprime[X]; /* front plate */
485  k2 = (-1.0 - pprime[X]) / dprime[X]; /* back plate */
486 
487  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
488 
489  if ((hitp->hit_vpriv[Z] + x + 1.0)
490  * (hitp->hit_vpriv[Z] + x + 1.0)
491  - hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y]
492  * (1.0 + 2 * x) >= x * x
493  && hitp->hit_vpriv[Z] >= -1.0
494  && hitp->hit_vpriv[Z] <= 0.0) {
495  hitp->hit_magic = RT_HIT_MAGIC;
496  hitp->hit_dist = k1;
497  hitp->hit_surfno = RHC_NORM_FRT; /* -H */
498  hitp++;
499  }
500 
501  VJOIN1(hitp->hit_vpriv, pprime, k2, dprime); /* hit' */
502 
503  if ((hitp->hit_vpriv[Z] + x + 1.0)
504  * (hitp->hit_vpriv[Z] + x + 1.0)
505  - hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y]
506  * (1.0 + 2 * x) >= x * x
507  && hitp->hit_vpriv[Z] >= -1.0
508  && hitp->hit_vpriv[Z] <= 0.0) {
509  hitp->hit_magic = RT_HIT_MAGIC;
510  hitp->hit_dist = k2;
511  hitp->hit_surfno = RHC_NORM_BACK; /* +H */
512  hitp++;
513  }
514  }
515 
516  /* check top plate */
517  if (hitp == &hits[1] && !ZERO(dprime[Z])) {
518  /* 0 or 1 hits so far, this is worthwhile */
519  k1 = -pprime[Z] / dprime[Z]; /* top plate */
520 
521  VJOIN1(hitp->hit_vpriv, pprime, k1, dprime); /* hit' */
522 
523  if (hitp->hit_vpriv[X] >= -1.0 && hitp->hit_vpriv[X] <= 0.0
524  && hitp->hit_vpriv[Y] >= -1.0
525  && hitp->hit_vpriv[Y] <= 1.0) {
526  hitp->hit_magic = RT_HIT_MAGIC;
527  hitp->hit_dist = k1;
528  hitp->hit_surfno = RHC_NORM_TOP; /* -B */
529  hitp++;
530  }
531  }
532 
533  if (hitp != &hits[2]) {
534  return 0; /* MISS */
535  }
536 
537  if (hits[0].hit_dist < hits[1].hit_dist) {
538  /* entry is [0], exit is [1] */
539  struct seg *segp;
540 
541  RT_GET_SEG(segp, ap->a_resource);
542  segp->seg_stp = stp;
543  segp->seg_in = hits[0]; /* struct copy */
544  segp->seg_out = hits[1]; /* struct copy */
545  BU_LIST_INSERT(&(seghead->l), &(segp->l));
546  } else {
547  /* entry is [1], exit is [0] */
548  struct seg *segp;
549 
550  RT_GET_SEG(segp, ap->a_resource);
551  segp->seg_stp = stp;
552  segp->seg_in = hits[1]; /* struct copy */
553  segp->seg_out = hits[0]; /* struct copy */
554  BU_LIST_INSERT(&(seghead->l), &(segp->l));
555  }
556 
557  return 2; /* HIT */
558 }
559 
560 
561 /**
562  * Given ONE ray distance, return the normal and entry/exit point.
563  */
564 void
565 rt_rhc_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
566 {
567  fastf_t c;
568  vect_t can_normal; /* normal to canonical rhc */
569  struct rhc_specific *rhc =
570  (struct rhc_specific *)stp->st_specific;
571 
572  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
573 
574  switch (hitp->hit_surfno) {
575  case RHC_NORM_BODY:
576  c = rhc->rhc_cprime;
577  VSET(can_normal,
578  0.0,
579  hitp->hit_vpriv[Y] * (1.0 + 2.0 * c),
580  -hitp->hit_vpriv[Z] - c - 1.0);
581  MAT4X3VEC(hitp->hit_normal, rhc->rhc_invRoS, can_normal);
582  VUNITIZE(hitp->hit_normal);
583  break;
584 
585  case RHC_NORM_TOP:
586  VREVERSE(hitp->hit_normal, rhc->rhc_Bunit);
587  break;
588 
589  case RHC_NORM_FRT:
590  VREVERSE(hitp->hit_normal, rhc->rhc_Hunit);
591  break;
592 
593  case RHC_NORM_BACK:
594  VMOVE(hitp->hit_normal, rhc->rhc_Hunit);
595  break;
596 
597  default:
598  bu_log("rt_rhc_norm: surfno=%d bad\n", hitp->hit_surfno);
599  break;
600  }
601 }
602 
603 
604 /**
605  * Return the curvature of the rhc.
606  */
607 void
608 rt_rhc_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
609 {
610  fastf_t b, c, rsq, y;
611  fastf_t zp1_sq, zp2; /* 1st deriv sqr, 2nd deriv */
612  struct rhc_specific *rhc =
613  (struct rhc_specific *)stp->st_specific;
614 
615  switch (hitp->hit_surfno) {
616  case RHC_NORM_BODY:
617  /* most nearly flat direction */
618  VMOVE(cvp->crv_pdir, rhc->rhc_Hunit);
619  cvp->crv_c1 = 0;
620  /* k = z'' / (1 + z'^2) ^ 3/2 */
621  b = rhc->rhc_b;
622  c = rhc->rhc_c;
623  y = hitp->hit_point[Y];
624  rsq = rhc->rhc_rsq;
625  zp1_sq = y * (b * b + 2 * b * c) / rsq;
626  zp1_sq *= zp1_sq / (c * c + y * y * (b * b + 2 * b * c) / rsq);
627  zp2 = c * c / (rsq * c * c + y * y * (b * b + 2 * b * c));
628  cvp->crv_c2 = zp2 / pow((1 + zp1_sq), 1.5);
629  break;
630 
631  case RHC_NORM_BACK:
632  case RHC_NORM_FRT:
633  case RHC_NORM_TOP:
634  /* any tangent direction */
635  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
636  cvp->crv_c1 = cvp->crv_c2 = 0;
637  break;
638  }
639 }
640 
641 
642 /**
643  * For a hit on the surface of an rhc, return the (u, v) coordinates
644  * of the hit point, 0 <= u, v <= 1
645  * u = azimuth
646  * v = elevation
647  */
648 void
649 rt_rhc_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
650 {
651  struct rhc_specific *rhc = (struct rhc_specific *)stp->st_specific;
652 
653  vect_t work;
654  vect_t pprime;
655  fastf_t len;
656 
657  if (ap) {
658  RT_CK_APPLICATION(ap);
659  }
660 
661  /*
662  * hit_point is on surface; project back to unit rhc,
663  * creating a vector from vertex to hit point.
664  */
665  VSUB2(work, hitp->hit_point, rhc->rhc_V);
666  MAT4X3VEC(pprime, rhc->rhc_SoR, work);
667 
668  switch (hitp->hit_surfno) {
669  case RHC_NORM_BODY:
670  /* Skin. x, y coordinates define rotation. radius = 1 */
671  len = sqrt(pprime[Y] * pprime[Y] + pprime[Z] * pprime[Z]);
672  uvp->uv_u = acos(pprime[Y] / len) * M_1_PI;
673  uvp->uv_v = -pprime[X]; /* height */
674  break;
675 
676  case RHC_NORM_FRT:
677  case RHC_NORM_BACK:
678  /* end plates - circular mapping, not seamless w/body, top */
679  len = sqrt(pprime[Y] * pprime[Y] + pprime[Z] * pprime[Z]);
680  uvp->uv_u = acos(pprime[Y] / len) * M_1_PI;
681  uvp->uv_v = len; /* rim v = 1 for both plates */
682  break;
683 
684  case RHC_NORM_TOP:
685  uvp->uv_u = 1.0 - (pprime[Y] + 1.0) / 2.0;
686  uvp->uv_v = -pprime[X]; /* height */
687  break;
688  }
689 
690  /* uv_du should be relative to rotation, uv_dv relative to height */
691  uvp->uv_du = uvp->uv_dv = 0;
692 }
693 
694 
695 void
696 rt_rhc_free(struct soltab *stp)
697 {
698  struct rhc_specific *rhc =
699  (struct rhc_specific *)stp->st_specific;
700 
701  BU_PUT(rhc, struct rhc_specific);
702 }
703 
704 
705 /* Our canonical hyperbola in the Y-Z plane has equation
706  * z = +- (a/b) * sqrt(b^2 + y^2), and opens toward +Z and -Z with asymptote
707  * origin at the origin.
708  *
709  * The contour of an rhc in the plane B-R is the positive half of a hyperbola
710  * with asymptote origin at ((|B| + c)Bu), opening toward -B. We can transform
711  * this hyperbola to get an equivalent canonical hyperbola in the Y-Z plane,
712  * opening toward +Z (-B) with asymptote origin at the origin.
713  *
714  * This hyperbola passes through the point (r, |B| + a) (a = c). If we plug the
715  * point (r, |H| + a) into our canonical equation, we can derive b from |B|, a,
716  * and r.
717  *
718  * |B| + a = (a/b) * sqrt(b^2 + r^2)
719  * (|B| + a) / a = b * sqrt(b^2 + r^2)
720  * (|B| + a)^2 / a^2 = 1 + (r^2 / b^2)
721  * ((|B| + a)^2 - a^2) / a^2 = r^2 / b^2
722  * (a^2 * r^2) / ((|B| + a)^2 - a^2) = b^2
723  * (ar) / sqrt((|B| + a)^2 - a^2) = b
724  * (ar) / sqrt(|B| (|B| + 2a)) = b
725  */
726 static fastf_t
727 rhc_hyperbola_b(fastf_t mag_b, fastf_t c, fastf_t r)
728 {
729  return (c * r) / sqrt(mag_b * (mag_b + 2.0 * c));
730 }
731 
732 
733 static fastf_t
734 rhc_hyperbola_y(fastf_t b, fastf_t c, fastf_t z)
735 {
736  return sqrt(b * b * (((z * z) / (c * c)) - 1.0));
737 }
738 
739 
740 /* The contour of an rhc in the plane B-R is the positive half of a hyperbola
741  * with asymptote origin at ((|B| + c)Bu), opening toward -B. We can transform
742  * this hyperbola to get an equivalent hyperbola in the Y-Z plane, opening
743  * toward +Z (-B) with asymptote origin at (0, -(|B| + c)).
744  *
745  * The part of this hyperbola that passes between (0, -(|B| + c)) and (r, 0) is
746  * approximated by num_points points (including (0, -|B|) and (r, 0)).
747  *
748  * The constructed point list is returned (NULL returned on error). Because the
749  * above transformation puts the rhc vertex at the origin and the hyperbola
750  * asymptote origin at (0, -|B| + c), multiplying the z values by -1 gives
751  * corresponding distances along the rhc breadth vector B.
752  */
753 static struct rt_pt_node *
754 rhc_hyperbolic_curve(fastf_t mag_b, fastf_t c, fastf_t r, int num_points)
755 {
756  int count;
757  struct rt_pt_node *curve;
758 
759  if (num_points < 2) {
760  return NULL;
761  }
762 
763  BU_ALLOC(curve, struct rt_pt_node);
764  BU_ALLOC(curve->next, struct rt_pt_node);
765 
766  curve->next->next = NULL;
767  VSET(curve->p, 0, 0, -mag_b);
768  VSET(curve->next->p, 0, r, 0);
769 
770  if (num_points < 3) {
771  return curve;
772  }
773 
774  count = approximate_hyperbolic_curve(curve, c, rhc_hyperbola_b(mag_b, c, r), num_points - 2);
775 
776  if (count != (num_points - 2)) {
777  return NULL;
778  }
779 
780  return curve;
781 }
782 
783 
784 /* plot half of a hyperbolic contour curve using the given (r, b) points (pts),
785  * translation along H (rhc_H), and multiplier for r (rscale)
786  */
787 static void
788 rhc_plot_hyperbolic_curve(
789  struct bu_list *vhead,
790  struct rhc_specific *rhc,
791  struct rt_pt_node *pts,
792  vect_t rhc_H,
793  fastf_t rscale)
794 {
795  vect_t t, Ru, Bu;
796  point_t p;
797  struct rt_pt_node *node;
798 
799  VADD2(t, rhc->rhc_V, rhc_H);
800  VMOVE(Ru, rhc->rhc_Runit);
801  VMOVE(Bu, rhc->rhc_Bunit);
802 
803  VJOIN2(p, t, rscale * pts->p[Y], Ru, -pts->p[Z], Bu);
804  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_MOVE);
805 
806  node = pts->next;
807  while (node != NULL) {
808  VJOIN2(p, t, rscale * node->p[Y], Ru, -node->p[Z], Bu);
809  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_DRAW);
810 
811  node = node->next;
812  }
813 }
814 
815 
816 static void
817 rhc_plot_hyperbolas(
818  struct bu_list *vhead,
819  struct rt_rhc_internal *rhc,
820  struct rt_pt_node *pts)
821 {
822  vect_t rhc_H;
823  struct rhc_specific rhc_s;
824 
825  VMOVE(rhc_s.rhc_V, rhc->rhc_V);
826 
827  VMOVE(rhc_s.rhc_Bunit, rhc->rhc_B);
828  VUNITIZE(rhc_s.rhc_Bunit);
829 
830  VCROSS(rhc_s.rhc_Runit, rhc_s.rhc_Bunit, rhc->rhc_H);
831  VUNITIZE(rhc_s.rhc_Runit);
832 
833  /* plot hyperbolic contour curve of face containing V */
834  VSETALL(rhc_H, 0.0);
835  rhc_plot_hyperbolic_curve(vhead, &rhc_s, pts, rhc_H, 1.0);
836  rhc_plot_hyperbolic_curve(vhead, &rhc_s, pts, rhc_H, -1.0);
837 
838  /* plot hyperbolic contour curve of opposing face */
839  VMOVE(rhc_H, rhc->rhc_H);
840  rhc_plot_hyperbolic_curve(vhead, &rhc_s, pts, rhc_H, 1.0);
841  rhc_plot_hyperbolic_curve(vhead, &rhc_s, pts, rhc_H, -1.0);
842 }
843 
844 
845 static void
846 rhc_plot_curve_connections(
847  struct bu_list *vhead,
848  const struct rt_rhc_internal *rhc,
849  int num_connections)
850 {
851  point_t p;
852  vect_t Yu, Zu;
853  fastf_t b, y, z, z_step;
854  fastf_t mag_Z, zmax;
855  int connections_per_half;
856 
857  if (num_connections < 1) {
858  return;
859  }
860 
861  VMOVE(Zu, rhc->rhc_B);
862  VCROSS(Yu, rhc->rhc_H, Zu);
863  VUNITIZE(Yu);
864  VUNITIZE(Zu);
865 
866  mag_Z = MAGNITUDE(rhc->rhc_B);
867 
868  b = rhc_hyperbola_b(mag_Z, rhc->rhc_c, rhc->rhc_r);
869 
870  connections_per_half = .5 + num_connections / 2.0;
871  z_step = mag_Z / (connections_per_half + 1.0);
872 
873  zmax = mag_Z + rhc->rhc_c;
874  for (z = 0.0; z <= mag_Z; z += z_step) {
875  y = rhc_hyperbola_y(b, rhc->rhc_c, zmax - z);
876 
877  /* connect faces on one side of the curve */
878  VJOIN2(p, rhc->rhc_V, z, Zu, -y, Yu);
879  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_MOVE);
880 
881  VADD2(p, p, rhc->rhc_H);
882  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_DRAW);
883 
884  /* connect the faces on the other side */
885  VJOIN2(p, rhc->rhc_V, z, Zu, y, Yu);
886  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_MOVE);
887 
888  VADD2(p, p, rhc->rhc_H);
889  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_DRAW);
890  }
891 }
892 
893 
894 static int
895 rhc_curve_points(
896  struct rt_rhc_internal *rhc,
897  const struct rt_view_info *info)
898 {
899  fastf_t height, halfwidth, est_curve_length;
900  point_t p0, p1;
901 
902  height = -MAGNITUDE(rhc->rhc_B);
903  halfwidth = rhc->rhc_r;
904 
905  VSET(p0, 0, 0, height);
906  VSET(p1, 0, halfwidth, 0);
907 
908  est_curve_length = 2.0 * DIST_PT_PT(p0, p1);
909 
910  return est_curve_length / info->point_spacing;
911 }
912 
913 
914 int
915 rt_rhc_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
916 {
917  point_t p;
918  vect_t rhc_R;
919  int num_curve_points, num_connections;
920  struct rt_rhc_internal *rhc;
921  struct rt_pt_node *pts, *node, *tmp;
922 
923  BU_CK_LIST_HEAD(info->vhead);
924  RT_CK_DB_INTERNAL(ip);
925 
926  rhc = (struct rt_rhc_internal *)ip->idb_ptr;
927  if (!rhc_is_valid(rhc)) {
928  return -2;
929  }
930 
931  num_curve_points = rhc_curve_points(rhc, info);
932 
933  if (num_curve_points < 3) {
934  num_curve_points = 3;
935  }
936 
937  VCROSS(rhc_R, rhc->rhc_B, rhc->rhc_H);
938  VUNITIZE(rhc_R);
939  VSCALE(rhc_R, rhc_R, rhc->rhc_r);
940 
941  pts = rhc_hyperbolic_curve(MAGNITUDE(rhc->rhc_B), rhc->rhc_c, rhc->rhc_r, num_curve_points);
942  rhc_plot_hyperbolas(info->vhead, rhc, pts);
943 
944  node = pts;
945  while (node != NULL) {
946  tmp = node;
947  node = node->next;
948 
949  bu_free(tmp, "rt_pt_node");
950  }
951 
952  /* connect both halves of the hyperbolic contours of the opposing faces */
953  num_connections = primitive_curve_count(ip, info);
954  if (num_connections < 2) {
955  num_connections = 2;
956  }
957 
958  rhc_plot_curve_connections(info->vhead, rhc, num_connections);
959 
960  /* plot rectangular face */
961  VADD2(p, rhc->rhc_V, rhc_R);
963 
964  VADD2(p, p, rhc->rhc_H);
966 
967  VJOIN1(p, p, -2.0, rhc_R);
969 
970  VJOIN1(p, p, -1.0, rhc->rhc_H);
972 
973  VJOIN1(p, p, 2.0, rhc_R);
975 
976  return 0;
977 }
978 
979 
980 int
981 rt_rhc_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))
982 {
983  int i, n;
984  fastf_t b, c, *back, *front, rh;
985  fastf_t dtol, ntol;
986  vect_t Bu, Hu, Ru;
987  mat_t R;
988  mat_t invR;
989  struct rt_rhc_internal *xip;
990  struct rt_pt_node *old, *pos, *pts;
991 
992  BU_CK_LIST_HEAD(vhead);
993  RT_CK_DB_INTERNAL(ip);
994 
995  xip = (struct rt_rhc_internal *)ip->idb_ptr;
996  if (!rhc_is_valid(xip)) {
997  return -2;
998  }
999 
1000  /* compute |B| |H| */
1001  b = MAGNITUDE(xip->rhc_B); /* breadth */
1002  rh = xip->rhc_r; /* rectangular halfwidth */
1003  c = xip->rhc_c; /* dist to asympt origin */
1004 
1005  /* make unit vectors in B, H, and BxH directions */
1006  VMOVE(Hu, xip->rhc_H);
1007  VUNITIZE(Hu);
1008  VMOVE(Bu, xip->rhc_B);
1009  VUNITIZE(Bu);
1010  VCROSS(Ru, Bu, Hu);
1011 
1012  /* Compute R and Rinv matrices */
1013  MAT_IDN(R);
1014  VREVERSE(&R[0], Hu);
1015  VMOVE(&R[4], Ru);
1016  VREVERSE(&R[8], Bu);
1017  bn_mat_trn(invR, R); /* inv of rot mat is trn */
1018 
1019  if (rh < b) {
1020  dtol = primitive_get_absolute_tolerance(ttol, 2.0 * rh);
1021  } else {
1022  dtol = primitive_get_absolute_tolerance(ttol, 2.0 * b);
1023  }
1024 
1025  /* To ensure normal tolerance, remain below this angle */
1026  if (ttol->norm > 0.0) {
1027  ntol = ttol->norm;
1028  } else {
1029  /* tolerate everything */
1030  ntol = M_PI;
1031  }
1032 
1033  /* initial hyperbola approximation is a single segment */
1034  BU_ALLOC(pts, struct rt_pt_node);
1035  BU_ALLOC(pts->next, struct rt_pt_node);
1036 
1037  pts->next->next = NULL;
1038  VSET(pts->p, 0, -rh, 0);
1039  VSET(pts->next->p, 0, rh, 0);
1040  /* 2 endpoints in 1st approximation */
1041  n = 2;
1042  /* recursively break segment 'til within error tolerances */
1043  n += rt_mk_hyperbola(pts, rh, b, c, dtol, ntol);
1044 
1045  /* get mem for arrays */
1046  front = (fastf_t *)bu_malloc(3 * n * sizeof(fastf_t), "fast_t");
1047  back = (fastf_t *)bu_malloc(3 * n * sizeof(fastf_t), "fast_t");
1048 
1049  /* generate front & back plates in world coordinates */
1050  pos = pts;
1051  i = 0;
1052 
1053  while (pos) {
1054  /* rotate back to original position */
1055  MAT4X3VEC(&front[i], invR, pos->p);
1056  /* move to origin vertex origin */
1057  VADD2(&front[i], &front[i], xip->rhc_V);
1058  /* extrude front to create back plate */
1059  VADD2(&back[i], &front[i], xip->rhc_H);
1060  i += 3;
1061  old = pos;
1062  pos = pos->next;
1063  bu_free ((char *)old, "rt_pt_node");
1064  }
1065 
1066  /* Draw the front */
1067  RT_ADD_VLIST(vhead, &front[(n - 1)*ELEMENTS_PER_VECT],
1069 
1070  for (i = 0; i < n; i++) {
1071  RT_ADD_VLIST(vhead, &front[i * ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1072  }
1073 
1074  /* Draw the back */
1075  RT_ADD_VLIST(vhead, &back[(n - 1)*ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1076 
1077  for (i = 0; i < n; i++) {
1078  RT_ADD_VLIST(vhead, &back[i * ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1079  }
1080 
1081  /* Draw connections */
1082  for (i = 0; i < n; i++) {
1083  RT_ADD_VLIST(vhead, &front[i * ELEMENTS_PER_VECT], BN_VLIST_LINE_MOVE);
1084  RT_ADD_VLIST(vhead, &back[i * ELEMENTS_PER_VECT], BN_VLIST_LINE_DRAW);
1085  }
1086 
1087  /* free mem */
1088  bu_free((char *)front, "fastf_t");
1089  bu_free((char *)back, "fastf_t");
1090 
1091  return 0;
1092 }
1093 
1094 
1095 /*
1096  r: rectangular halfwidth
1097  b: breadth
1098  c: distance to asymptote origin
1099 */
1100 int
1102 {
1103  fastf_t A, B, C, discr, dist, intr, j, k, m, theta0, theta1, z0;
1104  int n;
1105  point_t mpt, p0, p1;
1106  vect_t norm_line, norm_hyperb;
1107  struct rt_pt_node *newpt;
1108 
1109 #define RHC_TOL .0001
1110  /* endpoints of segment approximating hyperbola */
1111  VMOVE(p0, pts->p);
1112  VMOVE(p1, pts->next->p);
1113  /* slope and intercept of segment */
1114  m = (p1[Z] - p0[Z]) / (p1[Y] - p0[Y]);
1115  intr = p0[Z] - m * p0[Y];
1116  /* find point on hyperbola with max dist between hyperbola and line */
1117  j = b + c;
1118  k = 1.0 - m * m * r * r / (b * (b + 2.0 * c));
1119  A = k;
1120  B = 2.0 * j * k;
1121  C = j * j * k - c * c;
1122  discr = sqrt(B * B - 4.0 * A * C);
1123  z0 = (-B + discr) / (2.0 * A);
1124 
1125  if (z0 + RHC_TOL >= -b) {
1126  /* use top sheet of hyperboloid */
1127  mpt[Z] = z0;
1128  } else {
1129  mpt[Z] = (-B - discr) / (2.0 * A);
1130  }
1131 
1132  if (NEAR_ZERO(mpt[Z], RHC_TOL)) {
1133  mpt[Z] = 0.;
1134  }
1135 
1136  mpt[X] = 0;
1137  mpt[Y] = ((mpt[Z] + b + c) * (mpt[Z] + b + c) - c * c) / (b * (b + 2.0 * c));
1138 
1139  if (NEAR_ZERO(mpt[Y], RHC_TOL)) {
1140  mpt[Y] = 0.0;
1141  }
1142 
1143  mpt[Y] = r * sqrt(mpt[Y]);
1144 
1145  if (p0[Y] < 0.0) {
1146  mpt[Y] = -mpt[Y];
1147  }
1148 
1149  /* max distance between that point and line */
1150  dist = fabs(m * mpt[Y] - mpt[Z] + intr) / sqrt(m * m + 1);
1151  /* angles between normal of line and of hyperbola at line endpoints */
1152  VSET(norm_line, m, -1.0, 0.0);
1153  VSET(norm_hyperb, 0., (2.0 * c + 1.0) / (p0[Z] + c + 1.0), -1.0);
1154  VUNITIZE(norm_line);
1155  VUNITIZE(norm_hyperb);
1156  theta0 = fabs(acos(VDOT(norm_line, norm_hyperb)));
1157  VSET(norm_hyperb, 0.0, (2.0 * c + 1.0) / (p1[Z] + c + 1.0), -1.0);
1158  VUNITIZE(norm_hyperb);
1159  theta1 = fabs(acos(VDOT(norm_line, norm_hyperb)));
1160 
1161  /* split segment at widest point if not within error tolerances */
1162  if (dist > dtol || theta0 > ntol || theta1 > ntol) {
1163  /* split segment */
1164  BU_ALLOC(newpt, struct rt_pt_node);
1165  VMOVE(newpt->p, mpt);
1166  newpt->next = pts->next;
1167  pts->next = newpt;
1168  /* keep track of number of pts added */
1169  n = 1;
1170  /* recurse on first new segment */
1171  n += rt_mk_hyperbola(pts, r, b, c, dtol, ntol);
1172  /* recurse on second new segment */
1173  n += rt_mk_hyperbola(newpt, r, b, c, dtol, ntol);
1174  } else {
1175  n = 0;
1176  }
1177 
1178  return n;
1179 }
1180 
1181 
1182 /**
1183  * Returns -
1184  * -1 failure
1185  * 0 OK. *r points to nmgregion that holds this tessellation.
1186  */
1187 int
1188 rt_rhc_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
1189 {
1190  int i, j, n;
1191  fastf_t b, c, *back, *front, rh;
1192  fastf_t dtol, ntol;
1193  vect_t Bu, Hu, Ru;
1194  mat_t R;
1195  mat_t invR;
1196  struct rt_rhc_internal *xip;
1197  struct rt_pt_node *old, *pos, *pts;
1198  struct shell *s;
1199  struct faceuse **outfaceuses;
1200  struct vertex **vfront, **vback, **vtemp, *vertlist[4];
1201  vect_t *norms;
1202  fastf_t bb_plus_2bc, b_plus_c, r_sq;
1203  int failure = 0;
1204 
1205  NMG_CK_MODEL(m);
1206  BN_CK_TOL(tol);
1207  RT_CK_TESS_TOL(ttol);
1208 
1209  RT_CK_DB_INTERNAL(ip);
1210 
1211  xip = (struct rt_rhc_internal *)ip->idb_ptr;
1212  if (!rhc_is_valid(xip)) {
1213  return -2;
1214  }
1215 
1216  /* compute |B| |H| */
1217  b = MAGNITUDE(xip->rhc_B); /* breadth */
1218  rh = xip->rhc_r; /* rectangular halfwidth */
1219  c = xip->rhc_c; /* dist to asympt origin */
1220 
1221  /* make unit vectors in B, H, and BxH directions */
1222  VMOVE(Hu, xip->rhc_H);
1223  VUNITIZE(Hu);
1224  VMOVE(Bu, xip->rhc_B);
1225  VUNITIZE(Bu);
1226  VCROSS(Ru, Bu, Hu);
1227 
1228  /* Compute R and Rinv matrices */
1229  MAT_IDN(R);
1230  VREVERSE(&R[0], Hu);
1231  VMOVE(&R[4], Ru);
1232  VREVERSE(&R[8], Bu);
1233  bn_mat_trn(invR, R); /* inv of rot mat is trn */
1234 
1235  if (rh < b) {
1236  dtol = primitive_get_absolute_tolerance(ttol, 2.0 * rh);
1237  } else {
1238  dtol = primitive_get_absolute_tolerance(ttol, 2.0 * b);
1239  }
1240 
1241  /* To ensure normal tolerance, remain below this angle */
1242  if (ttol->norm > 0.0) {
1243  ntol = ttol->norm;
1244  } else {
1245  /* tolerate everything */
1246  ntol = M_PI;
1247  }
1248 
1249  /* initial hyperbola approximation is a single segment */
1250  BU_ALLOC(pts, struct rt_pt_node);
1251  BU_ALLOC(pts->next, struct rt_pt_node);
1252 
1253  pts->next->next = NULL;
1254  VSET(pts->p, 0, -rh, 0);
1255  VSET(pts->next->p, 0, rh, 0);
1256  /* 2 endpoints in 1st approximation */
1257  n = 2;
1258  /* recursively break segment 'til within error tolerances */
1259  n += rt_mk_hyperbola(pts, rh, b, c, dtol, ntol);
1260 
1261  /* get mem for arrays */
1262  front = (fastf_t *)bu_malloc(3 * n * sizeof(fastf_t), "fastf_t");
1263  back = (fastf_t *)bu_malloc(3 * n * sizeof(fastf_t), "fastf_t");
1264  norms = (vect_t *)bu_calloc(n, sizeof(vect_t), "rt_rhc_tess: norms");
1265  vfront = (struct vertex **)bu_malloc((n + 1) * sizeof(struct vertex *), "vertex *");
1266  vback = (struct vertex **)bu_malloc((n + 1) * sizeof(struct vertex *), "vertex *");
1267  vtemp = (struct vertex **)bu_malloc((n + 1) * sizeof(struct vertex *), "vertex *");
1268  outfaceuses =
1269  (struct faceuse **)bu_malloc((n + 2) * sizeof(struct faceuse *), "faceuse *");
1270 
1271  if (!front || !back || !vfront || !vback || !vtemp || !outfaceuses) {
1272  fprintf(stderr, "rt_rhc_tess: no memory!\n");
1273  goto fail;
1274  }
1275 
1276  /* generate front & back plates in world coordinates */
1277  bb_plus_2bc = b * b + 2.0 * b * c;
1278  b_plus_c = b + c;
1279  r_sq = rh * rh;
1280  pos = pts;
1281  i = 0;
1282  j = 0;
1283 
1284  while (pos) {
1285  vect_t tmp_norm;
1286 
1287  /* calculate normal for 2D hyperbola */
1288  VSET(tmp_norm, 0.0, pos->p[Y]*bb_plus_2bc, (-r_sq * (pos->p[Z] + b_plus_c)));
1289  MAT4X3VEC(norms[j], invR, tmp_norm);
1290  VUNITIZE(norms[j]);
1291  /* rotate back to original position */
1292  MAT4X3VEC(&front[i], invR, pos->p);
1293  /* move to origin vertex origin */
1294  VADD2(&front[i], &front[i], xip->rhc_V);
1295  /* extrude front to create back plate */
1296  VADD2(&back[i], &front[i], xip->rhc_H);
1297  i += 3;
1298  j++;
1299  old = pos;
1300  pos = pos->next;
1301  bu_free ((char *)old, "rt_pt_node");
1302  }
1303 
1304  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
1305  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
1306 
1307  for (i = 0; i < n; i++) {
1308  vfront[i] = vtemp[i] = (struct vertex *)0;
1309  }
1310 
1311  /* Front face topology. Verts are considered to go CCW */
1312  outfaceuses[0] = nmg_cface(s, vfront, n);
1313 
1314  (void)nmg_mark_edges_real(&outfaceuses[0]->l.magic);
1315 
1316  /* Back face topology. Verts must go in opposite dir (CW) */
1317  outfaceuses[1] = nmg_cface(s, vtemp, n);
1318 
1319  for (i = 0; i < n; i++) {
1320  vback[i] = vtemp[n - 1 - i];
1321  }
1322 
1323  (void)nmg_mark_edges_real(&outfaceuses[1]->l.magic);
1324 
1325  /* Duplicate [0] as [n] to handle loop end condition, below */
1326  vfront[n] = vfront[0];
1327  vback[n] = vback[0];
1328 
1329  /* Build topology for all the rectangular side faces (n of them)
1330  * connecting the front and back faces.
1331  * increasing indices go towards counter-clockwise (CCW).
1332  */
1333  for (i = 0; i < n; i++) {
1334  vertlist[0] = vfront[i]; /* from top, */
1335  vertlist[1] = vback[i]; /* straight down, */
1336  vertlist[2] = vback[i + 1]; /* to left, */
1337  vertlist[3] = vfront[i + 1]; /* straight up. */
1338  outfaceuses[2 + i] = nmg_cface(s, vertlist, 4);
1339  }
1340 
1341  (void)nmg_mark_edges_real(&outfaceuses[n + 1]->l.magic);
1342 
1343  for (i = 0; i < n; i++) {
1344  NMG_CK_VERTEX(vfront[i]);
1345  NMG_CK_VERTEX(vback[i]);
1346  }
1347 
1348  /* Associate the vertex geometry, CCW */
1349  for (i = 0; i < n; i++) {
1350  nmg_vertex_gv(vfront[i], &front[3 * (i)]);
1351  }
1352 
1353  for (i = 0; i < n; i++) {
1354  nmg_vertex_gv(vback[i], &back[3 * (i)]);
1355  }
1356 
1357  /* Associate the face geometry */
1358  for (i = 0; i < n + 2; i++) {
1359  if (nmg_fu_planeeqn(outfaceuses[i], tol) < 0) {
1360  failure = (-1);
1361  goto fail;
1362  }
1363  }
1364 
1365  /* Associate vertexuse normals */
1366  for (i = 0; i < n; i++) {
1367  struct vertexuse *vu;
1368  struct faceuse *fu;
1369  vect_t rev_norm;
1370 
1371  VREVERSE(rev_norm, norms[i]);
1372 
1373  /* do "front" vertices */
1374  NMG_CK_VERTEX(vfront[i]);
1375 
1376  for (BU_LIST_FOR(vu, vertexuse, &vfront[i]->vu_hd)) {
1377  NMG_CK_VERTEXUSE(vu);
1378  fu = nmg_find_fu_of_vu(vu);
1379  NMG_CK_FACEUSE(fu);
1380 
1381  if (fu->f_p == outfaceuses[0]->f_p ||
1382  fu->f_p == outfaceuses[1]->f_p ||
1383  fu->f_p == outfaceuses[n + 1]->f_p) {
1384  continue; /* skip flat faces */
1385  }
1386 
1387  if (fu->orientation == OT_SAME) {
1388  nmg_vertexuse_nv(vu, norms[i]);
1389  } else if (fu->orientation == OT_OPPOSITE) {
1390  nmg_vertexuse_nv(vu, rev_norm);
1391  }
1392  }
1393 
1394  /* and "back" vertices */
1395  NMG_CK_VERTEX(vback[i]);
1396 
1397  for (BU_LIST_FOR(vu, vertexuse, &vback[i]->vu_hd)) {
1398  NMG_CK_VERTEXUSE(vu);
1399  fu = nmg_find_fu_of_vu(vu);
1400  NMG_CK_FACEUSE(fu);
1401 
1402  if (fu->f_p == outfaceuses[0]->f_p ||
1403  fu->f_p == outfaceuses[1]->f_p ||
1404  fu->f_p == outfaceuses[n + 1]->f_p) {
1405  continue; /* skip flat faces */
1406  }
1407 
1408  if (fu->orientation == OT_SAME) {
1409  nmg_vertexuse_nv(vu, norms[i]);
1410  } else if (fu->orientation == OT_OPPOSITE) {
1411  nmg_vertexuse_nv(vu, rev_norm);
1412  }
1413  }
1414  }
1415 
1416  /* Glue the edges of different outward pointing face uses together */
1417  nmg_gluefaces(outfaceuses, n + 2, tol);
1418 
1419  /* Compute "geometry" for region and shell */
1420  nmg_region_a(*r, tol);
1421 
1422 fail:
1423  /* free mem */
1424  bu_free((char *)front, "fastf_t");
1425  bu_free((char *)back, "fastf_t");
1426  bu_free((char *)vfront, "vertex *");
1427  bu_free((char *)vback, "vertex *");
1428  bu_free((char *)vtemp, "vertex *");
1429  bu_free((char *)norms, "rt_rhc_tess: norms");
1430  bu_free((char *)outfaceuses, "faceuse *");
1431 
1432  return failure;
1433 }
1434 
1435 
1436 /**
1437  * Import an RHC from the database format to the internal format.
1438  * Apply modeling transformations as well.
1439  */
1440 int
1441 rt_rhc_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
1442 {
1443  struct rt_rhc_internal *xip;
1444  union record *rp;
1445  vect_t v1, v2, v3;
1446 
1447  if (dbip) {
1448  RT_CK_DBI(dbip);
1449  }
1450 
1451  BU_CK_EXTERNAL(ep);
1452  rp = (union record *)ep->ext_buf;
1453 
1454  /* Check record type */
1455  if (rp->u_id != ID_SOLID) {
1456  bu_log("rt_rhc_import4: defective record\n");
1457  return -1;
1458  }
1459 
1460  RT_CK_DB_INTERNAL(ip);
1461  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1462  ip->idb_type = ID_RHC;
1463  ip->idb_meth = &OBJ[ID_RHC];
1464  BU_ALLOC(ip->idb_ptr, struct rt_rhc_internal);
1465 
1466  xip = (struct rt_rhc_internal *)ip->idb_ptr;
1467  xip->rhc_magic = RT_RHC_INTERNAL_MAGIC;
1468 
1469  /* Warning: type conversion */
1470  if (mat == NULL) {
1471  mat = bn_mat_identity;
1472  }
1473 
1474  if (dbip->dbi_version < 0) {
1475  flip_fastf_float(v1, &rp->s.s_values[0 * 3], 1, 1);
1476  flip_fastf_float(v2, &rp->s.s_values[1 * 3], 1, 1);
1477  flip_fastf_float(v3, &rp->s.s_values[2 * 3], 1, 1);
1478  } else {
1479  VMOVE(v1, &rp->s.s_values[0 * 3]);
1480  VMOVE(v2, &rp->s.s_values[1 * 3]);
1481  VMOVE(v3, &rp->s.s_values[2 * 3]);
1482  }
1483 
1484  MAT4X3PNT(xip->rhc_V, mat, v1);
1485  MAT4X3VEC(xip->rhc_H, mat, v2);
1486  MAT4X3VEC(xip->rhc_B, mat, v3);
1487 
1488  if (dbip->dbi_version < 0) {
1489  v1[X] = flip_dbfloat(rp->s.s_values[3 * 3 + 0]);
1490  v1[Y] = flip_dbfloat(rp->s.s_values[3 * 3 + 1]);
1491  } else {
1492  v1[X] = rp->s.s_values[3 * 3 + 0];
1493  v1[Y] = rp->s.s_values[3 * 3 + 1];
1494  }
1495 
1496  xip->rhc_r = v1[X] / mat[15];
1497  xip->rhc_c = v2[Y] / mat[15];
1498 
1499  if (xip->rhc_r <= SMALL_FASTF || xip->rhc_c <= SMALL_FASTF) {
1500  bu_log("rt_rhc_import4: r or c are zero\n");
1501  bu_free((char *)ip->idb_ptr, "rt_rhc_import4: ip->idb_ptr");
1502  return -1;
1503  }
1504 
1505  return 0; /* OK */
1506 }
1507 
1508 
1509 /**
1510  * The name is added by the caller, in the usual place.
1511  */
1512 int
1513 rt_rhc_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1514 {
1515  struct rt_rhc_internal *xip;
1516  union record *rhc;
1517 
1518  if (dbip) {
1519  RT_CK_DBI(dbip);
1520  }
1521 
1522  RT_CK_DB_INTERNAL(ip);
1523 
1524  if (ip->idb_type != ID_RHC) {
1525  return -1;
1526  }
1527 
1528  xip = (struct rt_rhc_internal *)ip->idb_ptr;
1529  if (!rhc_is_valid(xip)) {
1530  return -1;
1531  }
1532 
1533  BU_CK_EXTERNAL(ep);
1534  ep->ext_nbytes = sizeof(union record);
1535  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "rhc external");
1536  rhc = (union record *)ep->ext_buf;
1537 
1538  rhc->s.s_id = ID_SOLID;
1539  rhc->s.s_type = RHC;
1540 
1541  /* Warning: type conversion */
1542  VSCALE(&rhc->s.s_values[0 * 3], xip->rhc_V, local2mm);
1543  VSCALE(&rhc->s.s_values[1 * 3], xip->rhc_H, local2mm);
1544  VSCALE(&rhc->s.s_values[2 * 3], xip->rhc_B, local2mm);
1545  rhc->s.s_values[3 * 3] = xip->rhc_r * local2mm;
1546  rhc->s.s_values[3 * 3 + 1] = xip->rhc_c * local2mm;
1547 
1548  return 0;
1549 }
1550 
1551 
1552 /**
1553  * Import an RHC from the database format to the internal format.
1554  * Apply modeling transformations as well.
1555  */
1556 int
1557 rt_rhc_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
1558 {
1559  struct rt_rhc_internal *xip;
1560 
1561  /* must be double for import and export */
1562  double vec[11];
1563 
1564  if (dbip) {
1565  RT_CK_DBI(dbip);
1566  }
1567 
1568  BU_CK_EXTERNAL(ep);
1570 
1571  RT_CK_DB_INTERNAL(ip);
1572  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1573  ip->idb_type = ID_RHC;
1574  ip->idb_meth = &OBJ[ID_RHC];
1575  BU_ALLOC(ip->idb_ptr, struct rt_rhc_internal);
1576 
1577  xip = (struct rt_rhc_internal *)ip->idb_ptr;
1578  xip->rhc_magic = RT_RHC_INTERNAL_MAGIC;
1579 
1580  /* Convert from database (network) to internal (host) format */
1581  bu_cv_ntohd((unsigned char *)vec, ep->ext_buf, 11);
1582 
1583  /* Apply modeling transformations */
1584  if (mat == NULL) {
1585  mat = bn_mat_identity;
1586  }
1587 
1588  MAT4X3PNT(xip->rhc_V, mat, &vec[0 * 3]);
1589  MAT4X3VEC(xip->rhc_H, mat, &vec[1 * 3]);
1590  MAT4X3VEC(xip->rhc_B, mat, &vec[2 * 3]);
1591  xip->rhc_r = vec[3 * 3] / mat[15];
1592  xip->rhc_c = vec[3 * 3 + 1] / mat[15];
1593 
1594  if (xip->rhc_r <= SMALL_FASTF || xip->rhc_c <= SMALL_FASTF) {
1595  bu_log("rt_rhc_import4: r or c are zero\n");
1596  bu_free((char *)ip->idb_ptr, "rt_rhc_import4: ip->idb_ptr");
1597  return -1;
1598  }
1599 
1600  return 0; /* OK */
1601 }
1602 
1603 
1604 /**
1605  * The name is added by the caller, in the usual place.
1606  */
1607 int
1608 rt_rhc_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1609 {
1610  struct rt_rhc_internal *xip;
1611 
1612  /* must be double for import and export */
1613  double vec[11];
1614 
1615  if (dbip) {
1616  RT_CK_DBI(dbip);
1617  }
1618 
1619  RT_CK_DB_INTERNAL(ip);
1620 
1621  if (ip->idb_type != ID_RHC) {
1622  return -1;
1623  }
1624 
1625  xip = (struct rt_rhc_internal *)ip->idb_ptr;
1626  if (!rhc_is_valid(xip)) {
1627  return -1;
1628  }
1629 
1630  BU_CK_EXTERNAL(ep);
1631  ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * 11;
1632  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "rhc external");
1633 
1634  /* scale 'em into local buffer */
1635  VSCALE(&vec[0 * 3], xip->rhc_V, local2mm);
1636  VSCALE(&vec[1 * 3], xip->rhc_H, local2mm);
1637  VSCALE(&vec[2 * 3], xip->rhc_B, local2mm);
1638  vec[3 * 3] = xip->rhc_r * local2mm;
1639  vec[3 * 3 + 1] = xip->rhc_c * local2mm;
1640 
1641  /* Convert from internal (host) to database (network) format */
1642  bu_cv_htond(ep->ext_buf, (unsigned char *)vec, 11);
1643 
1644  return 0;
1645 }
1646 
1647 
1648 /**
1649  * Make human-readable formatted presentation of this solid.
1650  * First line describes type of solid.
1651  * Additional lines are indented one tab, and give parameter values.
1652  */
1653 int
1654 rt_rhc_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1655 {
1656  struct rt_rhc_internal *xip =
1657  (struct rt_rhc_internal *)ip->idb_ptr;
1658  char buf[256];
1659 
1660  RT_RHC_CK_MAGIC(xip);
1661  bu_vls_strcat(str, "Right Hyperbolic Cylinder (RHC)\n");
1662 
1663  if (!verbose) {
1664  return 0;
1665  }
1666 
1667  sprintf(buf, "\tV (%g, %g, %g)\n",
1668  INTCLAMP(xip->rhc_V[X] * mm2local),
1669  INTCLAMP(xip->rhc_V[Y] * mm2local),
1670  INTCLAMP(xip->rhc_V[Z] * mm2local));
1671  bu_vls_strcat(str, buf);
1672 
1673  sprintf(buf, "\tB (%g, %g, %g) mag=%g\n",
1674  INTCLAMP(xip->rhc_B[X] * mm2local),
1675  INTCLAMP(xip->rhc_B[Y] * mm2local),
1676  INTCLAMP(xip->rhc_B[Z] * mm2local),
1677  INTCLAMP(MAGNITUDE(xip->rhc_B) * mm2local));
1678  bu_vls_strcat(str, buf);
1679 
1680  sprintf(buf, "\tH (%g, %g, %g) mag=%g\n",
1681  INTCLAMP(xip->rhc_H[X] * mm2local),
1682  INTCLAMP(xip->rhc_H[Y] * mm2local),
1683  INTCLAMP(xip->rhc_H[Z] * mm2local),
1684  INTCLAMP(MAGNITUDE(xip->rhc_H) * mm2local));
1685  bu_vls_strcat(str, buf);
1686 
1687  sprintf(buf, "\tr=%g\n", INTCLAMP(xip->rhc_r * mm2local));
1688  bu_vls_strcat(str, buf);
1689 
1690  sprintf(buf, "\tc=%g\n", INTCLAMP(xip->rhc_c * mm2local));
1691  bu_vls_strcat(str, buf);
1692 
1693  return 0;
1694 }
1695 
1696 
1697 /**
1698  * Free the storage associated with the rt_db_internal version of this solid.
1699  */
1700 void
1702 {
1703  struct rt_rhc_internal *xip;
1704 
1705  RT_CK_DB_INTERNAL(ip);
1706 
1707  xip = (struct rt_rhc_internal *)ip->idb_ptr;
1708  RT_RHC_CK_MAGIC(xip);
1709  xip->rhc_magic = 0; /* sanity */
1710 
1711  bu_free((char *)xip, "rhc ifree");
1712  ip->idb_ptr = ((void *)0); /* sanity */
1713 }
1714 
1715 
1716 int
1717 rt_rhc_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
1718 {
1719  if (ip) {
1720  RT_CK_DB_INTERNAL(ip);
1721  }
1722 
1723  return 0; /* OK */
1724 }
1725 
1726 
1727 static int
1728 rhc_is_valid(struct rt_rhc_internal *rhc)
1729 {
1730  fastf_t mag_b, mag_h, cos_angle_bh;
1731 
1732  RT_RHC_CK_MAGIC(rhc);
1733 
1734  mag_b = MAGNITUDE(rhc->rhc_B);
1735  mag_h = MAGNITUDE(rhc->rhc_H);
1736 
1737  /* check for |H| > 0, |B| > 0, |R| > 0, c > 0 */
1738  if (NEAR_ZERO(mag_h, RT_LEN_TOL)
1739  || NEAR_ZERO(mag_b, RT_LEN_TOL)
1740  || NEAR_ZERO(rhc->rhc_r, RT_LEN_TOL)
1741  || NEAR_ZERO(rhc->rhc_c, RT_LEN_TOL))
1742  {
1743  return 0;
1744  }
1745 
1746  /* check B orthogonal to H */
1747  cos_angle_bh = VDOT(rhc->rhc_B, rhc->rhc_H) / (mag_b * mag_h);
1748  if (!NEAR_ZERO(cos_angle_bh, RT_DOT_TOL)) {
1749  return 0;
1750  }
1751 
1752  return 1;
1753 }
1754 
1755 
1756 void
1757 rt_rhc_surf_area(fastf_t *area, const struct rt_db_internal *ip)
1758 {
1759  struct rt_rhc_internal *rip;
1760  fastf_t A, arclen, integralArea, a, b, magB, sqrt_ra, height;
1761 
1762  fastf_t h;
1763  fastf_t sumodds = 0, sumevens = 0, x = 0;
1764  int i, j;
1765 
1766  /**
1767  * n is the number of divisions to use when using Simpson's
1768  * composite rule below to approximate the integral.
1769  *
1770  * A value of n = 1000000 should be enough to ensure that the
1771  * approximation is accurate to at least 10 decimal places. The
1772  * accuracy of the approximation increases by about 2 d.p with
1773  * each added 0 onto the end of the number (i.e. multiply by 10),
1774  * so there is a compromise between accuracy and performance,
1775  * although performance might only be an issue on old slow
1776  * hardware.
1777  *
1778  * I wouldn't recommend setting this less than about 10000,
1779  * because this might cause accuracy to be unsuitable for
1780  * professional or mission critical use.
1781  */
1782  int n = 1000000;
1783 
1784  if (area == NULL || ip == NULL) {
1785  return;
1786  }
1787 
1788  RT_CK_DB_INTERNAL(ip);
1789  rip = (struct rt_rhc_internal *)ip->idb_ptr;
1790  RT_RHC_CK_MAGIC(rip);
1791 
1792  b = rip->rhc_c;
1793  magB = MAGNITUDE(rip->rhc_B);
1794  height = MAGNITUDE(rip->rhc_H);
1795  a = (rip->rhc_r * b) / sqrt(magB * (2.0 * rip->rhc_c + magB));
1796  sqrt_ra = sqrt(rip->rhc_r * rip->rhc_r + b * b);
1797  integralArea = (b / a) * ((2.0 * rip->rhc_r * sqrt_ra) / 2.0 + ((a * a) / 2.0) * (log(sqrt_ra + rip->rhc_r) - log(sqrt_ra - rip->rhc_r)));
1798  A = 2.0 * rip->rhc_r * (rip->rhc_c + magB) - integralArea;
1799 
1800  h = (2.0 * rip->rhc_r) / n;
1801  for (i = 1; i <= (n / 2) - 1; i++) {
1802  x = -rip->rhc_r + 2.0 * i * h;
1803  sumodds += sqrt((b * b * x * x) / (a * a * x * x + pow(a, 4.0)) + 1.0);
1804  }
1805  for (j = 1; j <= (n / 2); j++) {
1806  x = -rip->rhc_r + (2.0 * j - 1.0) * h;
1807  sumevens += sqrt((b * b * x * x) / (a * a * x * x + pow(a, 4.0)) + 1.0);
1808  }
1809  arclen = (h / 3.0) * (sqrt((b * b * rip->rhc_r * rip->rhc_r) / (a * a * rip->rhc_r * rip->rhc_r + pow(a, 4.0)) + 1.0) + 2.0 * sumodds + 4.0 * sumevens + sqrt((b * b * rip->rhc_r * rip->rhc_r) / (a * a * rip->rhc_r * rip->rhc_r + pow(a, 4.0)) + 1.0));
1810 
1811  *area = 2.0 * A + 2.0 * rip->rhc_r * height + arclen * height;
1812 }
1813 
1814 /**
1815  * Computes centroid of a right hyperbolic cylinder
1816  */
1817 void
1818 rt_rhc_centroid(point_t *cent, const struct rt_db_internal *ip)
1819 {
1820  if (cent != NULL && ip != NULL) {
1821  struct rt_rhc_internal *rip;
1822  fastf_t totalArea, guessArea, a, b, magB, sqrt_xa, sqrt_ga, xf, epsilon, high, low, scale_factor;
1823  fastf_t guess = 0.0;
1824  vect_t shift_h;
1825 
1826  RT_CK_DB_INTERNAL(ip);
1827  rip = (struct rt_rhc_internal *)ip->idb_ptr;
1828  RT_RHC_CK_MAGIC(rip);
1829 
1830  magB = MAGNITUDE(rip->rhc_B);
1831  b = rip->rhc_c;
1832  a = (rip->rhc_r * b) / sqrt(magB * (2 * rip->rhc_c + magB));
1833  xf = magB + a;
1834 
1835  /* epsilon is an upperbound on the error */
1836 
1837  epsilon = 0.0001;
1838 
1839  sqrt_xa = sqrt((xf * xf) - (a * a));
1840  totalArea = (b / a) * ((xf * sqrt_xa) - ((a * a) * log(sqrt_xa + xf)) - ((a * a) * log(xf)));
1841 
1842  low = a;
1843  high = xf;
1844 
1845  while (abs(high - low) > epsilon) {
1846  guess = (high + low) / 2.0;
1847  sqrt_ga = sqrt((guess * guess) - (a * a));
1848  guessArea = (b / a) * ((guess * sqrt_ga) - ((a * a) * log(sqrt_ga + guess)) - ((a * a) * log(guess)));
1849 
1850  if (guessArea > totalArea / 2.0) {
1851  high = guess;
1852  } else {
1853  low = guess;
1854  }
1855  }
1856 
1857  scale_factor = 1.0 - ((guess - a) / magB);
1858 
1859  VSCALE(shift_h, rip->rhc_H, 0.5);
1860  VSCALE(*cent, rip->rhc_B, scale_factor);
1861  VADD2(*cent, shift_h, *cent);
1862  }
1863 }
1864 
1865 
1866 /*
1867  * Local Variables:
1868  * mode: C
1869  * tab-width: 8
1870  * indent-tabs-mode: t
1871  * c-file-style: "stroustrup"
1872  * End:
1873  * ex: shiftwidth=4 tabstop=8
1874  */
int rt_rhc_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: rhc.c:1513
fastf_t rhc_b
Definition: rhc.c:188
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
Definition: raytrace.h:800
ustring back
#define RT_LEN_TOL
Definition: raytrace.h:169
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
struct faceuse * nmg_cface(struct shell *s, struct vertex **verts, int n)
Definition: nmg_mod.c:1130
struct hit seg_in
IN information.
Definition: raytrace.h:370
int approximate_hyperbolic_curve(struct rt_pt_node *pts, fastf_t a, fastf_t b, int num_new_points)
fastf_t C[2 *MAX_CNT+1][2 *MAX_CNT+1]
Definition: dsp_brep.cpp:38
const struct bu_structparse rt_rhc_parse[]
Definition: rhc.c:195
Definition: list.h:118
int nmg_fu_planeeqn(struct faceuse *fu, const struct bn_tol *tol)
Definition: nmg_mod.c:1311
int rt_rhc_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: rhc.c:1188
#define RT_DOT_TOL
Definition: raytrace.h:170
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
#define RT_CK_RTI(_p)
Definition: raytrace.h:1833
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
struct soltab * seg_stp
pointer back to soltab
Definition: raytrace.h:372
if lu s
Definition: nmg_mod.c:3860
int rt_rhc_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: rhc.c:209
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
Definition: raytrace.h:368
#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 RHC_NORM_TOP
Definition: rhc.c:370
#define SMALL_FASTF
Definition: defines.h:342
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
int rt_rhc_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: rhc.c:1608
Header file for the BRL-CAD common definitions.
void nmg_vertexuse_nv(struct vertexuse *vu, const fastf_t *norm)
Definition: nmg_mk.c:1719
void flip_fastf_float(fastf_t *ff, const dbfloat_t *fp, int n, int flip)
Definition: db_flip.c:74
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)
int rt_rhc_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
Definition: rhc.c:915
fastf_t rhc_rsq
Definition: rhc.c:191
struct bu_list l
Definition: raytrace.h:369
int rt_mk_hyperbola(struct rt_pt_node *pts, fastf_t r, fastf_t b, fastf_t c, fastf_t dtol, fastf_t ntol)
Definition: rhc.c:1101
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
int rt_rhc_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
Definition: rhc.c:384
if(share_geom)
Definition: nmg_mod.c:3829
int idb_major_type
Definition: raytrace.h:192
mat_t rhc_invRoS
Definition: rhc.c:187
Definition: color.c:49
fastf_t rhc_c
Definition: rhc.c:189
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 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
int rt_rhc_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: rhc.c:1441
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
int rt_rhc_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: rhc.c:273
void rt_rhc_print(const struct soltab *stp)
Definition: rhc.c:354
void rt_rhc_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
Definition: rhc.c:608
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
fastf_t uv_dv
delta in v
Definition: raytrace.h:344
fastf_t primitive_get_absolute_tolerance(const struct rt_tess_tol *ttol, fastf_t rel_to_abs)
void rt_rhc_free(struct soltab *stp)
Definition: rhc.c:696
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
#define RT_CK_TESS_TOL(_p)
Definition: raytrace.h:184
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
struct hit seg_out
OUT information.
Definition: raytrace.h:371
vect_t rhc_Runit
Definition: rhc.c:185
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
#define RHC_NORM_BODY
Definition: rhc.c:369
fastf_t flip_dbfloat(dbfloat_t d)
Definition: db_flip.c:58
void bn_mat_trn(mat_t om, const mat_t im)
#define UNUSED(parameter)
Definition: common.h:239
int nmg_mark_edges_real(const uint32_t *magic_p)
Definition: nmg_misc.c:846
void rt_rhc_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
Definition: rhc.c:565
#define BU_PUT(_ptr, _type)
Definition: malloc.h:215
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
#define BN_CK_TOL(_p)
Definition: tol.h:82
int rt_rhc_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: rhc.c:1654
#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
#define RHC_NORM_BACK
Definition: rhc.c:372
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
int rt_rhc_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: rhc.c:1717
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
#define ID_RHC
Right Hyperbolic Cylinder.
Definition: raytrace.h:476
struct rt_pt_node * next
ptr to next pt
Definition: raytrace.h:1912
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
#define RT_GET_SEG(p, res)
Definition: raytrace.h:379
fastf_t point_spacing
Definition: raytrace.h:1934
#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
mat_t rhc_SoR
Definition: rhc.c:186
int rt_rhc_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: rhc.c:981
void bu_cv_ntohd(unsigned char *out, const unsigned char *in, size_t count)
const struct rt_functab OBJ[]
Definition: table.c:159
void bn_vec_ortho(vect_t out, const vect_t in)
#define R
Definition: msr.c:55
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t uv_du
delta in u
Definition: raytrace.h:343
point_t p
a point
Definition: raytrace.h:1911
void rt_rhc_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: rhc.c:1757
void rt_rhc_centroid(point_t *cent, const struct rt_db_internal *ip)
Definition: rhc.c:1818
point_t rhc_V
Definition: rhc.c:182
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
ustring front
int rt_rhc_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: rhc.c:1557
double norm
normal tol
Definition: raytrace.h:182
int hit_surfno
solid-specific surface indicator
Definition: raytrace.h:255
#define RHC_NORM_FRT
Definition: rhc.c:371
vect_t rhc_Bunit
Definition: rhc.c:183
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
void nmg_gluefaces(struct faceuse **fulist, int n, const struct bn_tol *tol)
Definition: nmg_mod.c:1408
const struct rt_functab * st_meth
pointer to per-solid methods
Definition: raytrace.h:428
fastf_t rhc_cprime
Definition: rhc.c:190
void rt_rhc_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
Definition: rhc.c:649
#define RT_PCOEF_TOL
Definition: raytrace.h:171
void rt_rhc_ifree(struct rt_db_internal *ip)
Definition: rhc.c:1701
int st_id
Solid ident.
Definition: raytrace.h:431
size_t ext_nbytes
Definition: parse.h:210
vect_t rhc_Hunit
Definition: rhc.c:184
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
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
#define RHC_TOL
#define RT_RHC_INTERNAL_MAGIC
Definition: magic.h:106
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 nmg_region_a(struct nmgregion *r, const struct bn_tol *tol)
Definition: nmg_mk.c:2557