BRL-CAD
ehy.c
Go to the documentation of this file.
1 /* E H Y . 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/ehy/ehy.c
23  *
24  * Intersect a ray with an Elliptical Hyperboloid.
25  *
26  * Algorithm -
27  *
28  * Given V, H, R, and B, there is a set of points on this ehy
29  *
30  * { (x, y, z) | (x, y, z) is on ehy }
31  *
32  * Through a series of Affine Transformations, this set of points will
33  * be transformed into a set of points on an ehy located at the origin
34  * with a semi-major axis R1 along the +Y axis, a semi-minor axis R2
35  * along the -X axis, a height H along the -Z axis, a vertex V at the
36  * origin, and a distance c between the tip of the hyperboloid and
37  * vertex of the asymptotic cone.
38  *
39  *
40  * { (x', y', z') | (x', y', z') is on ehy 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) = (R2/(-|R2|))
47  * (R1/(|R1|)) . X
48  * (H /(-|H |))
49  *
50  * and S(X) = (1/|R2| 0 0)
51  * (0 1/|R1| 0) . X
52  * (0 0 1/|H |)
53  *
54  * To find the intersection of a line with the surface of the ehy,
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 ehy. Let
60  * W' be the point of intersection between L' and the unit ehy.
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 ehy; 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 ehy at W' when
73  *
74  * A * k**2 + B * k + C = 0
75  *
76  * where
77  *
78  * A = Dz'**2 - (2*c' + 1)*(Dx'**2 + Dy'**2)
79  * B = 2*((Pz' + c' + 1)*Dz' - (2*c' + 1)*(Dx'*Px' + Dy'Py'))
80  * C = Pz'**2 - (2*c' + 1)*(Px'**2 + Py'**2 - 1) + 2*(c' + 1)*Pz'
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.0 * 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 ehy IFF
104  * -1 <= Wz' <= 0.
105  *
106  * NORMALS. Given the point W on the surface of the ehy, what is the
107  * vector normal to the tangent plane at that point?
108  *
109  * Map W onto the unit ehy, i.e.: W' = S(R(W - V)).
110  *
111  * Plane on unit ehy at W' has a normal vector N' where
112  *
113  * N' = <Wx', Wy', -(z + c + 1) / (2*c + 1)>
114  *
115  * The plane transforms back to the tangent plane at W, and this new
116  * plane (on the original ehy) 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 Q
121  * to Q', then inverse[ transpose(M) ] (H) is perpendicular to Q'.
122  * Here, H and Q are in "prime space" with the unit sphere. [Somehow,
123  * the notation here is backwards]. So, the mapping matrix M =
124  * inverse(S o R), because S o R maps from normal space to the unit
125  * 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), and S =
134  * transpose(S).
135  *
136  * Note that the normal vector produced above will not have unit
137  * length.
138  *
139  * THE TOP PLATE.
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 with
143  * k = (0 - Pz')/Dz'. The solution is within the top plate IFF Wx'**2
144  * + Wy'**2 <= 1.
145  *
146  * The normal for a hit on the top plate is -Hunit.
147  *
148  */
149 
150 #include "common.h"
151 
152 #include <stddef.h>
153 #include <string.h>
154 #include <math.h>
155 #include "bio.h"
156 
157 #include "bu/cv.h"
158 #include "vmath.h"
159 #include "db.h"
160 #include "nmg.h"
161 #include "rtgeom.h"
162 #include "raytrace.h"
163 
164 #include "../../librt_private.h"
165 
166 
167 struct ehy_specific {
168  point_t ehy_V; /* vector to ehy origin */
169  vect_t ehy_Hunit; /* unit H vector */
170  vect_t ehy_Aunit; /* unit vector along semi-major axis */
171  vect_t ehy_Bunit; /* unit vector, A x H, semi-minor axis */
172  mat_t ehy_SoR; /* Scale(Rot(vect)) */
173  mat_t ehy_invRoS; /* invRot(Scale(vect)) */
174  fastf_t ehy_cprime; /* c / |H| */
175 };
176 
177 
178 const struct bu_structparse rt_ehy_parse[] = {
179  { "%f", 3, "V", bu_offsetofarray(struct rt_ehy_internal, ehy_V, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
180  { "%f", 3, "H", bu_offsetofarray(struct rt_ehy_internal, ehy_H, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
181  { "%f", 3, "A", bu_offsetofarray(struct rt_ehy_internal, ehy_Au, fastf_t, X), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
182  { "%f", 1, "r_1", bu_offsetof(struct rt_ehy_internal, ehy_r1), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
183  { "%f", 1, "r_2", bu_offsetof(struct rt_ehy_internal, ehy_r2), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
184  { "%f", 1, "c", bu_offsetof(struct rt_ehy_internal, ehy_c), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
185  { {'\0', '\0', '\0', '\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
186 };
187 
188 
189 static int ehy_is_valid(struct rt_ehy_internal *ehy);
190 
191 /**
192  * Create a bounding RPP for an ehy
193  */
194 int
195 rt_ehy_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
196  struct rt_ehy_internal *xip;
197  vect_t ehy_A, ehy_B, ehy_An, ehy_Bn, ehy_H;
198  vect_t pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8;
199  RT_CK_DB_INTERNAL(ip);
200  xip = (struct rt_ehy_internal *)ip->idb_ptr;
201  RT_EHY_CK_MAGIC(xip);
202 
203  VMOVE(ehy_H, xip->ehy_H);
204  VUNITIZE(ehy_H);
205  VMOVE(ehy_A, xip->ehy_Au);
206  VCROSS(ehy_B, ehy_A, ehy_H);
207 
208  VSETALL((*min), INFINITY);
209  VSETALL((*max), -INFINITY);
210 
211  VSCALE(ehy_A, ehy_A, xip->ehy_r1);
212  VSCALE(ehy_B, ehy_B, xip->ehy_r2);
213  VREVERSE(ehy_An, ehy_A);
214  VREVERSE(ehy_Bn, ehy_B);
215 
216  VADD3(pt1, xip->ehy_V, ehy_A, ehy_B);
217  VADD3(pt2, xip->ehy_V, ehy_A, ehy_Bn);
218  VADD3(pt3, xip->ehy_V, ehy_An, ehy_B);
219  VADD3(pt4, xip->ehy_V, ehy_An, ehy_Bn);
220  VADD4(pt5, xip->ehy_V, ehy_A, ehy_B, xip->ehy_H);
221  VADD4(pt6, xip->ehy_V, ehy_A, ehy_Bn, xip->ehy_H);
222  VADD4(pt7, xip->ehy_V, ehy_An, ehy_B, xip->ehy_H);
223  VADD4(pt8, xip->ehy_V, ehy_An, ehy_Bn, xip->ehy_H);
224 
225  /* Find the RPP of the rotated axis-aligned ehy bbox - that is,
226  * the bounding box the given ehy would have if its height
227  * vector were in the positive Z direction. This does not give
228  * us an optimal bbox except in the case where the ehy is
229  * actually axis aligned to start with, but it's usually
230  * at least a bit better than the bounding sphere RPP. */
231  VMINMAX((*min), (*max), pt1);
232  VMINMAX((*min), (*max), pt2);
233  VMINMAX((*min), (*max), pt3);
234  VMINMAX((*min), (*max), pt4);
235  VMINMAX((*min), (*max), pt5);
236  VMINMAX((*min), (*max), pt6);
237  VMINMAX((*min), (*max), pt7);
238  VMINMAX((*min), (*max), pt8);
239  return 0;
240 }
241 
242 
243 /**
244  * Given a pointer to a GED database record, and a transformation
245  * matrix, determine if this is a valid EHY, and if so, precompute
246  * various terms of the formula.
247  *
248  * Returns -
249  * 0 EHY is OK
250  * !0 Error in description
251  *
252  * Implicit return -
253  * A struct ehy_specific is created, and its address is stored in
254  * stp->st_specific for use by ehy_shot().
255  */
256 int
257 rt_ehy_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
258 {
259  struct rt_ehy_internal *xip;
260  struct ehy_specific *ehy;
261 
262  fastf_t magsq_h;
263  fastf_t mag_h;
264  fastf_t c, r1, r2;
265  mat_t R;
266  mat_t Rinv;
267  mat_t S;
268 
269  RT_CK_DB_INTERNAL(ip);
270 
271  xip = (struct rt_ehy_internal *)ip->idb_ptr;
272 
273  if (!ehy_is_valid(xip)) {
274  return -2;
275  }
276 
277  magsq_h = MAGSQ(xip->ehy_H);
278  mag_h = sqrt(magsq_h);
279  r1 = xip->ehy_r1;
280  r2 = xip->ehy_r2;
281  c = xip->ehy_c;
282 
283  /*
284  * EHY is ok
285  */
286  stp->st_id = ID_EHY; /* set soltab ID */
287  stp->st_meth = &OBJ[ID_EHY];
288 
289  BU_GET(ehy, struct ehy_specific);
290  stp->st_specific = (void *)ehy;
291 
292  /* make unit vectors in A, H, and BxH directions */
293  VMOVE(ehy->ehy_Hunit, xip->ehy_H);
294  VUNITIZE(ehy->ehy_Hunit);
295  VMOVE(ehy->ehy_Aunit, xip->ehy_Au);
296  VCROSS(ehy->ehy_Bunit, ehy->ehy_Aunit, ehy->ehy_Hunit);
297 
298  VMOVE(ehy->ehy_V, xip->ehy_V);
299  ehy->ehy_cprime = c / mag_h;
300 
301  /* Compute R and Rinv matrices */
302  MAT_IDN(R);
303  VREVERSE(&R[0], ehy->ehy_Bunit);
304  VMOVE(&R[4], ehy->ehy_Aunit);
305  VREVERSE(&R[8], ehy->ehy_Hunit);
306  bn_mat_trn(Rinv, R); /* inv of rot mat is trn */
307 
308  /* Compute S */
309  MAT_IDN(S);
310  S[ 0] = 1.0/r2;
311  S[ 5] = 1.0/r1;
312  S[10] = 1.0/mag_h;
313 
314  /* Compute SoR and invRoS */
315  bn_mat_mul(ehy->ehy_SoR, S, R);
316  bn_mat_mul(ehy->ehy_invRoS, Rinv, S);
317 
318  /* Compute bounding sphere */
319  /* bounding sphere center */
320  VJOIN1(stp->st_center, ehy->ehy_V, mag_h / 2.0, ehy->ehy_Hunit);
321  /* bounding radius */
322  stp->st_bradius = sqrt(0.25*magsq_h + r2*r2 + r1*r1);
323  /* approximate bounding radius */
324  stp->st_aradius = stp->st_bradius;
325 
326  if (rt_ehy_bbox(ip, &(stp->st_min), &(stp->st_max), &rtip->rti_tol)) return 1;
327 
328  return 0; /* OK */
329 }
330 
331 
332 void
333 rt_ehy_print(const struct soltab *stp)
334 {
335  const struct ehy_specific *ehy =
336  (struct ehy_specific *)stp->st_specific;
337 
338  VPRINT("V", ehy->ehy_V);
339  VPRINT("Hunit", ehy->ehy_Hunit);
340  VPRINT("Aunit", ehy->ehy_Aunit);
341  VPRINT("Bunit", ehy->ehy_Bunit);
342  fprintf(stderr, "c = %g\n", ehy->ehy_cprime);
343  bn_mat_print("S o R", ehy->ehy_SoR);
344  bn_mat_print("invR o S", ehy->ehy_invRoS);
345 }
346 
347 
348 /* hit_surfno is set to one of these */
349 #define EHY_NORM_BODY (1) /* compute normal */
350 #define EHY_NORM_TOP (2) /* copy ehy_N */
351 
352 
353 /**
354  * Intersect a ray with a ehy. If an intersection occurs, a struct
355  * seg will be acquired and filled in.
356  *
357  * Returns -
358  * 0 MISS
359  * >0 HIT
360  */
361 int
362 rt_ehy_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
363 {
364  struct ehy_specific *ehy =
365  (struct ehy_specific *)stp->st_specific;
366  vect_t dp; /* D' */
367  vect_t pp; /* P' */
368  fastf_t k1, k2; /* distance constants of solution */
369  fastf_t cp; /* c' */
370  vect_t xlated; /* translated vector */
371  struct hit hits[3]; /* 2 potential hit points */
372  struct hit *hitp; /* pointer to hit point */
373 
374  /* for finding roots */
375  fastf_t a, b, c; /* coeffs of polynomial */
376  fastf_t disc; /* discriminant */
377 
378  hitp = &hits[0];
379 
380  /* out, Mat, vect */
381  MAT4X3VEC(dp, ehy->ehy_SoR, rp->r_dir);
382  VSUB2(xlated, rp->r_pt, ehy->ehy_V);
383  MAT4X3VEC(pp, ehy->ehy_SoR, xlated);
384 
385  cp = ehy->ehy_cprime;
386 
387  /* Find roots of the equation, using formula for quadratic */
388 
389  a = dp[Z] * dp[Z]
390  - (2 * cp + 1) * (dp[X] * dp[X] + dp[Y] * dp[Y]);
391  b = 2.0 * (dp[Z] * (pp[Z] + cp + 1)
392  - (2 * cp + 1) * (dp[X] * pp[X] + dp[Y] * pp[Y]));
393  c = pp[Z] * pp[Z]
394  - (2 * cp + 1) * (pp[X] * pp[X] + pp[Y] * pp[Y] - 1.0)
395  + 2 * (cp + 1) * pp[Z];
396  if (!NEAR_ZERO(a, RT_PCOEF_TOL)) {
397  disc = b*b - 4 * a * c;
398  if (disc <= 0)
399  goto check_plates;
400  disc = sqrt(disc);
401 
402  k1 = (-b + disc) / (2.0 * a);
403  k2 = (-b - disc) / (2.0 * a);
404 
405  /*
406  * k1 and k2 are potential solutions to intersection with
407  * side. See if they fall in range.
408  */
409  VJOIN1(hitp->hit_vpriv, pp, k1, dp); /* hit' */
410  if (hitp->hit_vpriv[Z] >= -1.0
411  && hitp->hit_vpriv[Z] <= 0.0) {
412  hitp->hit_magic = RT_HIT_MAGIC;
413  hitp->hit_dist = k1;
414  hitp->hit_surfno = EHY_NORM_BODY; /* compute N */
415  hitp++;
416  }
417 
418  VJOIN1(hitp->hit_vpriv, pp, k2, dp); /* hit' */
419  if (hitp->hit_vpriv[Z] >= -1.0
420  && hitp->hit_vpriv[Z] <= 0.0) {
421  hitp->hit_magic = RT_HIT_MAGIC;
422  hitp->hit_dist = k2;
423  hitp->hit_surfno = EHY_NORM_BODY; /* compute N */
424  hitp++;
425  }
426  } else if (!NEAR_ZERO(b, RT_PCOEF_TOL)) {
427  k1 = -c/b;
428  VJOIN1(hitp->hit_vpriv, pp, k1, dp); /* hit' */
429  if (hitp->hit_vpriv[Z] >= -1.0
430  && hitp->hit_vpriv[Z] <= 0.0) {
431  hitp->hit_magic = RT_HIT_MAGIC;
432  hitp->hit_dist = k1;
433  hitp->hit_surfno = EHY_NORM_BODY; /* compute N */
434  hitp++;
435  }
436  }
437 
438 
439  /*
440  * Check for hitting the top plate.
441  */
442 check_plates:
443  /* check top plate */
444  if (hitp == &hits[1] && !ZERO(dp[Z])) {
445  /* 1 hit so far, this is worthwhile */
446  k1 = -pp[Z] / dp[Z]; /* top plate */
447 
448  VJOIN1(hitp->hit_vpriv, pp, k1, dp); /* hit' */
449  if (hitp->hit_vpriv[X] * hitp->hit_vpriv[X] +
450  hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y] <= 1.0) {
451  hitp->hit_magic = RT_HIT_MAGIC;
452  hitp->hit_dist = k1;
453  hitp->hit_surfno = EHY_NORM_TOP; /* -H */
454  hitp++;
455  }
456  }
457 
458  if (hitp != &hits[2])
459  return 0; /* MISS */
460 
461  if (hits[0].hit_dist < hits[1].hit_dist) {
462  /* entry is [0], exit is [1] */
463  struct seg *segp;
464 
465  RT_GET_SEG(segp, ap->a_resource);
466  segp->seg_stp = stp;
467  segp->seg_in = hits[0]; /* struct copy */
468  segp->seg_out = hits[1]; /* struct copy */
469  BU_LIST_INSERT(&(seghead->l), &(segp->l));
470  } else {
471  /* entry is [1], exit is [0] */
472  struct seg *segp;
473 
474  RT_GET_SEG(segp, ap->a_resource);
475  segp->seg_stp = stp;
476  segp->seg_in = hits[1]; /* struct copy */
477  segp->seg_out = hits[0]; /* struct copy */
478  BU_LIST_INSERT(&(seghead->l), &(segp->l));
479  }
480  return 2; /* HIT */
481 }
482 
483 
484 /**
485  * Given ONE ray distance, return the normal and entry/exit point.
486  */
487 void
488 rt_ehy_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
489 {
490  vect_t can_normal; /* normal to canonical ehy */
491  fastf_t cp, scale;
492  struct ehy_specific *ehy =
493  (struct ehy_specific *)stp->st_specific;
494 
495  VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
496  switch (hitp->hit_surfno) {
497  case EHY_NORM_BODY:
498  cp = ehy->ehy_cprime;
499  VSET(can_normal,
500  hitp->hit_vpriv[X] * (2 * cp + 1),
501  hitp->hit_vpriv[Y] * (2 * cp + 1),
502  -(hitp->hit_vpriv[Z] + cp + 1));
503  MAT4X3VEC(hitp->hit_normal, ehy->ehy_invRoS, can_normal);
504  scale = 1.0 / MAGNITUDE(hitp->hit_normal);
505  VSCALE(hitp->hit_normal, hitp->hit_normal, scale);
506 
507  /* tuck away this scale for the curvature routine */
508  hitp->hit_vpriv[X] = scale;
509  break;
510  case EHY_NORM_TOP:
511  VREVERSE(hitp->hit_normal, ehy->ehy_Hunit);
512  break;
513  default:
514  bu_log("rt_ehy_norm: surfno=%d bad\n", hitp->hit_surfno);
515  break;
516  }
517 }
518 
519 
520 /**
521  * Return the curvature of the ehy.
522  */
523 void
524 rt_ehy_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
525 {
526  fastf_t a, b, c, scale;
527  mat_t M1, M2;
528  struct ehy_specific *ehy =
529  (struct ehy_specific *)stp->st_specific;
530  vect_t u, v; /* basis vectors (with normal) */
531  vect_t vec1, vec2; /* eigen vectors */
532  vect_t tmp;
533 
534  switch (hitp->hit_surfno) {
535  case EHY_NORM_BODY:
536  /*
537  * choose a tangent plane coordinate system (u, v, normal)
538  * form a right-handed triple
539  */
540  bn_vec_ortho(u, hitp->hit_normal);
541  VCROSS(v, hitp->hit_normal, u);
542 
543  /* get the saved away scale factor */
544  scale = - hitp->hit_vpriv[X];
545 
546  MAT_IDN(M1);
547  M1[0] = M1[5] = (4 * ehy->ehy_cprime + 2)
548  / (ehy->ehy_cprime * ehy->ehy_cprime);
549  M1[10] = -1.;
550  /* M1 = invR * S * M1 * S * R */
551  bn_mat_mul(M2, ehy->ehy_invRoS, M1);
552  bn_mat_mul(M1, M2, ehy->ehy_SoR);
553 
554  /* find the second fundamental form */
555  MAT4X3VEC(tmp, ehy->ehy_invRoS, u);
556  a = VDOT(u, tmp) * scale;
557  b = VDOT(v, tmp) * scale;
558  MAT4X3VEC(tmp, ehy->ehy_invRoS, v);
559  c = VDOT(v, tmp) * scale;
560 
561  bn_eigen2x2(&cvp->crv_c1, &cvp->crv_c2, vec1, vec2, a, b, c);
562  VCOMB2(cvp->crv_pdir, vec1[X], u, vec1[Y], v);
563  VUNITIZE(cvp->crv_pdir);
564  break;
565  case EHY_NORM_TOP:
566  cvp->crv_c1 = cvp->crv_c2 = 0;
567  /* any tangent direction */
568  bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
569  break;
570  }
571 
572 }
573 
574 
575 /**
576  * For a hit on the surface of an ehy, return the (u, v) coordinates
577  * of the hit point, 0 <= u, v <= 1.
578  *
579  * u = azimuth
580  * v = elevation
581  */
582 void
583 rt_ehy_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
584 {
585  struct ehy_specific *ehy =
586  (struct ehy_specific *)stp->st_specific;
587  vect_t work;
588  vect_t pprime;
589  fastf_t len;
590 
591  if (ap) RT_CK_APPLICATION(ap);
592 
593  /*
594  * hit_point is on surface; project back to unit ehy, creating a
595  * vector from vertex to hit point.
596  */
597  VSUB2(work, hitp->hit_point, ehy->ehy_V);
598  MAT4X3VEC(pprime, ehy->ehy_SoR, work);
599 
600  switch (hitp->hit_surfno) {
601  case EHY_NORM_BODY:
602  /* top plate, polar coords */
603  if (ZERO(pprime[Z] + 1.0)) { /* i.e., == -1.0 */
604  /* bottom pt of body */
605  uvp->uv_u = 0;
606  } else {
607  len = sqrt(pprime[X]*pprime[X] + pprime[Y]*pprime[Y]);
608  uvp->uv_u = acos(pprime[X]/len) * M_1_2PI;
609  }
610  uvp->uv_v = -pprime[Z];
611  break;
612  case EHY_NORM_TOP:
613  /* top plate, polar coords */
614  len = sqrt(pprime[X]*pprime[X] + pprime[Y]*pprime[Y]);
615  uvp->uv_u = acos(pprime[X]/len) * M_1_2PI;
616  uvp->uv_v = 1.0 - len;
617  break;
618  }
619  /* Handle other half of acos() domain */
620  if (pprime[Y] < 0)
621  uvp->uv_u = 1.0 - uvp->uv_u;
622 
623  /* uv_du should be relative to rotation, uv_dv relative to height
624  */
625  uvp->uv_du = uvp->uv_dv = 0;
626 }
627 
628 
629 void
630 rt_ehy_free(struct soltab *stp)
631 {
632  struct ehy_specific *ehy =
633  (struct ehy_specific *)stp->st_specific;
634 
635  BU_PUT(ehy, struct ehy_specific);
636 }
637 
638 
639 /* Our canonical hyperbola in the Y-Z plane has equation
640  * z = +- (a/b) * sqrt(b^2 + y^2), and opens toward +Z and -Z with asymptote
641  * origin at the origin.
642  *
643  * The contour of an ehy in the plane H-R (where R is one of the ehy axes A or
644  * B) is the positive half of a hyperbola with asymptote origin at
645  * ((|H| + c)Hu), opening toward -H. We can transform this hyperbola to get an
646  * equivalent canonical hyperbola in the Y-Z plane, opening toward +Z (-H) with
647  * asymptote origin at the origin.
648  *
649  * This hyperbola passes through the point (r, |H| + a) (where r = |A| or |B|,
650  * and a = c). If we plug the point (r, |H| + a) into our canonical equation,
651  * we can derive b from |H|, a, and r.
652  *
653  * |H| + a = (a/b) * sqrt(b^2 + r^2)
654  * (|H| + a) / a = b * sqrt(b^2 + r^2)
655  * (|H| + a)^2 / a^2 = 1 + (r^2 / b^2)
656  * ((|H| + a)^2 - a^2) / a^2 = r^2 / b^2
657  * (a^2 * r^2) / ((|H| + a)^2 - a^2) = b^2
658  * (ar) / sqrt((|H| + a)^2 - a^2) = b
659  * (ar) / sqrt(|H| (|H| + 2a)) = b
660  */
661 static fastf_t
662 ehy_hyperbola_b(fastf_t mag_h, fastf_t c, fastf_t r)
663 {
664  return (c * r) / sqrt(mag_h * (mag_h + 2.0 * c));
665 }
666 
667 
668 /* The contour of an ehy in the plane H-R (where R is one of the ehy axes A or
669  * B) is the positive half of a hyperbola with asymptote origin at
670  * ((|H| + c)Hu), opening toward -H. We can transform this hyperbola to get an
671  * equivalent hyperbola in the Y-Z plane, opening toward +Z (-H) with asymptote
672  * origin at (0, -(|H| + c)).
673  *
674  * The part of this hyperbola that passes between (0, -(|H| + c)) and (r, 0)
675  * (r = |A| or |B|) is approximated by num_points points (including (0, -|H|)
676  * and (r, 0)).
677  *
678  * The constructed point list is returned (NULL returned on error). Because the
679  * above transformation puts the ehy vertex at the origin and the hyperbola
680  * asymptote origin at (0, -|H| + c), multiplying the z values by -1 gives
681  * corresponding distances along the ehy height vector H.
682  */
683 static struct rt_pt_node *
684 ehy_hyperbolic_curve(fastf_t mag_h, fastf_t c, fastf_t r, int num_points)
685 {
686  int count;
687  struct rt_pt_node *curve;
688 
689  BU_ALLOC(curve, struct rt_pt_node);
690  BU_ALLOC(curve->next, struct rt_pt_node);
691 
692  curve->next->next = NULL;
693  VSET(curve->p, 0, 0, -mag_h);
694  VSET(curve->next->p, 0, r, 0);
695 
696  count = approximate_hyperbolic_curve(curve, c, ehy_hyperbola_b(mag_h, c, r), num_points - 2);
697 
698  if (count != (num_points - 2)) {
699  return NULL;
700  }
701 
702  return curve;
703 }
704 
705 
706 /* The contour of an ehy in the plane H-R (where R is one of the ehy axes A or
707  * B) is the positive half of a hyperbola with asymptote origin at
708  * ((|H| + c)Hu), opening toward -H. We can transform this hyperbola to get an
709  * equivalent hyperbola in the Y-Z plane, with asymptote origin at
710  * (0, |H| + a) (a = c) opening toward +Z.
711  *
712  * The equation for this hyperbola is a variant of the equation for our
713  * canonical hyperbola in the Y-Z plane (z = (a/b) * sqrt(y^2 + b^2)):
714  * z = (|H| + a) - (a/b) * sqrt(y^2 + b^2)
715  *
716  * Solving this equation for y yields:
717  * y = (b/a) * sqrt((|H| + a - z)^2 - a^2)
718  *
719  * Substituting b = (ar) / sqrt(|H| (|H| + 2a)) (see above comment):
720  *
721  * y = (r / sqrt(|H| (|H| + 2a))) * sqrt((|H| + a - z)^2 - a^2)
722  * = r * sqrt(((|H| + a - z)^2 - a^2) / (|H| (|H| + 2a))))
723  */
724 static fastf_t
725 ehy_hyperbola_y(fastf_t mag_H, fastf_t c, fastf_t r, fastf_t z)
726 {
727  fastf_t n, d;
728 
729  n = pow(mag_H + c - z, 2) - c * c;
730  d = mag_H * (mag_H + 2.0 * c);
731 
732  return r * sqrt(n / d);
733 }
734 
735 
736 /* Plot the elliptical cross section of the given ehy at distance h along the
737  * ehy height vector (h >= 0, h <= |H|) consisting of num_points points.
738  */
739 static void
740 ehy_plot_ellipse(
741  struct bu_list *vhead,
742  struct rt_ehy_internal *ehy,
743  fastf_t h,
744  fastf_t num_points)
745 {
746  fastf_t mag_H;
747  vect_t V, Hu, Au, Bu, A, B, cross_section_plane;
748 
749  VMOVE(V, ehy->ehy_V);
750 
751  mag_H = MAGNITUDE(ehy->ehy_H);
752  VSCALE(Hu, ehy->ehy_H, 1.0 / mag_H);
753 
754  VMOVE(Au, ehy->ehy_Au);
755  VCROSS(Bu, Au, Hu);
756 
757  /* calculate semi-major and semi-minor axis for the elliptical
758  * cross-section at distance h along H
759  */
760  VSCALE(A, Au, ehy_hyperbola_y(mag_H, ehy->ehy_c, ehy->ehy_r1, h));
761  VSCALE(B, Bu, ehy_hyperbola_y(mag_H, ehy->ehy_c, ehy->ehy_r2, h));
762  VJOIN1(cross_section_plane, V, h, Hu);
763 
764  plot_ellipse(vhead, cross_section_plane, A, B, num_points);
765 }
766 
767 
768 static void
769 ehy_plot_hyperbola(
770  struct bu_list *vhead,
771  struct rt_ehy_internal *ehy,
772  struct rt_pt_node *pts,
773  vect_t Ru,
774  fastf_t r)
775 {
776  point_t p;
777  vect_t ehy_V, Hu;
778  fastf_t mag_H, c, z;
779  struct rt_pt_node *node;
780 
781  VMOVE(ehy_V, ehy->ehy_V);
782  mag_H = MAGNITUDE(ehy->ehy_H);
783  VSCALE(Hu, ehy->ehy_H, 1.0 / mag_H);
784  c = ehy->ehy_c;
785 
786  z = pts->p[Z];
787  VJOIN2(p, ehy_V, ehy_hyperbola_y(mag_H, c, r, -z), Ru, -z, Hu);
788  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_MOVE);
789 
790  node = pts->next;
791  while (node != NULL) {
792  z = node->p[Z];
793  VJOIN2(p, ehy_V, ehy_hyperbola_y(mag_H, c, r, -z), Ru, -z, Hu);
794 
795  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_DRAW);
796 
797  node = node->next;
798  }
799 }
800 
801 
802 static int
803 ehy_curve_points(
804  const struct rt_ehy_internal *ehy,
805  const struct rt_view_info *info)
806 {
807  fastf_t avg_r, approx_curve_len;
808  point_t p1, p2;
809 
810  avg_r = (ehy->ehy_r1 + ehy->ehy_r2) / 2.0;
811 
812  VADD2(p1, ehy->ehy_V, ehy->ehy_H);
813  VJOIN1(p2, ehy->ehy_V, avg_r, ehy->ehy_Au);
814 
815  approx_curve_len = 2.0 * DIST_PT_PT(p1, p2);
816 
817  return approx_curve_len / info->point_spacing;
818 }
819 
820 
821 static int
822 ehy_ellipse_points(
823  const struct rt_ehy_internal *ehy,
824  const struct rt_view_info *info)
825 {
826  fastf_t avg_radius, avg_circumference;
827 
828  avg_radius = (ehy->ehy_r1 + ehy->ehy_r2) / 2.0;
829  avg_circumference = M_2PI * avg_radius;
830 
831  return avg_circumference / info->point_spacing;
832 }
833 
834 
835 int
836 rt_ehy_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
837 {
838  vect_t ehy_H, Hu, Au, Bu;
839  fastf_t mag_H, z, z_step, c, r1, r2;
840  int i, num_curve_points, num_ellipse_points, num_curves;
841  struct rt_ehy_internal *ehy;
842  struct rt_pt_node *pts_r1, *pts_r2, *node, *node1, *node2;
843 
844  BU_CK_LIST_HEAD(info->vhead);
845  RT_CK_DB_INTERNAL(ip);
846  ehy = (struct rt_ehy_internal *)ip->idb_ptr;
847  RT_EHY_CK_MAGIC(ehy);
848 
849  num_curve_points = ehy_curve_points(ehy, info);
850 
851  if (num_curve_points < 3) {
852  num_curve_points = 3;
853  }
854 
855  num_ellipse_points = ehy_ellipse_points(ehy, info);
856 
857  if (num_ellipse_points < 6) {
858  num_ellipse_points = 6;
859  }
860 
861  VMOVE(ehy_H, ehy->ehy_H);
862 
863  mag_H = MAGNITUDE(ehy_H);
864  VSCALE(Hu, ehy->ehy_H, 1.0 / mag_H);
865 
866  VMOVE(Au, ehy->ehy_Au);
867  VCROSS(Bu, Au, Hu);
868 
869  r1 = ehy->ehy_r1;
870  r2 = ehy->ehy_r2;
871  c = ehy->ehy_c;
872 
873  pts_r1 = ehy_hyperbolic_curve(mag_H, c, r1, num_curve_points);
874  pts_r2 = ehy_hyperbolic_curve(mag_H, c, r2, num_curve_points);
875 
876  if (pts_r1 == NULL || pts_r2 == NULL) {
877  return -1;
878  }
879 
880  num_curves = mag_H / info->curve_spacing;
881  if (num_curves < 2) {
882  num_curves = 2;
883  }
884 
885  z_step = mag_H / num_curves;
886  z = 0.0;
887  for (i = 0; i < num_curves; ++i) {
888  ehy_plot_ellipse(info->vhead, ehy, z, num_ellipse_points);
889  z += z_step;
890  }
891 
892  ehy_plot_hyperbola(info->vhead, ehy, pts_r1, Au, r1);
893  ehy_plot_hyperbola(info->vhead, ehy, pts_r1, Au, -r1);
894  ehy_plot_hyperbola(info->vhead, ehy, pts_r1, Bu, r2);
895  ehy_plot_hyperbola(info->vhead, ehy, pts_r1, Bu, -r2);
896 
897  node1 = pts_r1;
898  node2 = pts_r2;
899  for (i = 0; i < num_curve_points; ++i) {
900  node = node1;
901  node1 = node1->next;
902  bu_free(node, "rt_pt_node");
903 
904  node = node2;
905  node2 = node2->next;
906  bu_free(node, "rt_pt_node");
907  }
908 
909  return 0;
910 }
911 
912 
913 int
914 rt_ehy_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))
915 {
916  fastf_t c, dtol, mag_h, ntol, r1, r2;
917  fastf_t **ellipses, theta_prev, theta_new;
918  int *pts_dbl, i, j, nseg;
919  int jj, na, nb, nell, recalc_b;
920  mat_t R;
921  mat_t invR;
922  point_t p1;
923  struct rt_pt_node *pos_a, *pos_b, *pts_a, *pts_b;
924  struct rt_ehy_internal *xip;
925  vect_t A, Au, B, Bu, Hu, V, Work;
926 
927  BU_CK_LIST_HEAD(vhead);
928  RT_CK_DB_INTERNAL(ip);
929  xip = (struct rt_ehy_internal *)ip->idb_ptr;
930 
931  if (!ehy_is_valid(xip)) {
932  return -2;
933  }
934 
935  mag_h = MAGNITUDE(xip->ehy_H);
936  r1 = xip->ehy_r1;
937  r2 = xip->ehy_r2;
938  c = xip->ehy_c;
939 
940  /* make unit vectors in A, H, and BxH directions */
941  VMOVE(Hu, xip->ehy_H);
942  VUNITIZE(Hu);
943  VMOVE(Au, xip->ehy_Au);
944  VCROSS(Bu, Au, Hu);
945 
946  /* Compute R and Rinv matrices */
947  MAT_IDN(R);
948  VREVERSE(&R[0], Bu);
949  VMOVE(&R[4], Au);
950  VREVERSE(&R[8], Hu);
951  bn_mat_trn(invR, R); /* inv of rot mat is trn */
952 
953  dtol = primitive_get_absolute_tolerance(ttol, 2.0 * xip->ehy_r2);
954 
955  /* stay below ntol to ensure normal tolerance */
956  ntol = M_PI;
957  if (ttol->norm > 0.0) {
958  ntol = ttol->norm;
959  }
960 
961  /*
962  * build ehy from 2 hyperbolas
963  */
964 
965  /* approximate positive half of hyperbola along semi-minor axis */
966  BU_ALLOC(pts_b, struct rt_pt_node);
967  BU_ALLOC(pts_b->next, struct rt_pt_node);
968 
969  pts_b->next->next = NULL;
970  VSET(pts_b->p, 0, 0, -mag_h);
971  VSET(pts_b->next->p, 0, r2, 0);
972  /* 2 endpoints in 1st approximation */
973  nb = 2;
974  /* recursively break segment 'til within error tolerances */
975  nb += rt_mk_hyperbola(pts_b, r2, mag_h, c, dtol, ntol);
976  nell = nb - 1; /* # of ellipses needed */
977 
978  /* construct positive half of hyperbola along semi-major axis of
979  * ehy using same z coords as hyperbola along semi-minor axis
980  */
981  BU_ALLOC(pts_a, struct rt_pt_node);
982  VMOVE(pts_a->p, pts_b->p); /* 1st pt is the apex */
983  pts_a->next = NULL;
984  pos_b = pts_b->next;
985  pos_a = pts_a;
986  while (pos_b) {
987  /* copy node from b_hyperbola to a_hyperbola */
988  BU_ALLOC(pos_a->next, struct rt_pt_node);
989  pos_a = pos_a->next;
990  pos_a->p[Z] = pos_b->p[Z];
991  /* at given z, find y on hyperbola */
992  pos_a->p[Y] = r1
993  * sqrt(
994  ((pos_b->p[Z] + mag_h + c)
995  * (pos_b->p[Z] + mag_h + c) - c*c)
996  / (mag_h*(mag_h + 2*c))
997  );
998  pos_a->p[X] = 0;
999  pos_b = pos_b->next;
1000  }
1001  pos_a->next = NULL;
1002 
1003  /* see if any segments need further breaking up */
1004  recalc_b = 0;
1005  pos_a = pts_a;
1006  while (pos_a->next) {
1007  na = rt_mk_hyperbola(pos_a, r1, mag_h, c, dtol, ntol);
1008  if (na != 0) {
1009  recalc_b = 1;
1010  nell += na;
1011  }
1012  pos_a = pos_a->next;
1013  }
1014  /* if any were broken, recalculate hyperbola 'a' */
1015  if (recalc_b) {
1016  /* free mem for old approximation of hyperbola */
1017  pos_b = pts_b;
1018  while (pos_b) {
1019  struct rt_pt_node *next;
1020 
1021  /* get next node before freeing */
1022  next = pos_b->next;
1023  bu_free((char *)pos_b, "rt_pt_node");
1024  pos_b = next;
1025  }
1026  /* construct hyperbola along semi-major axis of ehy using same
1027  * z coords as parab along semi-minor axis
1028  */
1029  BU_ALLOC(pts_b, struct rt_pt_node);
1030  pts_b->p[Z] = pts_a->p[Z];
1031  pts_b->next = NULL;
1032  pos_a = pts_a->next;
1033  pos_b = pts_b;
1034  while (pos_a) {
1035  /* copy node from a_hyperbola to b_hyperbola */
1036  BU_ALLOC(pos_b->next, struct rt_pt_node);
1037  pos_b = pos_b->next;
1038  pos_b->p[Z] = pos_a->p[Z];
1039  /* at given z, find y on hyperbola */
1040  pos_b->p[Y] = r2
1041  * sqrt(
1042  ((pos_a->p[Z] + mag_h + c)
1043  * (pos_a->p[Z] + mag_h + c) - c*c)
1044  / (mag_h*(mag_h + 2*c))
1045  );
1046  pos_b->p[X] = 0;
1047  pos_a = pos_a->next;
1048  }
1049  pos_b->next = NULL;
1050  }
1051 
1052  /* make array of ptrs to ehy ellipses */
1053  ellipses = (fastf_t **)bu_malloc(nell * sizeof(fastf_t *), "fastf_t ell[]");
1054  /* keep track of whether pts in each ellipse are doubled or not */
1055  pts_dbl = (int *)bu_malloc(nell * sizeof(int), "dbl ints");
1056 
1057  /* make ellipses at each z level */
1058  i = 0;
1059  nseg = 0;
1060  theta_prev = M_2PI;
1061  pos_a = pts_a->next; /* skip over apex of ehy */
1062  pos_b = pts_b->next;
1063  while (pos_a) {
1064  VSCALE(A, Au, pos_a->p[Y]); /* semimajor axis */
1065  VSCALE(B, Bu, pos_b->p[Y]); /* semiminor axis */
1066  VJOIN1(V, xip->ehy_V, -pos_a->p[Z], Hu);
1067 
1068  VSET(p1, 0., pos_b->p[Y], 0.);
1069  theta_new = ell_angle(p1, pos_a->p[Y], pos_b->p[Y], dtol, ntol);
1070  if (nseg == 0) {
1071  nseg = (int)(M_2PI / theta_new) + 1;
1072  pts_dbl[i] = 0;
1073  } else if (theta_new < theta_prev) {
1074  nseg *= 2;
1075  pts_dbl[i] = 1;
1076  } else
1077  pts_dbl[i] = 0;
1078  theta_prev = theta_new;
1079 
1080  ellipses[i] = (fastf_t *)bu_malloc(3*(nseg+1)*sizeof(fastf_t),
1081  "pts ell");
1082  rt_ell(ellipses[i], V, A, B, nseg);
1083 
1084  i++;
1085  pos_a = pos_a->next;
1086  pos_b = pos_b->next;
1087  }
1088 
1089  /* Draw the top ellipse */
1090  RT_ADD_VLIST(vhead,
1091  &ellipses[nell-1][(nseg-1)*ELEMENTS_PER_VECT],
1093  for (i = 0; i < nseg; i++) {
1094  RT_ADD_VLIST(vhead,
1095  &ellipses[nell-1][i*ELEMENTS_PER_VECT],
1097  }
1098 
1099  /* connect ellipses */
1100  for (i = nell-2; i >= 0; i--) {
1101  /* skip top ellipse */
1102  int bottom, top;
1103 
1104  top = i + 1;
1105  bottom = i;
1106  if (pts_dbl[top])
1107  nseg /= 2; /* # segs in 'bottom' ellipse */
1108 
1109  /* Draw the current ellipse */
1110  RT_ADD_VLIST(vhead,
1111  &ellipses[bottom][(nseg-1)*ELEMENTS_PER_VECT],
1113  for (j = 0; j < nseg; j++) {
1114  RT_ADD_VLIST(vhead,
1115  &ellipses[bottom][j*ELEMENTS_PER_VECT],
1117  }
1118 
1119  /* make connections between ellipses */
1120  for (j = 0; j < nseg; j++) {
1121  if (pts_dbl[top])
1122  jj = j + j; /* top ellipse index */
1123  else
1124  jj = j;
1125  RT_ADD_VLIST(vhead,
1126  &ellipses[bottom][j*ELEMENTS_PER_VECT],
1128  RT_ADD_VLIST(vhead,
1129  &ellipses[top][jj*ELEMENTS_PER_VECT],
1131  }
1132  }
1133 
1134  VADD2(Work, xip->ehy_V, xip->ehy_H);
1135  for (i = 0; i < nseg; i++) {
1136  /* Draw connector */
1137  RT_ADD_VLIST(vhead, Work, BN_VLIST_LINE_MOVE);
1138  RT_ADD_VLIST(vhead,
1139  &ellipses[0][i*ELEMENTS_PER_VECT],
1141  }
1142 
1143  /* free mem */
1144  for (i = 0; i < nell; i++) {
1145  bu_free((char *)ellipses[i], "pts ell");
1146  }
1147  bu_free((char *)ellipses, "fastf_t ell[]");
1148  bu_free((char *)pts_dbl, "dbl ints");
1149 
1150  return 0;
1151 }
1152 
1153 
1154 /**
1155  * Returns -
1156  * -1 failure
1157  * 0 OK. *r points to nmgregion that holds this tessellation.
1158  */
1159 int
1160 rt_ehy_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
1161 {
1162  fastf_t c, dtol, mag_h, ntol, r1, r2, cprime;
1163  fastf_t **ellipses, theta_prev, theta_new;
1164  int *pts_dbl, face, i, j, nseg;
1165  int jj, na, nb, nell, recalc_b;
1166  mat_t R;
1167  mat_t invR;
1168  mat_t invRoS;
1169  mat_t S;
1170  mat_t SoR;
1171  struct rt_ehy_internal *xip;
1172  point_t p1;
1173  struct rt_pt_node *pos_a, *pos_b, *pts_a, *pts_b;
1174  struct shell *s;
1175  struct faceuse **outfaceuses = NULL;
1176  struct faceuse *fu_top;
1177  struct loopuse *lu;
1178  struct edgeuse *eu;
1179  struct vertex *vertp[3];
1180  struct vertex ***vells = (struct vertex ***)NULL;
1181  vect_t A, Au, B, Bu, Hu, V;
1182  struct bu_ptbl vert_tab;
1183 
1184  RT_CK_DB_INTERNAL(ip);
1185  xip = (struct rt_ehy_internal *)ip->idb_ptr;
1186 
1187  if (!ehy_is_valid(xip)) {
1188  return 1;
1189  }
1190 
1191  mag_h = MAGNITUDE(xip->ehy_H);
1192  c = xip->ehy_c;
1193  cprime = c / mag_h;
1194  r1 = xip->ehy_r1;
1195  r2 = xip->ehy_r2;
1196 
1197  /* make unit vectors in A, H, and BxH directions */
1198  VMOVE(Hu, xip->ehy_H);
1199  VUNITIZE(Hu);
1200  VMOVE(Au, xip->ehy_Au);
1201  VCROSS(Bu, Au, Hu);
1202 
1203  /* Compute R and Rinv matrices */
1204  MAT_IDN(R);
1205  VREVERSE(&R[0], Bu);
1206  VMOVE(&R[4], Au);
1207  VREVERSE(&R[8], Hu);
1208  bn_mat_trn(invR, R); /* inv of rot mat is trn */
1209 
1210  /* Compute S */
1211  MAT_IDN(S);
1212  S[ 0] = 1.0/r2;
1213  S[ 5] = 1.0/r1;
1214  S[10] = 1.0/mag_h;
1215 
1216  /* Compute SoR and invRoS */
1217  bn_mat_mul(SoR, S, R);
1218  bn_mat_mul(invRoS, invR, S);
1219 
1220  dtol = primitive_get_absolute_tolerance(ttol, 2.0 * xip->ehy_r2);
1221 
1222  /* stay below ntol to ensure normal tolerance */
1223  ntol = M_PI;
1224  if (ttol->norm > 0.0) {
1225  ntol = ttol->norm;
1226  }
1227 
1228  /*
1229  * build ehy from 2 hyperbolas
1230  */
1231 
1232  /* approximate positive half of hyperbola along semi-minor axis */
1233  BU_ALLOC(pts_b, struct rt_pt_node);
1234  BU_ALLOC(pts_b->next, struct rt_pt_node);
1235 
1236  pts_b->next->next = NULL;
1237  VSET(pts_b->p, 0, 0, -mag_h);
1238  VSET(pts_b->next->p, 0, r2, 0);
1239  /* 2 endpoints in 1st approximation */
1240  nb = 2;
1241  /* recursively break segment 'til within error tolerances */
1242  nb += rt_mk_hyperbola(pts_b, r2, mag_h, c, dtol, ntol);
1243  nell = nb - 1; /* # of ellipses needed */
1244 
1245  /* construct positive half of hyperbola along semi-major axis of
1246  * ehy using same z coords as parab along semi-minor axis
1247  */
1248  BU_ALLOC(pts_a, struct rt_pt_node);
1249  VMOVE(pts_a->p, pts_b->p); /* 1st pt is the apex */
1250  pts_a->next = NULL;
1251  pos_b = pts_b->next;
1252  pos_a = pts_a;
1253  while (pos_b) {
1254  /* copy node from b_hyperbola to a_hyperbola */
1255  BU_ALLOC(pos_a->next, struct rt_pt_node);
1256  pos_a = pos_a->next;
1257  pos_a->p[Z] = pos_b->p[Z];
1258  /* at given z, find y on hyperbola */
1259  pos_a->p[Y] = r1
1260  * sqrt(
1261  ((pos_b->p[Z] + mag_h + c)
1262  * (pos_b->p[Z] + mag_h + c) - c*c)
1263  / (mag_h*(mag_h + 2*c))
1264  );
1265  pos_a->p[X] = 0;
1266  pos_b = pos_b->next;
1267  }
1268  pos_a->next = NULL;
1269 
1270  /* see if any segments need further breaking up */
1271  recalc_b = 0;
1272  pos_a = pts_a;
1273  while (pos_a->next) {
1274  na = rt_mk_hyperbola(pos_a, r1, mag_h, c, dtol, ntol);
1275  if (na != 0) {
1276  recalc_b = 1;
1277  nell += na;
1278  }
1279  pos_a = pos_a->next;
1280  }
1281  /* if any were broken, recalculate hyperbola 'a' */
1282  if (recalc_b) {
1283  /* free mem for old approximation of hyperbola */
1284  pos_b = pts_b;
1285  while (pos_b) {
1286  struct rt_pt_node *tmp;
1287 
1288  tmp = pos_b->next;
1289  bu_free((char *)pos_b, "rt_pt_node");
1290  pos_b = tmp;
1291  }
1292  /* construct hyperbola along semi-major axis of ehy using same
1293  * z coords as parab along semi-minor axis
1294  */
1295  BU_ALLOC(pts_b, struct rt_pt_node);
1296  pts_b->p[Z] = pts_a->p[Z];
1297  pts_b->next = NULL;
1298  pos_a = pts_a->next;
1299  pos_b = pts_b;
1300  while (pos_a) {
1301  /* copy node from a_hyperbola to b_hyperbola */
1302  BU_ALLOC(pos_b->next, struct rt_pt_node);
1303  pos_b = pos_b->next;
1304  pos_b->p[Z] = pos_a->p[Z];
1305  /* at given z, find y on hyperbola */
1306  pos_b->p[Y] = r2
1307  * sqrt(
1308  ((pos_a->p[Z] + mag_h + c)
1309  * (pos_a->p[Z] + mag_h + c) - c*c)
1310  / (mag_h*(mag_h + 2*c))
1311  );
1312  pos_b->p[X] = 0;
1313  pos_a = pos_a->next;
1314  }
1315  pos_b->next = NULL;
1316  }
1317 
1318  /* make array of ptrs to ehy ellipses */
1319  ellipses = (fastf_t **)bu_malloc(nell * sizeof(fastf_t *), "fastf_t ell[]");
1320 
1321  /* keep track of whether pts in each ellipse are doubled or not */
1322  pts_dbl = (int *)bu_malloc(nell * sizeof(int), "dbl ints");
1323 
1324  /* make ellipses at each z level */
1325  i = 0;
1326  nseg = 0;
1327  theta_prev = M_2PI;
1328  pos_a = pts_a->next; /* skip over apex of ehy */
1329  pos_b = pts_b->next;
1330  while (pos_a) {
1331  VSCALE(A, Au, pos_a->p[Y]); /* semimajor axis */
1332  VSCALE(B, Bu, pos_b->p[Y]); /* semiminor axis */
1333  VJOIN1(V, xip->ehy_V, -pos_a->p[Z], Hu);
1334 
1335  VSET(p1, 0., pos_b->p[Y], 0.);
1336  theta_new = ell_angle(p1, pos_a->p[Y], pos_b->p[Y], dtol, ntol);
1337  if (nseg == 0) {
1338  nseg = (int)(M_2PI / theta_new) + 1;
1339  pts_dbl[i] = 0;
1340  /* maximum number of faces needed for ehy */
1341  face = nseg*(1 + 3*((1 << (nell-1)) - 1));
1342  /* array for each triangular face */
1343  outfaceuses = (struct faceuse **)
1344  bu_malloc((face+1) * sizeof(struct faceuse *), "ehy: *outfaceuses[]");
1345  } else if (theta_new < theta_prev) {
1346  nseg *= 2;
1347  pts_dbl[i] = 1;
1348  } else {
1349  pts_dbl[i] = 0;
1350  }
1351  theta_prev = theta_new;
1352 
1353  ellipses[i] = (fastf_t *)bu_malloc(3*(nseg+1)*sizeof(fastf_t),
1354  "pts ell");
1355  rt_ell(ellipses[i], V, A, B, nseg);
1356 
1357  i++;
1358  pos_a = pos_a->next;
1359  pos_b = pos_b->next;
1360  }
1361 
1362  /*
1363  * put ehy geometry into nmg data structures
1364  */
1365 
1366  *r = nmg_mrsv(m); /* Make region, empty shell, vertex */
1367  s = BU_LIST_FIRST(shell, &(*r)->s_hd);
1368 
1369  /* vertices of ellipses of ehy */
1370  vells = (struct vertex ***)
1371  bu_malloc(nell*sizeof(struct vertex **), "vertex [][]");
1372  j = nseg;
1373  for (i = nell-1; i >= 0; i--) {
1374  vells[i] = (struct vertex **)
1375  bu_malloc(j*sizeof(struct vertex *), "vertex []");
1376  if (i && pts_dbl[i])
1377  j /= 2;
1378  }
1379 
1380  /* top face of ehy */
1381  for (i = 0; i < nseg; i++)
1382  vells[nell-1][i] = (struct vertex *)NULL;
1383  face = 0;
1384  BU_ASSERT_PTR(outfaceuses, !=, NULL);
1385  if ((outfaceuses[face++] = nmg_cface(s, vells[nell-1], nseg)) == 0) {
1386  bu_log("rt_ehy_tess() failure, top face\n");
1387  goto fail;
1388  }
1389  fu_top = outfaceuses[0];
1390 
1391  /* Mark edges of this face as real, this is the only real edge */
1392  for (BU_LIST_FOR(lu, loopuse, &outfaceuses[0]->lu_hd)) {
1393  NMG_CK_LOOPUSE(lu);
1394 
1395  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
1396  continue;
1397  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
1398  struct edge *e;
1399 
1400  NMG_CK_EDGEUSE(eu);
1401  e = eu->e_p;
1402  NMG_CK_EDGE(e);
1403  e->is_real = 1;
1404  }
1405  }
1406 
1407  for (i = 0; i < nseg; i++) {
1408  NMG_CK_VERTEX(vells[nell-1][i]);
1409  nmg_vertex_gv(vells[nell-1][i], &ellipses[nell-1][3*i]);
1410  }
1411 
1412  /* connect ellipses with triangles */
1413  for (i = nell-2; i >= 0; i--) {
1414  /* skip top ellipse */
1415  int bottom, top;
1416 
1417  top = i + 1;
1418  bottom = i;
1419  if (pts_dbl[top])
1420  nseg /= 2; /* # segs in 'bottom' ellipse */
1421  vertp[0] = (struct vertex *)0;
1422 
1423  /* make triangular faces */
1424  for (j = 0; j < nseg; j++) {
1425  jj = j + j; /* top ellipse index */
1426 
1427  if (pts_dbl[top]) {
1428  /* first triangle */
1429  vertp[1] = vells[top][jj+1];
1430  vertp[2] = vells[top][jj];
1431  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1432  bu_log("rt_ehy_tess() failure\n");
1433  goto fail;
1434  }
1435  if (j == 0)
1436  vells[bottom][j] = vertp[0];
1437 
1438  /* second triangle */
1439  vertp[2] = vertp[1];
1440  if (j == nseg-1)
1441  vertp[1] = vells[bottom][0];
1442  else
1443  vertp[1] = (struct vertex *)0;
1444  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1445  bu_log("rt_ehy_tess() failure\n");
1446  goto fail;
1447  }
1448  if (j != nseg-1)
1449  vells[bottom][j+1] = vertp[1];
1450 
1451  /* third triangle */
1452  vertp[0] = vertp[1];
1453  if (j == nseg-1)
1454  vertp[1] = vells[top][0];
1455  else
1456  vertp[1] = vells[top][jj+2];
1457  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1458  bu_log("rt_ehy_tess() failure\n");
1459  goto fail;
1460  }
1461  } else {
1462  /* first triangle */
1463  if (j == nseg-1)
1464  vertp[1] = vells[top][0];
1465  else
1466  vertp[1] = vells[top][j+1];
1467  vertp[2] = vells[top][j];
1468  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1469  bu_log("rt_ehy_tess() failure\n");
1470  goto fail;
1471  }
1472  if (j == 0)
1473  vells[bottom][j] = vertp[0];
1474 
1475  /* second triangle */
1476  vertp[2] = vertp[0];
1477  if (j == nseg-1)
1478  vertp[0] = vells[bottom][0];
1479  else
1480  vertp[0] = (struct vertex *)0;
1481  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1482  bu_log("rt_ehy_tess() failure\n");
1483  goto fail;
1484  }
1485  if (j != nseg-1)
1486  vells[bottom][j+1] = vertp[0];
1487  }
1488  }
1489 
1490  /* associate geometry with each vertex */
1491  for (j = 0; j < nseg; j++) {
1492  NMG_CK_VERTEX(vells[bottom][j]);
1493  nmg_vertex_gv(vells[bottom][j],
1494  &ellipses[bottom][3*j]);
1495  }
1496  }
1497 
1498  /* connect bottom of ellipse to apex of ehy */
1499  VADD2(V, xip->ehy_V, xip->ehy_H);
1500  vertp[0] = (struct vertex *)0;
1501  vertp[1] = vells[0][1];
1502  vertp[2] = vells[0][0];
1503  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1504  bu_log("rt_ehy_tess() failure\n");
1505  goto fail;
1506  }
1507  /* associate geometry with topology */
1508  NMG_CK_VERTEX(vertp[0]);
1509  nmg_vertex_gv(vertp[0], (fastf_t *)V);
1510  /* create rest of faces around apex */
1511  for (i = 1; i < nseg; i++) {
1512  vertp[2] = vertp[1];
1513  if (i == nseg-1)
1514  vertp[1] = vells[0][0];
1515  else
1516  vertp[1] = vells[0][i+1];
1517  if ((outfaceuses[face++] = nmg_cface(s, vertp, 3)) == 0) {
1518  bu_log("rt_ehy_tess() failure\n");
1519  goto fail;
1520  }
1521  }
1522 
1523  /* Associate the face geometry */
1524  for (i = 0; i < face; i++) {
1525  if (nmg_fu_planeeqn(outfaceuses[i], tol) < 0)
1526  goto fail;
1527  }
1528 
1529  /* Glue the edges of different outward pointing face uses together */
1530  nmg_gluefaces(outfaceuses, face, tol);
1531 
1532  /* Compute "geometry" for region and shell */
1533  nmg_region_a(*r, tol);
1534 
1535  /* XXX just for testing, to make up for loads of triangles ... */
1536  nmg_shell_coplanar_face_merge(s, tol, 1);
1537 
1538  /* free mem */
1539  bu_free((char *)outfaceuses, "faceuse []");
1540  for (i = 0; i < nell; i++) {
1541  bu_free((char *)ellipses[i], "pts ell");
1542  bu_free((char *)vells[i], "vertex []");
1543  }
1544  bu_free((char *)ellipses, "fastf_t ell[]");
1545  bu_free((char *)vells, "vertex [][]");
1546 
1547  /* Assign vertexuse normals */
1548  nmg_vertex_tabulate(&vert_tab, &s->l.magic);
1549  for (i = 0; i < BU_PTBL_END(&vert_tab); i++) {
1550  point_t pt_prime, tmp_pt;
1551  vect_t norm, rev_norm, tmp_vect;
1552  struct vertex_g *vg;
1553  struct vertex *v;
1554  struct vertexuse *vu;
1555 
1556  v = (struct vertex *)BU_PTBL_GET(&vert_tab, i);
1557  NMG_CK_VERTEX(v);
1558  vg = v->vg_p;
1559  NMG_CK_VERTEX_G(vg);
1560 
1561  VSUB2(tmp_pt, vg->coord, xip->ehy_V);
1562  MAT4X3VEC(pt_prime, SoR, tmp_pt);
1563  VSET(tmp_vect, pt_prime[X]*(2*cprime+1), pt_prime[Y]*(2*cprime+1), -(pt_prime[Z]+cprime+1));
1564  MAT4X3VEC(norm, invRoS, tmp_vect);
1565  VUNITIZE(norm);
1566  VREVERSE(rev_norm, norm);
1567 
1568  for (BU_LIST_FOR(vu, vertexuse, &v->vu_hd)) {
1569  struct faceuse *fu;
1570 
1571  NMG_CK_VERTEXUSE(vu);
1572  fu = nmg_find_fu_of_vu(vu);
1573 
1574  /* don't assign vertexuse normals to top face (flat) */
1575  if (fu == fu_top || fu->fumate_p == fu_top)
1576  continue;
1577 
1578  NMG_CK_FACEUSE(fu);
1579  if (fu->orientation == OT_SAME)
1580  nmg_vertexuse_nv(vu, norm);
1581  else if (fu->orientation == OT_OPPOSITE)
1582  nmg_vertexuse_nv(vu, rev_norm);
1583  }
1584  }
1585 
1586  bu_ptbl_free(&vert_tab);
1587  return 0;
1588 
1589 fail:
1590  /* free mem */
1591  bu_free((char *)outfaceuses, "faceuse []");
1592  for (i = 0; i < nell; i++) {
1593  bu_free((char *)ellipses[i], "pts ell");
1594  bu_free((char *)vells[i], "vertex []");
1595  }
1596  bu_free((char *)ellipses, "fastf_t ell[]");
1597  bu_free((char *)vells, "vertex [][]");
1598 
1599  return -1;
1600 }
1601 
1602 
1603 /**
1604  * Import an EHY from the database format to the internal format.
1605  * Apply modeling transformations as well.
1606  */
1607 int
1608 rt_ehy_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
1609 {
1610  struct rt_ehy_internal *xip;
1611  union record *rp;
1612  vect_t v1, v2, v3;
1613 
1614  if (dbip) RT_CK_DBI(dbip);
1615 
1616  BU_CK_EXTERNAL(ep);
1617  rp = (union record *)ep->ext_buf;
1618  /* Check record type */
1619  if (rp->u_id != ID_SOLID) {
1620  bu_log("rt_ehy_import4: defective record\n");
1621  return -1;
1622  }
1623 
1624  RT_CK_DB_INTERNAL(ip);
1625  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1626  ip->idb_type = ID_EHY;
1627  ip->idb_meth = &OBJ[ID_EHY];
1628  BU_ALLOC(ip->idb_ptr, struct rt_ehy_internal);
1629 
1630  xip = (struct rt_ehy_internal *)ip->idb_ptr;
1631  xip->ehy_magic = RT_EHY_INTERNAL_MAGIC;
1632 
1633  /* Warning: type conversion */
1634  if (mat == NULL) mat = bn_mat_identity;
1635 
1636  if (dbip->dbi_version < 0) {
1637  flip_fastf_float(v1, &rp->s.s_values[0*3], 1, 1);
1638  flip_fastf_float(v2, &rp->s.s_values[1*3], 1, 1);
1639  flip_fastf_float(v3, &rp->s.s_values[2*3], 1, 1);
1640  } else {
1641  VMOVE(v1, &rp->s.s_values[0*3]);
1642  VMOVE(v2, &rp->s.s_values[1*3]);
1643  VMOVE(v3, &rp->s.s_values[2*3]);
1644  }
1645 
1646  MAT4X3PNT(xip->ehy_V, mat, v1);
1647  MAT4X3VEC(xip->ehy_H, mat, v2);
1648  MAT4X3VEC(xip->ehy_Au, mat, v3);
1649 
1650  VUNITIZE(xip->ehy_Au);
1651 
1652  if (dbip->dbi_version < 0) {
1653  v1[X] = flip_dbfloat(rp->s.s_values[3*3+0]);
1654  v1[Y] = flip_dbfloat(rp->s.s_values[3*3+1]);
1655  v1[Z] = flip_dbfloat(rp->s.s_values[3*3+2]);
1656  } else {
1657  v1[X] = rp->s.s_values[3*3+0];
1658  v1[Y] = rp->s.s_values[3*3+1];
1659  v1[Z] = rp->s.s_values[3*3+2];
1660  }
1661 
1662  xip->ehy_r1 = v1[X] / mat[15];
1663  xip->ehy_r2 = v1[Y] / mat[15];
1664  xip->ehy_c = v1[Z] / mat[15];
1665 
1666  if (xip->ehy_r1 <= SMALL_FASTF || xip->ehy_r2 <= SMALL_FASTF || xip->ehy_c <= SMALL_FASTF) {
1667  bu_log("rt_ehy_import4: r1, r2, or c are zero\n");
1668  bu_free((char *)ip->idb_ptr, "rt_ehy_import4: ip->idb_ptr");
1669  return -1;
1670  }
1671 
1672  return 0; /* OK */
1673 }
1674 
1675 
1676 /**
1677  * The name is added by the caller, in the usual place.
1678  */
1679 int
1680 rt_ehy_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1681 {
1682  struct rt_ehy_internal *xip;
1683  union record *ehy;
1684 
1685  if (dbip) RT_CK_DBI(dbip);
1686 
1687  RT_CK_DB_INTERNAL(ip);
1688  if (ip->idb_type != ID_EHY) return -1;
1689  xip = (struct rt_ehy_internal *)ip->idb_ptr;
1690  RT_EHY_CK_MAGIC(xip);
1691 
1692  BU_CK_EXTERNAL(ep);
1693  ep->ext_nbytes = sizeof(union record);
1694  ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "ehy external");
1695  ehy = (union record *)ep->ext_buf;
1696 
1697  ehy->s.s_id = ID_SOLID;
1698  ehy->s.s_type = EHY;
1699 
1700  if (!NEAR_EQUAL(MAGNITUDE(xip->ehy_Au), 1.0, RT_LEN_TOL)) {
1701  bu_log("rt_ehy_export4: Au not a unit vector!\n");
1702  return -1;
1703  }
1704 
1705  if (MAGNITUDE(xip->ehy_H) < RT_LEN_TOL
1706  || xip->ehy_c < RT_LEN_TOL
1707  || xip->ehy_r1 < RT_LEN_TOL
1708  || xip->ehy_r2 < RT_LEN_TOL) {
1709  bu_log("rt_ehy_export4: not all dimensions positive!\n");
1710  return -1;
1711  }
1712 
1713  if (!NEAR_ZERO(VDOT(xip->ehy_Au, xip->ehy_H), RT_DOT_TOL)) {
1714  bu_log("rt_ehy_export4: Au and H are not perpendicular!\n");
1715  return -1;
1716  }
1717 
1718  if (xip->ehy_r2 > xip->ehy_r1) {
1719  bu_log("rt_ehy_export4: semi-minor axis cannot be longer than semi-major axis!\n");
1720  return -1;
1721  }
1722 
1723  /* Warning: type conversion */
1724  VSCALE(&ehy->s.s_values[0*3], xip->ehy_V, local2mm);
1725  VSCALE(&ehy->s.s_values[1*3], xip->ehy_H, local2mm);
1726  /* don't scale ehy_Au (unit vector!!) */
1727  VMOVE(&ehy->s.s_values[2*3], xip->ehy_Au);
1728  ehy->s.s_values[3*3] = xip->ehy_r1 * local2mm;
1729  ehy->s.s_values[3*3+1] = xip->ehy_r2 * local2mm;
1730  ehy->s.s_values[3*3+2] = xip->ehy_c * local2mm;
1731 
1732  return 0;
1733 }
1734 
1735 
1736 /**
1737  * Import an EHY from the database format to the internal format.
1738  * Apply modeling transformations as well.
1739  */
1740 int
1741 rt_ehy_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
1742 {
1743  struct rt_ehy_internal *xip;
1744 
1745  /* must be double for import and export */
1746  double vec[3*4];
1747 
1748  BU_CK_EXTERNAL(ep);
1749  if (dbip) RT_CK_DBI(dbip);
1750 
1752 
1753  RT_CK_DB_INTERNAL(ip);
1754  ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
1755  ip->idb_type = ID_EHY;
1756  ip->idb_meth = &OBJ[ID_EHY];
1757  BU_ALLOC(ip->idb_ptr, struct rt_ehy_internal);
1758 
1759  xip = (struct rt_ehy_internal *)ip->idb_ptr;
1760  xip->ehy_magic = RT_EHY_INTERNAL_MAGIC;
1761 
1762  /* Convert from database (network) to internal (host) format */
1763  bu_cv_ntohd((unsigned char *)vec, ep->ext_buf, 3*4);
1764 
1765  /* Apply modeling transformations */
1766  if (mat == NULL) mat = bn_mat_identity;
1767  MAT4X3PNT(xip->ehy_V, mat, &vec[0*3]);
1768  MAT4X3VEC(xip->ehy_H, mat, &vec[1*3]);
1769  MAT4X3VEC(xip->ehy_Au, mat, &vec[2*3]);
1770  VUNITIZE(xip->ehy_Au);
1771  xip->ehy_r1 = vec[3*3] / mat[15];
1772  xip->ehy_r2 = vec[3*3+1] / mat[15];
1773  xip->ehy_c = vec[3*3+2] / mat[15];
1774 
1775  if (xip->ehy_r1 <= SMALL_FASTF || xip->ehy_r2 <= SMALL_FASTF || xip->ehy_c <= SMALL_FASTF) {
1776  bu_log("rt_ehy_import5: r1, r2, or c are zero\n");
1777  bu_free((char *)ip->idb_ptr, "rt_ehy_import5: ip->idb_ptr");
1778  return -1;
1779  }
1780 
1781  return 0; /* OK */
1782 }
1783 
1784 
1785 /**
1786  * The name is added by the caller, in the usual place.
1787  */
1788 int
1789 rt_ehy_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
1790 {
1791  struct rt_ehy_internal *xip;
1792 
1793  /* must be double for import and export */
1794  double vec[3*4];
1795 
1796  if (dbip) RT_CK_DBI(dbip);
1797 
1798  RT_CK_DB_INTERNAL(ip);
1799  if (ip->idb_type != ID_EHY) return -1;
1800  xip = (struct rt_ehy_internal *)ip->idb_ptr;
1801  RT_EHY_CK_MAGIC(xip);
1802 
1803  BU_CK_EXTERNAL(ep);
1804  ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * 3*4;
1805  ep->ext_buf = (uint8_t *)bu_malloc(ep->ext_nbytes, "ehy external");
1806 
1807  if (!NEAR_EQUAL(MAGNITUDE(xip->ehy_Au), 1.0, RT_LEN_TOL)) {
1808  bu_log("rt_ehy_export4: Au not a unit vector!\n");
1809  return -1;
1810  }
1811 
1812  if (MAGNITUDE(xip->ehy_H) < RT_LEN_TOL
1813  || xip->ehy_c < RT_LEN_TOL
1814  || xip->ehy_r1 < RT_LEN_TOL
1815  || xip->ehy_r2 < RT_LEN_TOL) {
1816  bu_log("rt_ehy_export4: not all dimensions positive!\n");
1817  return -1;
1818  }
1819 
1820  if (!NEAR_ZERO(VDOT(xip->ehy_Au, xip->ehy_H), RT_DOT_TOL)) {
1821  bu_log("rt_ehy_export4: Au and H are not perpendicular!\n");
1822  return -1;
1823  }
1824 
1825  if (xip->ehy_r2 > xip->ehy_r1) {
1826  bu_log("rt_ehy_export4: semi-minor axis cannot be longer than semi-major axis!\n");
1827  return -1;
1828  }
1829 
1830  /* Warning: type conversion */
1831  VSCALE(&vec[0*3], xip->ehy_V, local2mm);
1832  VSCALE(&vec[1*3], xip->ehy_H, local2mm);
1833  /* don't scale ehy_Au (unit vector!!) */
1834  VMOVE(&vec[2*3], xip->ehy_Au);
1835  vec[3*3] = xip->ehy_r1 * local2mm;
1836  vec[3*3+1] = xip->ehy_r2 * local2mm;
1837  vec[3*3+2] = xip->ehy_c * local2mm;
1838 
1839  /* Convert from internal (host) to database (network) format */
1840  bu_cv_htond(ep->ext_buf, (unsigned char *)vec, 3*4);
1841 
1842  return 0;
1843 }
1844 
1845 
1846 /**
1847  * Make human-readable formatted presentation of this solid. First
1848  * line describes type of solid. Additional lines are indented one
1849  * tab, and give parameter values.
1850  */
1851 int
1852 rt_ehy_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
1853 {
1854  struct rt_ehy_internal *xip = (struct rt_ehy_internal *)ip->idb_ptr;
1855  char buf[256];
1856 
1857  RT_EHY_CK_MAGIC(xip);
1858  bu_vls_strcat(str, "Elliptical Hyperboloid (EHY)\n");
1859 
1860  sprintf(buf, "\tV (%g, %g, %g)\n",
1861  INTCLAMP(xip->ehy_V[X] * mm2local),
1862  INTCLAMP(xip->ehy_V[Y] * mm2local),
1863  INTCLAMP(xip->ehy_V[Z] * mm2local));
1864  bu_vls_strcat(str, buf);
1865 
1866  sprintf(buf, "\tH (%g, %g, %g) mag=%g\n",
1867  INTCLAMP(xip->ehy_H[X] * mm2local),
1868  INTCLAMP(xip->ehy_H[Y] * mm2local),
1869  INTCLAMP(xip->ehy_H[Z] * mm2local),
1870  INTCLAMP(MAGNITUDE(xip->ehy_H) * mm2local));
1871  bu_vls_strcat(str, buf);
1872 
1873  if (!verbose)
1874  return 0;
1875 
1876  sprintf(buf, "\tA=%g\n", INTCLAMP(xip->ehy_r1 * mm2local));
1877  bu_vls_strcat(str, buf);
1878 
1879  sprintf(buf, "\tB=%g\n", INTCLAMP(xip->ehy_r2 * mm2local));
1880  bu_vls_strcat(str, buf);
1881 
1882  sprintf(buf, "\tc=%g\n", INTCLAMP(xip->ehy_c * mm2local));
1883  bu_vls_strcat(str, buf);
1884 
1885  return 0;
1886 }
1887 
1888 
1889 /**
1890  * Free the storage associated with the rt_db_internal version of this
1891  * solid.
1892  */
1893 void
1895 {
1896  struct rt_ehy_internal *xip;
1897 
1898  RT_CK_DB_INTERNAL(ip);
1899 
1900  xip = (struct rt_ehy_internal *)ip->idb_ptr;
1901  RT_EHY_CK_MAGIC(xip);
1902  xip->ehy_magic = 0; /* sanity */
1903 
1904  bu_free((char *)xip, "ehy ifree");
1905  ip->idb_ptr = ((void *)0); /* sanity */
1906 }
1907 
1908 
1909 int
1910 rt_ehy_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
1911 {
1912  if (!ps) return 0;
1913  if (ip) RT_CK_DB_INTERNAL(ip);
1914 
1915  return 0; /* OK */
1916 }
1917 
1918 
1919 static int
1920 ehy_is_valid(struct rt_ehy_internal *ehy)
1921 {
1922  fastf_t mag_h, cos_angle_ah;
1923  vect_t a, h;
1924 
1925  RT_EHY_CK_MAGIC(ehy);
1926 
1927  if (!(ehy->ehy_r1 > 0.0 && ehy->ehy_r2 > 0.0 && ehy->ehy_c > 0.0)) {
1928  return 0;
1929  }
1930 
1931  VMOVE(h, ehy->ehy_H);
1932  VMOVE(a, ehy->ehy_Au);
1933 
1934  /* Check that A is a unit vector. If it is, then it should be true that
1935  * |A| == |A|^2 == 1.0.
1936  */
1937  if (!NEAR_EQUAL(MAGSQ(a), 1.0, RT_LEN_TOL)) {
1938  return 0;
1939  }
1940 
1941  /* check that |H| > 0.0 */
1942  mag_h = MAGNITUDE(h);
1943  if (NEAR_ZERO(mag_h, RT_LEN_TOL)) {
1944  return 0;
1945  }
1946 
1947  /* check that A and H are orthogonal */
1948  cos_angle_ah = VDOT(a, h) / mag_h;
1949  if (!NEAR_ZERO(cos_angle_ah, RT_DOT_TOL)) {
1950  return 0;
1951  }
1952 
1953  return 1;
1954 }
1955 
1956 
1957 void
1958 rt_ehy_surf_area(fastf_t *area, const struct rt_db_internal *ip)
1959 {
1960  struct rt_ehy_internal *eip;
1961  fastf_t a, b, h, integralArea, sqrt_rb;
1962  RT_CK_DB_INTERNAL(ip);
1963  eip = (struct rt_ehy_internal *)ip->idb_ptr;
1964  RT_EHY_CK_MAGIC(eip);
1965 
1966  a = eip->ehy_c;
1967  h = MAGNITUDE(eip->ehy_H);
1968  b = (eip->ehy_r1 * a) / sqrt(h * (h - 2 * a));
1969 
1970  /** Formula taken from : https://docs.google.com/file/d/0BydeQ6BPlVejRWt6NlJLVDl0d28/edit
1971  * Area can be calculated by substracting integral of hyperbola from the area of the bounding rectangle
1972  */
1973  sqrt_rb = sqrt(eip->ehy_r1 * eip->ehy_r1 + b * b);
1974  integralArea = (a / b) * ((eip->ehy_r1 * sqrt_rb) + ((b * b / 2) * (log(sqrt_rb + eip->ehy_r1) - log(sqrt_rb - eip->ehy_r1))));
1975  *area = 2 * eip->ehy_r1 * (a + h) - integralArea;
1976 }
1977 
1978 
1979 /** @} */
1980 /*
1981  * Local Variables:
1982  * mode: C
1983  * tab-width: 8
1984  * indent-tabs-mode: t
1985  * c-file-style: "stroustrup"
1986  * End:
1987  * ex: shiftwidth=4 tabstop=8
1988  */
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
int rt_ehy_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
Definition: ehy.c:1160
Definition: raytrace.h:800
#define NMG_EDGEUSE_MAGIC
Definition: magic.h:120
#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
fastf_t ehy_cprime
Definition: ehy.c:174
#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)
Definition: list.h:118
int nmg_fu_planeeqn(struct faceuse *fu, const struct bn_tol *tol)
Definition: nmg_mod.c:1311
#define RT_DOT_TOL
Definition: raytrace.h:170
#define RT_CK_APPLICATION(_p)
Definition: raytrace.h:1675
int rt_ehy_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: ehy.c:1608
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
void bu_vls_strcat(struct bu_vls *vp, const char *s)
Definition: vls.c:368
lu
Definition: nmg_mod.c:3855
#define VSET(a, b, c, d)
Definition: color.c:53
#define VSETALL(a, s)
Definition: color.c:54
Definition: raytrace.h:215
void rt_ehy_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
Definition: ehy.c:488
#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
int rt_ehy_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
Definition: ehy.c:257
#define SMALL_FASTF
Definition: defines.h:342
int rt_ehy_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *tol)
Definition: ehy.c:195
fastf_t st_aradius
Radius of APPROXIMATING sphere.
Definition: raytrace.h:433
void rt_ehy_print(const struct soltab *stp)
Definition: ehy.c:333
fastf_t ell_angle(fastf_t *p1, fastf_t a, fastf_t b, fastf_t dtol, fastf_t ntol)
Definition: ell.c:1788
void rt_ehy_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
Definition: ehy.c:524
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 bn_eigen2x2(fastf_t *val1, fastf_t *val2, vect_t vec1, vect_t vec2, fastf_t a, fastf_t b, fastf_t c)
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_ehy_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
Definition: ehy.c:1741
vect_t ehy_Hunit
Definition: ehy.c:169
int rt_ehy_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: ehy.c:1789
NMG_CK_LOOPUSE(lu)
struct bu_list l
Definition: raytrace.h:369
fastf_t curve_spacing
Definition: raytrace.h:1940
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
Definition: ptbl.h:62
if(share_geom)
Definition: nmg_mod.c:3829
int idb_major_type
Definition: raytrace.h:192
Definition: color.c:49
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
#define ID_EHY
Elliptical Hyperboloid.
Definition: raytrace.h:478
void bn_mat_print(const char *title, const mat_t m)
Definition: mat.c:81
void rt_ehy_ifree(struct rt_db_internal *ip)
Definition: ehy.c:1894
void nmg_shell_coplanar_face_merge(struct shell *s, const struct bn_tol *tol, const int simplify)
Definition: nmg_mod.c:93
#define RT_CK_DB_INTERNAL(_p)
Definition: raytrace.h:207
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
fastf_t st_bradius
Radius of BOUNDING sphere.
Definition: raytrace.h:434
fastf_t crv_c2
curvature in other direction
Definition: raytrace.h:309
#define RT_EHY_INTERNAL_MAGIC
Definition: magic.h:90
int rt_ehy_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: ehy.c:914
#define BU_PTBL_GET(ptbl, i)
Definition: ptbl.h:108
#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)
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
uint8_t * ext_buf
Definition: parse.h:216
static void top()
#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
int rt_ehy_adaptive_plot(struct rt_db_internal *ip, const struct rt_view_info *info)
Definition: ehy.c:836
int rt_ehy_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
Definition: ehy.c:1852
struct hit seg_out
OUT information.
Definition: raytrace.h:371
#define EHY_NORM_TOP
Definition: ehy.c:350
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
void rt_ehy_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
Definition: ehy.c:583
struct bu_list l
Definition: ptbl.h:63
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
#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 BU_LIST_FIRST_MAGIC(hp)
Definition: list.h:416
vect_t ehy_Bunit
Definition: ehy.c:171
mat_t ehy_invRoS
Definition: ehy.c:173
uint32_t magic
Magic # for mem id/check.
Definition: list.h:119
#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
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
void bu_ptbl_free(struct bu_ptbl *b)
Definition: ptbl.c:226
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
struct rt_pt_node * next
ptr to next pt
Definition: raytrace.h:1912
#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
void bu_cv_ntohd(unsigned char *out, const unsigned char *in, size_t count)
int rt_ehy_params(struct pc_pc_set *ps, const struct rt_db_internal *ip)
Definition: ehy.c:1910
const struct rt_functab OBJ[]
Definition: table.c:159
#define BU_PTBL_END(ptbl)
Definition: ptbl.h:106
#define BU_ASSERT_PTR(_lhs, _relation, _rhs)
Definition: defines.h:227
void bn_vec_ortho(vect_t out, const vect_t in)
#define R
Definition: msr.c:55
int rt_ehy_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
Definition: ehy.c:362
const struct bu_structparse rt_ehy_parse[]
Definition: ehy.c:178
void rt_ell(fastf_t *ov, const fastf_t *V, const fastf_t *A, const fastf_t *B, int sides)
Definition: epa.c:1141
void * st_specific
-> ID-specific (private) struct
Definition: raytrace.h:435
fastf_t uv_du
delta in u
Definition: raytrace.h:343
void rt_ehy_free(struct soltab *stp)
Definition: ehy.c:630
point_t p
a point
Definition: raytrace.h:1911
mat_t ehy_SoR
Definition: ehy.c:172
fastf_t crv_c1
curvature in principle dir
Definition: raytrace.h:308
#define A
Definition: msr.c:51
Definition: color.c:51
int dbi_version
PRIVATE: use db_version()
Definition: raytrace.h:824
void plot_ellipse(struct bu_list *vhead, const vect_t t, const vect_t a, const vect_t b, int num_points)
void rt_ehy_surf_area(fastf_t *area, const struct rt_db_internal *ip)
Definition: ehy.c:1958
double norm
normal tol
Definition: raytrace.h:182
int hit_surfno
solid-specific surface indicator
Definition: raytrace.h:255
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
#define RT_PCOEF_TOL
Definition: raytrace.h:171
int st_id
Solid ident.
Definition: raytrace.h:431
size_t ext_nbytes
Definition: parse.h:210
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
#define EHY_NORM_BODY
Definition: ehy.c:349
Definition: vls.h:56
double fastf_t
Definition: defines.h:300
#define VPRINT(a, b)
Definition: raytrace.h:1881
point_t ehy_V
Definition: ehy.c:168
int rt_ehy_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
Definition: ehy.c:1680
Definition: color.c:50
#define BU_LIST_FIRST(structure, hp)
Definition: list.h:312
fastf_t uv_v
Range 0..1.
Definition: raytrace.h:342
void nmg_vertex_tabulate(struct bu_ptbl *tab, const uint32_t *magic_p)
Definition: nmg_info.c:1985
point_t st_center
Centroid of solid.
Definition: raytrace.h:432
vect_t ehy_Aunit
Definition: ehy.c:170
#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