g_eto.c

Go to the documentation of this file.
00001 /*                         G _ E T O . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 1992-2006 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * as published by the Free Software Foundation; either version 2 of
00010  * the License, or (at your option) any later version.
00011  *
00012  * This library is distributed in the hope that it will be useful, but
00013  * WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015  * Library General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU Lesser General Public
00018  * License along with this file; see the file named COPYING for more
00019  * information.
00020  */
00021 
00022 /** @addtogroup g_  */
00023 /*@{*/
00024 /** @file g_eto.c
00025  *      Intersect a ray with an Elliptical Torus.
00026  *
00027  * Authors -
00028  *      Michael Markowski       (Programming)
00029  *      ERIM GIFT code          (ETO Eqn)
00030  *
00031  *  Source -
00032  *      SECAD/VLD Computing Consortium, Bldg 394
00033  *      The U. S. Army Ballistic Research Laboratory
00034  *      Aberdeen Proving Ground, Maryland  21005
00035  *
00036  */
00037 
00038 #ifndef lint
00039 static const char RCSeto[] = "@(#)$Header: /cvsroot/brlcad/brlcad/src/librt/g_eto.c,v 14.14 2006/09/16 02:04:24 lbutler Exp $ (BRL)";
00040 #endif
00041 
00042 #include "common.h"
00043 
00044 #include <stddef.h>
00045 #include <stdio.h>
00046 #ifdef HAVE_STRING_H
00047 #  include <string.h>
00048 #else
00049 #  include <strings.h>
00050 #endif
00051 #include <math.h>
00052 
00053 #include "machine.h"
00054 #include "vmath.h"
00055 #include "db.h"
00056 #include "nmg.h"
00057 #include "raytrace.h"
00058 #include "rtgeom.h"
00059 #include "./debug.h"
00060 
00061 /*
00062  * The ETO has the following input fields:
00063  *      V       V from origin to center.
00064  *      N       Normal to plane of eto.
00065  *      C       Semi-major axis of elliptical cross section.
00066  *      r       Radius of revolution.
00067  *      rd      Semi-minor axis of elliptical cross section.
00068  *
00069  */
00070 
00071 /*
00072  *  Algorithm:
00073  *
00074  *  Given V, N, C, r, and rd, there is a set of points on this eto
00075  *
00076  *  { (x,y,z) | (x,y,z) is on eto defined by V, N, C, r, rd }
00077  *
00078  *  Through a series of  Transformations, this set will be
00079  *  transformed into a set of points on an eto centered at the origin
00080  *  which lies on the X-Y plane (ie, N is on the Z axis).
00081  *
00082  *  { (x',y',z') | (x',y',z') is an eto at origin }
00083  *
00084  *  The transformation from X to X' is accomplished by:
00085  *
00086  *  X' = R( X - V )
00087  *
00088  *  where R(X) =  ( B/(|B|) )
00089  *               (  A/(|A|)  ) . X
00090  *                ( N/(|N|) )
00091  *
00092  *  To find the intersection of a line with the eto, consider
00093  *  the parametric line L:
00094  *
00095  *      L : { P(n) | P + t(n) . D }
00096  *
00097  *  Call W the actual point of intersection between L and the eto.
00098  *  Let W' be the point of intersection between L' and the unit eto.
00099  *
00100  *      L' : { P'(n) | P' + t(n) . D' }
00101  *
00102  *  W = invR( W' ) + V
00103  *
00104  *  Where W' = k D' + P'.
00105  *
00106  *
00107  *  Given a line and a ratio, alpha, finds the equation of the
00108  *  unit eto in terms of the variable 't'.
00109  *
00110  *  Given that the equation for the eto is:
00111  *
00112  *            _______                           _______
00113  *           / 2    2              2           / 2    2               2
00114  *  [Eu(+- \/ x  + y   - R) + Ev z]  + [Fu(+-\/ x  + y   - R) + Fv z ]
00115  * --------------------------------    -------------------------------  = 1
00116  *               2                                      2
00117  *             Rc                                     Rd
00118  *
00119  *  First, find X, Y, and Z in terms of 't' for this line, then
00120  *  substitute them into the equation above.
00121  *
00122  *      Wx = Dx*t + Px
00123  *
00124  *      Wx**2 = Dx**2 * t**2  +  2 * Dx * Px  +  Px**2
00125  *
00126  *  The real roots of the equation in 't' are the intersect points
00127  *  along the parameteric line.
00128  *
00129  *  NORMALS.  Given the point W on the eto, what is the vector
00130  *  normal to the tangent plane at that point?
00131  *
00132  *  Map W onto the eto, ie:  W' = R( W - V ).
00133  *  In this case, we find W' by solving the parameteric line given k.
00134  *
00135  *  The gradient of the eto at W' is in fact the
00136  *  normal vector.
00137  *
00138  *  For f(x,y,z) = 0, the gradient of f() is ( df/dx, df/dy, df/dz ).
00139  *
00140  *  Note that the normal vector (gradient) produced above will not have
00141  *  unit length.  Also, to make this useful for the original eto, it will
00142  *  have to be rotated back to the orientation of the original eto.
00143  */
00144 
00145 struct eto_specific {
00146         vect_t  eto_V;          /* Vector to center of eto */
00147         vect_t  eto_N;          /* unit normal to plane of eto */
00148         vect_t  eto_C;          /* semi-major axis of ellipse */
00149         fastf_t eto_r;          /* radius of revolution */
00150         fastf_t eto_rc;         /* semi-major axis of ellipse */
00151         fastf_t eto_rd;         /* semi-minor axis of ellipse */
00152         mat_t   eto_R;          /* Rot(vect) */
00153         mat_t   eto_invR;       /* invRot(vect') */
00154         fastf_t eu, ev, fu, fv;
00155 };
00156 
00157 const struct bu_structparse rt_eto_parse[] = {
00158     { "%f", 3, "V",   bu_offsetof(struct rt_eto_internal, eto_V[X]), BU_STRUCTPARSE_FUNC_NULL },
00159     { "%f", 3, "N",   bu_offsetof(struct rt_eto_internal, eto_N[X]), BU_STRUCTPARSE_FUNC_NULL },
00160     { "%f", 3, "C",   bu_offsetof(struct rt_eto_internal, eto_C[X]), BU_STRUCTPARSE_FUNC_NULL },
00161     { "%f", 1, "r",   bu_offsetof(struct rt_eto_internal, eto_r),    BU_STRUCTPARSE_FUNC_NULL },
00162     { "%f", 1, "r_d", bu_offsetof(struct rt_eto_internal, eto_rd),   BU_STRUCTPARSE_FUNC_NULL },
00163     { {'\0','\0','\0','\0'}, 0, (char *)NULL, 0, BU_STRUCTPARSE_FUNC_NULL }
00164  };
00165 
00166 /**
00167  *                      R T _ E T O _ P R E P
00168  *
00169  *  Given a pointer to a GED database record, and a transformation matrix,
00170  *  determine if this is a valid eto, and if so, precompute various
00171  *  terms of the formula.
00172  *
00173  *  Returns -
00174  *      0       ETO is OK
00175  *      !0      Error in description
00176  *
00177  *  Implicit return -
00178  *      A struct eto_specific is created, and it's address is stored in
00179  *      stp->st_specific for use by rt_eto_shot().
00180  */
00181 int
00182 rt_eto_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
00183 {
00184         register struct eto_specific *eto;
00185         LOCAL vect_t    P, w1;  /* for RPP calculation */
00186         LOCAL vect_t    Au, Bu, Cu, Nu;
00187         FAST fastf_t    ch, cv, dh, f, phi;
00188         struct rt_eto_internal  *tip;
00189 
00190         tip = (struct rt_eto_internal *)ip->idb_ptr;
00191         RT_ETO_CK_MAGIC(tip);
00192 
00193         /* Solid is OK, compute constant terms now */
00194         BU_GETSTRUCT( eto, eto_specific );
00195         stp->st_specific = (genptr_t)eto;
00196 
00197         eto->eto_r = tip->eto_r;
00198         eto->eto_rd = tip->eto_rd;
00199         eto->eto_rc = MAGNITUDE( tip->eto_C );
00200         if ( NEAR_ZERO(eto->eto_r, 0.0001) || NEAR_ZERO(eto->eto_rd, 0.0001)
00201                 || NEAR_ZERO(eto->eto_rc, 0.0001)) {
00202                 bu_log("eto(%s): r, rd, or rc zero length\n", stp->st_name);
00203                 return(1);
00204         }
00205 
00206         VMOVE( eto->eto_V, tip->eto_V );
00207         VMOVE( eto->eto_N, tip->eto_N );
00208         VMOVE( eto->eto_C, tip->eto_C );
00209         VMOVE( Nu, tip->eto_N );
00210         VUNITIZE( Nu );         /* z axis of coord sys */
00211         bn_vec_ortho( Bu, Nu ); /* x axis */
00212         VUNITIZE( Bu );
00213         VCROSS( Au, Nu, Bu );   /* y axis */
00214         VMOVE( Cu, tip->eto_C );
00215         VUNITIZE( Cu );
00216 
00217         /* get horizontal and vertical components of C and Rd */
00218         cv = VDOT( eto->eto_C, Nu );
00219         ch = sqrt( VDOT( eto->eto_C, eto->eto_C ) - cv * cv );
00220         /* angle between C and Nu */
00221         phi = acos( cv / eto->eto_rc );
00222         dh = eto->eto_rd * cos(phi);
00223         /* make sure ellipse doesn't overlap itself when revolved */
00224         if (ch > eto->eto_r || dh > eto->eto_r) {
00225                 bu_log("eto(%s): revolved ellipse overlaps itself\n",
00226                         stp->st_name);
00227                 return(1);
00228         }
00229 
00230         eto->ev = fabs( VDOT( Cu, Nu ) );       /* vertical component of Cu */
00231         eto->eu = sqrt( 1.0 - eto->ev * eto->ev );      /* horiz component */
00232         eto->fu = -eto->ev;
00233         eto->fv =  eto->eu;
00234 
00235         /* Compute R and invR matrices */
00236         MAT_IDN( eto->eto_R );
00237         VMOVE( &eto->eto_R[0], Bu );
00238         VMOVE( &eto->eto_R[4], Au );
00239         VMOVE( &eto->eto_R[8], Nu );
00240         bn_mat_inv( eto->eto_invR, eto->eto_R );
00241 
00242         stp->st_aradius = stp->st_bradius = eto->eto_r + eto->eto_rc;
00243 
00244         /*
00245          *  Compute the bounding RPP planes for a circular eto.
00246          *
00247          *  Given a circular eto with vertex V, vector N, and
00248          *  radii r and |C|.  A bounding plane with direction
00249          *  vector P will touch the surface of the eto at the
00250          *  points:  V +/- [|C| + r * |N x P|] P
00251          */
00252 
00253         /* X */
00254         VSET( P, 1.0, 0, 0 );           /* bounding plane normal */
00255         VCROSS( w1, Nu, P );    /* for sin(angle N P) */
00256         f = eto->eto_rc + eto->eto_r * MAGNITUDE(w1);
00257         VSCALE( w1, P, f );
00258         f = fabs( w1[X] );
00259         stp->st_min[X] = eto->eto_V[X] - f;
00260         stp->st_max[X] = eto->eto_V[X] + f;
00261 
00262         /* Y */
00263         VSET( P, 0, 1.0, 0 );           /* bounding plane normal */
00264         VCROSS( w1, Nu, P );    /* for sin(angle N P) */
00265         f = eto->eto_rc + eto->eto_r * MAGNITUDE(w1);
00266         VSCALE( w1, P, f );
00267         f = fabs( w1[Y] );
00268         stp->st_min[Y] = eto->eto_V[Y] - f;
00269         stp->st_max[Y] = eto->eto_V[Y] + f;
00270 
00271         /* Z */
00272         VSET( P, 0, 0, 1.0 );           /* bounding plane normal */
00273         VCROSS( w1, Nu, P );    /* for sin(angle N P) */
00274         f = eto->eto_rc + eto->eto_r * MAGNITUDE(w1);
00275         VSCALE( w1, P, f );
00276         f = fabs( w1[Z] );
00277         stp->st_min[Z] = eto->eto_V[Z] - f;
00278         stp->st_max[Z] = eto->eto_V[Z] + f;
00279 
00280         return(0);                      /* OK */
00281 }
00282 
00283 /**
00284  *                      R T _ E T O _ P R I N T
00285  */
00286 void
00287 rt_eto_print(register const struct soltab *stp)
00288 {
00289         register const struct eto_specific *eto =
00290                 (struct eto_specific *)stp->st_specific;
00291 
00292         VPRINT("V", eto->eto_V);
00293         VPRINT("N", eto->eto_N);
00294         VPRINT("C", eto->eto_C);
00295         bu_log("r = %f\n", eto->eto_r);
00296         bu_log("rc = %f\n", eto->eto_rc);
00297         bu_log("rd = %f\n", eto->eto_rd);
00298         bn_mat_print("R", eto->eto_R );
00299         bn_mat_print("invR", eto->eto_invR );
00300         bu_log( "rpp: (%g, %g, %g) to (%g, %g, %g)\n",
00301                 stp->st_min[X], stp->st_min[Y], stp->st_min[Z],
00302                 stp->st_max[X], stp->st_max[Y], stp->st_max[Z]);
00303 }
00304 
00305 /**
00306  *                      R T _ E T O _ S H O T
00307  *
00308  *  Intersect a ray with an eto, where all constant terms have
00309  *  been precomputed by rt_eto_prep().  If an intersection occurs,
00310  *  one or two struct seg(s) will be acquired and filled in.
00311  *
00312  *  NOTE:  All lines in this function are represented parametrically
00313  *  by a point,  P( x0, y0, z0 ) and a direction normal,
00314  *  D = ax + by + cz.  Any point on a line can be expressed
00315  *  by one variable 't', where
00316  *
00317  *      X = a*t + x0,   eg,  X = Dx*t + Px
00318  *      Y = b*t + y0,
00319  *      Z = c*t + z0.
00320  *
00321  *  First, convert the line to the coordinate system of a "stan-
00322  *  dard" eto.  This is a eto which lies in the X-Y plane
00323  *  and circles the origin.  The semimajor axis is C.
00324  *
00325  *  Then find the equation of that line and the standard eto,
00326  *  which turns out to be a quartic equation in 't'.  Solve the
00327  *  equation using a general polynomial root finder.  Use those
00328  *  values of 't' to compute the points of intersection in the
00329  *  original coordinate system.
00330  *
00331  *  Returns -
00332  *      0       MISS
00333  *      >0      HIT
00334  */
00335 int
00336 rt_eto_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
00337 {
00338         register struct eto_specific *eto =
00339                 (struct eto_specific *)stp->st_specific;
00340         register struct seg *segp;
00341         LOCAL vect_t    dprime;         /* D' */
00342         LOCAL vect_t    pprime;         /* P' */
00343         LOCAL vect_t    work;           /* temporary vector */
00344         LOCAL bn_poly_t C;              /* The final equation */
00345         LOCAL bn_complex_t      val[4]; /* The complex roots */
00346         LOCAL double    k[4];           /* The real roots */
00347         register int    i;
00348         LOCAL int       j;
00349         LOCAL vect_t    cor_pprime;     /* new ray origin */
00350         LOCAL fastf_t   cor_proj;
00351         LOCAL fastf_t   A1,A2,A3,A4,A5,A6,A7,A8,B1,B2,B3,C1,C2,C3,D1,term;
00352 
00353         /* Convert vector into the space of the unit eto */
00354         MAT4X3VEC( dprime, eto->eto_R, rp->r_dir );
00355         VUNITIZE( dprime );
00356 
00357         VSUB2( work, rp->r_pt, eto->eto_V );
00358         MAT4X3VEC( pprime, eto->eto_R, work );
00359 
00360         /* normalize distance from eto.  substitute
00361          * corrected pprime which contains a translation along ray
00362          * direction to closest approach to vertex of eto.
00363          * Translating ray origin along direction of ray to closest pt. to
00364          * origin of solid's coordinate system, new ray origin is
00365          * 'cor_pprime'.
00366          */
00367         cor_proj = VDOT( pprime, dprime );
00368         VSCALE( cor_pprime, dprime, cor_proj );
00369         VSUB2( cor_pprime, pprime, cor_pprime );
00370 
00371         /*
00372          *  NOTE: The following code is based on code in
00373          *        eto.f by ERIM for GIFT.
00374          *
00375          *  Given a line, finds the equation of the
00376          *  eto in terms of the variable 't'.
00377          *
00378          *  The equation for the eto is:
00379          *
00380             _______                           ________
00381            / 2    2              2           / 2    2               2
00382   [Eu(+- \/ x  + y   - R) + Ev z]  + [Fu(+-\/ x  + y   - R) + Fv z ]
00383  --------------------------------    -------------------------------  = 1
00384                2                                      2
00385              Rc                                     Rd
00386          *
00387          *                  ^   ^
00388          *       where Ev = C . N
00389          *
00390          *                  ^   ^
00391          *             Eu = C . A
00392          *
00393          *                  ^   ^
00394          *             Fv = C . A
00395          *
00396          *                  ^   ^
00397          *             Fu =-C . N.
00398          *
00399          *  First, find X, Y, and Z in terms of 't' for this line, then
00400          *  substitute them into the equation above.
00401          *
00402          *      Wx = Dx*t + Px, etc.
00403          *
00404          *  Regrouping coefficients and constants, the equation can then
00405          *  be rewritten as:
00406          *
00407          *  [A1*sqrt(C1 + C2*t + C3*t^2) + A2 + A3*t]^2  +
00408          *  [B1*sqrt(C1 + C2*t + C3*t^2) + B2 + B3*t]^2  - D1 = 0
00409          *
00410          *  where, (variables defined in code below)
00411          */
00412         A1 = eto->eto_rd * eto->eu;
00413         B1 = eto->eto_rc * eto->fu;
00414         C1 = cor_pprime[X] * cor_pprime[X] + cor_pprime[Y] * cor_pprime[Y];
00415         C2 = 2 * (dprime[X] * cor_pprime[X] + dprime[Y] * cor_pprime[Y]);
00416         C3 = dprime[X] * dprime[X] + dprime[Y] * dprime[Y];
00417         A2 = -eto->eto_rd * eto->eto_r * eto->eu + eto->eto_rd * eto->ev * cor_pprime[Z];
00418         B2 = -eto->eto_rc * eto->eto_r * eto->fu + eto->eto_rc * eto->fv * cor_pprime[Z];
00419         A3 = eto->eto_rd * eto->ev * dprime[Z];
00420         B3 = eto->eto_rc * eto->fv * dprime[Z];
00421         D1 = eto->eto_rc * eto->eto_rc * eto->eto_rd * eto->eto_rd;
00422 
00423         /*
00424          *  Squaring the two terms above and again regrouping coefficients
00425          *  the equation now becomes:
00426          *
00427          *  A6*t^2 + A5*t + A4 = -(A8*t + A7)*sqrt(C1 + C2*t + C3*t^2)
00428          *
00429          *  where, (variables defined in code)
00430          */
00431         A4 = A1*A1*C1 + B1*B1*C1 + A2*A2 + B2*B2 - D1;
00432         A5 = A1*A1*C2 + B1*B1*C2 + 2*A2*A3 + 2*B2*B3;
00433         A6 = A1*A1*C3 + B1*B1*C3 + A3*A3 + B3*B3;
00434         A7 = 2*(A1*A2 + B1*B2);
00435         A8 = 2*(A1*A3 + B1*B3);
00436         term = A6*A6 - A8*A8*C3;
00437 
00438         /*
00439          *  Squaring both sides and subtracting RHS from LHS yields:
00440          */
00441         C.dgr=4;
00442         C.cf[4] = (A4*A4 - A7*A7*C1);                   /* t^0 */
00443         C.cf[3] = (2*A4*A5 - A7*A7*C2 - 2*A7*A8*C1);    /* t^1 */
00444         C.cf[2] = (2*A4*A6 + A5*A5 - A7*A7*C3 - 2*A7*A8*C2 - A8*A8*C1); /* t^2 */
00445         C.cf[1] = (2*A5*A6 - 2*A7*A8*C3 - A8*A8*C2);    /* t^3 */
00446         C.cf[0] = term;                                 /* t^4 */
00447         /* NOTE: End of ERIM based code */
00448 
00449         /*  It is known that the equation is 4th order.  Therefore,
00450          *  if the root finder returns other than 4 roots, error.
00451          */
00452         if ( (i = rt_poly_roots( &C, val, stp->st_dp->d_namep )) != 4 ){
00453                 if( i > 0 )  {
00454                         bu_log("eto:  rt_poly_roots() 4!=%d\n", i);
00455                         bn_pr_roots( stp->st_name, val, i );
00456                 } else if (i < 0) {
00457                     static int reported=0;
00458                     bu_log("The root solver failed to converge on a solution for %s\n", stp->st_dp->d_namep);
00459                     if (!reported) {
00460                         VPRINT("while shooting from:\t", rp->r_pt);
00461                         VPRINT("while shooting at:\t", rp->r_dir);
00462                         bu_log("Additional elliptical torus convergence failure details will be suppressed.\n");
00463                         reported=1;
00464                     }
00465                 }
00466                 return(0);              /* MISS */
00467         }
00468 
00469         /*  Only real roots indicate an intersection in real space.
00470          *
00471          *  Look at each root returned; if the imaginary part is zero
00472          *  or sufficiently close, then use the real part as one value
00473          *  of 't' for the intersections
00474          */
00475         for ( j=0, i=0; j < 4; j++ ){
00476                 if( NEAR_ZERO( val[j].im, 0.0001 ) )
00477                         k[i++] = val[j].re;
00478         }
00479 
00480         /* reverse above translation by adding distance to all 'k' values. */
00481         for( j = 0; j < i; ++j )
00482                 k[j] -= cor_proj;
00483 
00484         /* Here, 'i' is number of points found */
00485         switch( i )  {
00486         case 0:
00487                 return(0);              /* No hit */
00488 
00489         default:
00490                 bu_log("rt_eto_shot: reduced 4 to %d roots\n",i);
00491                 bn_pr_roots( stp->st_name, val, 4 );
00492                 return(0);              /* No hit */
00493 
00494         case 2:
00495                 {
00496                         /* Sort most distant to least distant. */
00497                         FAST fastf_t    u;
00498                         if( (u=k[0]) < k[1] )  {
00499                                 /* bubble larger towards [0] */
00500                                 k[0] = k[1];
00501                                 k[1] = u;
00502                         }
00503                 }
00504                 break;
00505         case 4:
00506                 {
00507                         register short  n;
00508                         register short  lim;
00509 
00510                         /*  Inline rt_pt_sort().  Sorts k[] into descending order. */
00511                         for( lim = i-1; lim > 0; lim-- )  {
00512                                 for( n = 0; n < lim; n++ )  {
00513                                         FAST fastf_t    u;
00514                                         if( (u=k[n]) < k[n+1] )  {
00515                                                 /* bubble larger towards [0] */
00516                                                 k[n] = k[n+1];
00517                                                 k[n+1] = u;
00518                                         }
00519                                 }
00520                         }
00521                 }
00522                 break;
00523         }
00524 
00525         /* Now, t[0] > t[npts-1] */
00526         /* k[1] is entry point, and k[0] is farthest exit point */
00527         RT_GET_SEG(segp, ap->a_resource);
00528         segp->seg_stp = stp;
00529         segp->seg_in.hit_dist = k[1];
00530         segp->seg_out.hit_dist = k[0];
00531         /* Set aside vector for rt_eto_norm() later */
00532         VJOIN1( segp->seg_in.hit_vpriv, pprime, k[1], dprime );
00533         VJOIN1( segp->seg_out.hit_vpriv, pprime, k[0], dprime );
00534         BU_LIST_INSERT( &(seghead->l), &(segp->l) );
00535 
00536         if( i == 2 )
00537                 return(2);                      /* HIT */
00538 
00539         /* 4 points */
00540         /* k[3] is entry point, and k[2] is exit point */
00541         RT_GET_SEG(segp, ap->a_resource);
00542         segp->seg_stp = stp;
00543         segp->seg_in.hit_dist = k[3];
00544         segp->seg_out.hit_dist = k[2];
00545         VJOIN1( segp->seg_in.hit_vpriv, pprime, k[3], dprime );
00546         VJOIN1( segp->seg_out.hit_vpriv, pprime, k[2], dprime );
00547         BU_LIST_INSERT( &(seghead->l), &(segp->l) );
00548         return(4);                      /* HIT */
00549 }
00550 
00551 #define SEG_MISS(SEG)           (SEG).seg_stp=(struct soltab *) 0;
00552 /**
00553  *                      R T _ E T O _ V S H O T
00554  *
00555  *  This is the Becker vector version
00556  */
00557 void
00558 rt_eto_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
00559                                /* An array of solid pointers */
00560                                /* An array of ray pointers */
00561                                /* array of segs (results returned) */
00562                                /* Number of ray/object pairs */
00563 
00564 {
00565 }
00566 
00567 /**
00568  *                      R T _ E T O _ N O R M
00569  *
00570  *  Compute the normal to the eto,
00571  *  given a point on the eto centered at the origin on the X-Y plane.
00572  *  The gradient of the eto at that point is in fact the
00573  *  normal vector, which will have to be given unit length.
00574  *  To make this useful for the original eto, it will have
00575  *  to be rotated back to the orientation of the original eto.
00576  *  The equation for the eto is:
00577  *
00578  *            _______                           ________
00579  *           / 2    2              2           / 2    2               2
00580  *  [Eu(+- \/ x  + y   - R) + Ev z]  + [Fu(+-\/ x  + y   - R) + Fv z ]
00581  * --------------------------------    -------------------------------  = 1
00582  *               2                                      2
00583  *             Rc                                     Rd
00584  *
00585  *  The normal is the gradient of f(x,y,z) = 0 or
00586  *
00587  *      (df/dx, df/dy, df/dz)
00588  */
00589 void
00590 rt_eto_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
00591 {
00592         register struct eto_specific *eto =
00593                 (struct eto_specific *)stp->st_specific;
00594         FAST fastf_t sqrt_x2y2, efact, ffact, xcomp, ycomp, zcomp;
00595         LOCAL vect_t normp;
00596 
00597         VJOIN1( hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir );
00598 
00599         sqrt_x2y2 = sqrt( hitp->hit_vpriv[X] * hitp->hit_vpriv[X]
00600                         + hitp->hit_vpriv[Y] * hitp->hit_vpriv[Y] );
00601 
00602         efact = 2 * eto->eto_rd * eto->eto_rd * (eto->eu *
00603                 (sqrt_x2y2 - eto->eto_r) + eto->ev * hitp->hit_vpriv[Z]);
00604 
00605         ffact = 2 * eto->eto_rc * eto->eto_rc * (eto->fu *
00606                 (sqrt_x2y2 - eto->eto_r) + eto->fv * hitp->hit_vpriv[Z]);
00607 
00608         xcomp = (efact * eto->eu + ffact * eto->fu) / sqrt_x2y2;
00609 
00610         ycomp = hitp->hit_vpriv[Y] * xcomp;
00611         xcomp = hitp->hit_vpriv[X] * xcomp;
00612         zcomp = efact * eto->ev + ffact * eto->fv;
00613 
00614         VSET( normp, xcomp, ycomp, zcomp );
00615         VUNITIZE( normp );
00616         MAT3X3VEC( hitp->hit_normal, eto->eto_invR, normp );
00617 }
00618 
00619 /**
00620  *                      R T _ E T O _ C U R V E
00621  *
00622  *  Return the curvature of the eto.
00623  */
00624 void
00625 rt_eto_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
00626 {
00627         fastf_t a, b, ch, cv, dh, dv, k_circ, k_ell, phi, rad, xp,
00628                 yp1, yp2, work;
00629         register struct eto_specific *eto =
00630                 (struct eto_specific *)stp->st_specific;
00631         vect_t  Cp, Dp, Hit_Ell, Nu, Radius, Ru;
00632 
00633         a = eto->eto_rc;
00634         b = eto->eto_rd;
00635         VMOVE( Nu, eto->eto_N );
00636         VUNITIZE( Nu );
00637 
00638         /* take elliptical slice of eto at hit point */
00639         VSET( Ru, hitp->hit_vpriv[X], hitp->hit_vpriv[Y], 0. );
00640         VUNITIZE( Ru );
00641         VSCALE( Radius, Ru, eto->eto_r );
00642 
00643         /* get horizontal and vertical components of C and Rd */
00644         cv = VDOT( eto->eto_C, Nu );
00645         ch = sqrt( VDOT( eto->eto_C, eto->eto_C ) - cv * cv );
00646         /* angle between C and Nu */
00647         phi = acos( cv / MAGNITUDE(eto->eto_C) );
00648         dv = eto->eto_rd * sin(phi);
00649         dh = -eto->eto_rd * cos(phi);
00650 
00651         /* build coord system for ellipse: x,y directions are Dp,Cp */
00652         VCOMB2( Cp, ch, Ru, cv, Nu );
00653         VCOMB2( Dp, dh, Ru, dv, Nu );
00654         VUNITIZE( Cp );
00655         VUNITIZE( Dp );
00656 
00657         /* put hit point in cross sectional coordinates */
00658         VSUB2( Hit_Ell, hitp->hit_vpriv, Radius );
00659         xp = VDOT( Hit_Ell, Dp );
00660         /* yp = VDOT( Hit_Ell, Cp ); */
00661 
00662         /* calculate curvature along ellipse */
00663         /* k = y'' / (1 + y'^2) ^ 3/2 */
00664         rad = 1. / sqrt(1. - xp*xp/(a*a));
00665         yp1 = -b/(a*a)*xp*rad;
00666         yp2 = -b/(a*a)*rad*(xp*xp*rad*rad + 1.);
00667         work = 1 + yp1*yp1;
00668         k_ell = yp2 / (work*sqrt(work));
00669 
00670         /* calculate curvature along radial circle */
00671         k_circ = -1. / MAGNITUDE(Radius);
00672 
00673         if (fabs(k_ell) < fabs(k_circ)) {
00674                 /* use 1st deriv for principle dir of curvature */
00675                 VCOMB2( cvp->crv_pdir, xp, Dp, yp1, Cp );
00676                 cvp->crv_c1 = k_ell;
00677                 cvp->crv_c2 = k_circ;
00678         } else {
00679                 VCROSS( cvp->crv_pdir, hitp->hit_normal, Radius );
00680                 cvp->crv_c1 = k_circ;
00681                 cvp->crv_c2 = k_ell;
00682         }
00683         VUNITIZE( cvp->crv_pdir );
00684 }
00685 
00686 /**
00687  *                      R T _ E T O _ U V
00688  */
00689 void
00690 rt_eto_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
00691 {
00692         fastf_t horz, theta_u, theta_v, vert;
00693         vect_t  Hit_Ell, Nu, Radius, Ru;
00694 
00695         register struct eto_specific *eto =
00696                 (struct eto_specific *)stp->st_specific;
00697 
00698         /* take elliptical slice of eto at hit point */
00699         VSET( Ru, hitp->hit_vpriv[X], hitp->hit_vpriv[Y], 0. );
00700         VUNITIZE( Ru );
00701         VSCALE( Radius, Ru, eto->eto_r );
00702         /* put cross sectional ellipse at origin */
00703         VSUB2( Hit_Ell, hitp->hit_vpriv, Radius );
00704 
00705         /* find angle between Ru and Hit_Ell
00706            (better for circle than ellipse...) */
00707         VMOVE( Nu, eto->eto_N );
00708         VUNITIZE( Nu );
00709         vert = VDOT( Hit_Ell, Nu );
00710         horz = VDOT( Hit_Ell, Ru );
00711         theta_u = atan2(vert, -horz);   /* tuck seam on eto inner diameter */
00712 
00713         /* find angle between hitp and x axis */
00714         theta_v = atan2(hitp->hit_vpriv[Y], hitp->hit_vpriv[X]);
00715 
00716         /* normalize to [0, 2pi] */
00717         if (theta_u < 0.)
00718                 theta_u += bn_twopi;
00719         if (theta_v < 0.)
00720                 theta_v += bn_twopi;
00721 
00722         /* normalize to [0, 1] */
00723         uvp->uv_u = theta_u/bn_twopi;
00724         uvp->uv_v = theta_v/bn_twopi;
00725         uvp->uv_du = uvp->uv_dv = 0;
00726 }
00727 
00728 /**
00729  *                      R T _ E T O _ F R E E
00730  */
00731 void
00732 rt_eto_free(struct soltab *stp)
00733 {
00734         register struct eto_specific *eto =
00735                 (struct eto_specific *)stp->st_specific;
00736 
00737         bu_free( (char *)eto, "eto_specific");
00738 }
00739 
00740 int
00741 rt_eto_class(void)
00742 {
00743         return(0);
00744 }
00745 
00746 /**
00747  *                      R T _ E T O _ P L O T
00748  *
00749  * The ETO has the following input fields:
00750  *      eto_V   V from origin to center
00751  *      eto_r   Radius scalar
00752  *      eto_N   Normal to plane of eto
00753  *      eto_C   Semimajor axis (vector) of eto cross section
00754  *      eto_rd  Semiminor axis length (scalar) of eto cross section
00755  *
00756  */
00757 int
00758 rt_eto_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
00759 {
00760         fastf_t         a, b;   /* axis lengths of ellipse */
00761         fastf_t         ang, ch, cv, dh, dv, ntol, dtol, phi, theta;
00762         fastf_t         *eto_ells;
00763         int             i, j, npts, nells;
00764         point_t         *ell;   /* array of ellipse points */
00765         point_t         Ell_V;  /* vertex of an ellipse */
00766         point_t         *rt_mk_ell();
00767         struct rt_eto_internal  *tip;
00768         vect_t          Au, Bu, Nu, Cp, Dp, Xu;
00769 
00770         RT_CK_DB_INTERNAL(ip);
00771         tip = (struct rt_eto_internal *)ip->idb_ptr;
00772         RT_ETO_CK_MAGIC(tip);
00773 
00774         a = MAGNITUDE( tip->eto_C );
00775         b = tip->eto_rd;
00776 
00777         if ( NEAR_ZERO(tip->eto_r, 0.0001) || NEAR_ZERO(b, 0.0001)
00778                 || NEAR_ZERO(a, 0.0001)) {
00779                 bu_log("eto_plot: r, rd, or rc zero length\n");
00780                 return(1);
00781         }
00782 
00783         /* Establish tolerances */
00784         if( ttol->rel <= 0.0 || ttol->rel >= 1.0 )  {
00785                 dtol = 0.0;             /* none */
00786         } else {
00787                 /*
00788                  * Convert relative to absolute by scaling smallest of
00789                  * radius and semi-minor axis
00790                  */
00791                 if (tip->eto_r < b)
00792                         dtol = ttol->rel * 2 * tip->eto_r;
00793                 else
00794                         dtol = ttol->rel * 2 * b;
00795         }
00796         if( ttol->abs <= 0.0 )  {
00797                 if( dtol <= 0.0 )  {
00798                         /* No tolerance given, use a default */
00799                         if (tip->eto_r < b)
00800                                 dtol = 2 * 0.10 * tip->eto_r;   /* 10% */
00801                         else
00802                                 dtol = 2 * 0.10 * b;    /* 10% */
00803                 } else {
00804                         /* Use absolute-ized relative tolerance */
00805                 }
00806         } else {
00807                 /* Absolute tolerance was given, pick smaller */
00808                 if( ttol->rel <= 0.0 || dtol > ttol->abs )
00809                         dtol = ttol->abs;
00810         }
00811         /* To ensure normal tolerance, remain below this angle */
00812         if( ttol->norm > 0.0 )
00813                 ntol = ttol->norm;
00814         else
00815                 /* tolerate everything */
00816                 ntol = bn_pi;
00817 
00818         /* (x, y) coords for an ellipse */
00819         ell = rt_mk_ell( &npts, a, b, dtol, ntol );
00820         /* generate coordinate axes */
00821         VMOVE( Nu, tip->eto_N );
00822         VUNITIZE( Nu );                 /* z axis of coord sys */
00823         bn_vec_ortho( Bu, Nu );         /* x axis */
00824         VUNITIZE( Bu );
00825         VCROSS( Au, Nu, Bu );           /* y axis */
00826 
00827         /* number of segments required in eto circles */
00828         nells = rt_num_circular_segments( dtol, tip->eto_r );
00829         theta = bn_twopi / nells;       /* put ellipse every theta rads */
00830         /* get horizontal and vertical components of C and Rd */
00831         cv = VDOT( tip->eto_C, Nu );
00832         ch = sqrt( VDOT( tip->eto_C, tip->eto_C ) - cv * cv );
00833         /* angle between C and Nu */
00834         phi = acos( cv / MAGNITUDE(tip->eto_C) );
00835         dv = tip->eto_rd * sin(phi);
00836         dh = -tip->eto_rd * cos(phi);
00837 
00838         /* make sure ellipse doesn't overlap itself when revolved */
00839         if (ch > tip->eto_r || dh > tip->eto_r) {
00840                 bu_log("eto_plot: revolved ellipse overlaps itself\n");
00841                 return(1);
00842         }
00843 
00844         /* get memory for nells ellipses */
00845         eto_ells = (fastf_t *)bu_malloc(nells * npts * sizeof(point_t), "ells[]");
00846 
00847         /* place each ellipse properly to make eto */
00848         for (i = 0, ang = 0.; i < nells; i++, ang += theta) {
00849                 /* direction of current ellipse */
00850                 VCOMB2( Xu, cos(ang), Bu, sin(ang), Au );
00851                 VUNITIZE( Xu );
00852                 /* vertex of ellipse */
00853                 VJOIN1( Ell_V, tip->eto_V, tip->eto_r, Xu );
00854                 /* coord system for ellipse: x,y directions are Dp,Cp */
00855                 VCOMB2( Cp, ch, Xu, cv, Nu );
00856                 VCOMB2( Dp, dh, Xu, dv, Nu );
00857                 VUNITIZE( Cp );
00858                 VUNITIZE( Dp );
00859 
00860 /* convert 2D address to index into 1D array */
00861 #define ETO_PT(www,lll) ((((www)%nells)*npts)+((lll)%npts))
00862 #define ETO_PTA(ww,ll)  (&eto_ells[ETO_PT(ww,ll)*3])
00863 #define ETO_NMA(ww,ll)  (norms[ETO_PT(ww,ll)])
00864 
00865                 /* make ellipse */
00866                 for (j = 0; j < npts; j++) {
00867                         VJOIN2( ETO_PTA(i,j),
00868                                 Ell_V, ell[j][X], Dp, ell[j][Y], Cp );
00869                 }
00870         }
00871 
00872         /* draw ellipses */
00873         for (i = 0; i < nells; i++) {
00874                 RT_ADD_VLIST( vhead, ETO_PTA(i,npts-1), BN_VLIST_LINE_MOVE );
00875                 for( j = 0; j < npts; j++ )
00876                         RT_ADD_VLIST( vhead, ETO_PTA(i,j), BN_VLIST_LINE_DRAW );
00877         }
00878 
00879         /* draw connecting circles */
00880         for (i = 0; i < npts; i++) {
00881                 RT_ADD_VLIST( vhead, ETO_PTA(nells-1,i), BN_VLIST_LINE_MOVE );
00882                 for( j = 0; j < nells; j++ )
00883                         RT_ADD_VLIST( vhead, ETO_PTA(j,i), BN_VLIST_LINE_DRAW );
00884         }
00885 
00886         bu_free( (char *)eto_ells, "ells[]" );
00887         return(0);
00888 }
00889 
00890 /**
00891  *      R T _ E L L 4
00892  *
00893  *      Approximate one fourth (1st quadrant) of an ellipse with line
00894  *      segments.  The initial single segment is broken at the point
00895  *      farthest from the ellipse if that point is not aleady within the
00896  *      distance and normal error tolerances.  The two resulting segments are
00897  *      passed recursively to this routine until each segment is within
00898  *      tolerance.
00899  */
00900 int
00901 rt_ell4(struct rt_pt_node *pts, fastf_t a, fastf_t b, fastf_t dtol, fastf_t ntol)
00902 {
00903         fastf_t dist, intr, m, theta0, theta1;
00904         int     n;
00905         point_t mpt, p0, p1;
00906         vect_t  norm_line, norm_ell;
00907         struct rt_pt_node *new, *rt_ptalloc(void);
00908 
00909         /* endpoints of segment approximating ellipse */
00910         VMOVE( p0, pts->p );
00911         VMOVE( p1, pts->next->p );
00912         /* slope and intercept of segment */
00913         m = ( p1[X] - p0[X] ) / ( p1[Y] - p0[Y] );
00914         intr = p0[X] - m * p0[Y];
00915         /* point on ellipse with max dist between ellipse and line */
00916         mpt[Y] = a / sqrt( b*b / (m*m*a*a) + 1 );
00917         mpt[X] = b * sqrt( 1 - mpt[Y] * mpt[Y] / (a*a) );
00918         mpt[Z] = 0;
00919         /* max distance between that point and line */
00920         dist = fabs( m * mpt[Y] - mpt[X] + intr ) / sqrt( m * m + 1 );
00921         /* angles between normal of line and of ellipse at line endpoints */
00922         VSET( norm_line, m, -1., 0.);
00923         VSET( norm_ell, b * b * p0[Y], a * a * p0[X], 0. );
00924         VUNITIZE( norm_line );
00925         VUNITIZE( norm_ell );
00926         theta0 = fabs( acos( VDOT( norm_line, norm_ell )));
00927         VSET( norm_ell, b * b * p1[Y], a * a * p1[X], 0. );
00928         VUNITIZE( norm_ell );
00929         theta1 = fabs( acos( VDOT( norm_line, norm_ell )));
00930         /* split segment at widest point if not within error tolerances */
00931         if ( dist > dtol || theta0 > ntol || theta1 > ntol ) {
00932                 /* split segment */
00933                 new = rt_ptalloc();
00934                 VMOVE( new->p, mpt );
00935                 new->next = pts->next;
00936                 pts->next = new;
00937                 /* keep track of number of pts added */
00938                 n = 1;
00939                 /* recurse on first new segment */
00940                 n += rt_ell4( pts, a, b, dtol, ntol );
00941                 /* recurse on second new segment */
00942                 n += rt_ell4( new, a, b, dtol, ntol );
00943         } else
00944                 n  = 0;
00945         return( n );
00946 }
00947 
00948 /**
00949  *      R T _ M K _ E L L
00950  *
00951  *      Return pointer an array of points approximating an ellipse
00952  *      with semi-major and semi-minor axes a and b.  The line
00953  *      segments fall within the normal and distance tolerances
00954  *      of ntol and dtol.
00955  */
00956 point_t *
00957 rt_mk_ell( n, a, b, dtol, ntol )
00958 int     *n;
00959 fastf_t a, b, dtol, ntol;
00960 {
00961         int             i;
00962         point_t         *ell;
00963         struct rt_pt_node       *ell_quad, *oldpos, *pos, *rt_ptalloc(void);
00964 
00965         ell_quad = rt_ptalloc();
00966         VSET( ell_quad->p, b, 0., 0. );
00967         ell_quad->next = rt_ptalloc();
00968         VSET( ell_quad->next->p, 0., a, 0. );
00969         ell_quad->next->next = NULL;
00970 
00971         *n = rt_ell4( ell_quad, a, b, dtol, ntol );
00972         ell = (point_t *)bu_malloc(4*(*n+1)*sizeof(point_t), "rt_mk_ell pts");
00973 
00974         /* put 1st quad of ellipse into an array */
00975         pos = ell_quad;
00976         for (i = 0; i < *n+2; i++) {
00977                 VMOVE( ell[i], pos->p );
00978                 oldpos = pos;
00979                 pos = pos->next;
00980                 bu_free( (char *)oldpos, "rt_pt_node" );
00981         }
00982         /* mirror 1st quad to make 2nd */
00983         for (i = (*n+1)+1; i < 2*(*n+1); i++) {
00984                 VMOVE( ell[i], ell[(*n*2+2)-i] );
00985                 ell[i][X] = -ell[i][X];
00986         }
00987         /* mirror 2nd quad to make 3rd */
00988         for (i = 2*(*n+1); i < 3*(*n+1); i++) {
00989                 VMOVE( ell[i], ell[i-(*n*2+2)] );
00990                 ell[i][X] = -ell[i][X];
00991                 ell[i][Y] = -ell[i][Y];
00992         }
00993         /* mirror 3rd quad to make 4th */
00994         for (i = 3*(*n+1); i < 4*(*n+1); i++) {
00995                 VMOVE( ell[i], ell[i-(*n*2+2)] );
00996                 ell[i][X] = -ell[i][X];
00997                 ell[i][Y] = -ell[i][Y];
00998         }
00999         *n = 4*(*n + 1);
01000         return(ell);
01001 }
01002 
01003 /**
01004  *                      R T _ E T O _ T E S S
01005  */
01006 int
01007 rt_eto_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct rt_tess_tol *ttol, const struct bn_tol *tol)
01008 {
01009         fastf_t         a, b;   /* axis lengths of ellipse */
01010         fastf_t         ang, ch, cv, dh, dv, ntol, dtol, phi, theta;
01011         fastf_t         *eto_ells = NULL;
01012         int             i, j, nfaces, npts, nells;
01013         point_t         *ell = NULL;    /* array of ellipse points */
01014         point_t         Ell_V;  /* vertex of an ellipse */
01015         point_t         *rt_mk_ell();
01016         struct rt_eto_internal  *tip;
01017         struct shell    *s;
01018         struct vertex   **verts = NULL;
01019         struct faceuse  **faces = NULL;
01020         struct vertex   **vertp[4];
01021         vect_t          Au, Bu, Nu, Cp, Dp, Xu;
01022         vect_t          *norms = NULL;  /* normal vectors for each vertex */
01023         int             fail=0;
01024 
01025         RT_CK_DB_INTERNAL(ip);
01026         tip = (struct rt_eto_internal *)ip->idb_ptr;
01027         RT_ETO_CK_MAGIC(tip);
01028 
01029         a = MAGNITUDE( tip->eto_C );
01030         b = tip->eto_rd;
01031 
01032         if ( NEAR_ZERO(tip->eto_r, 0.0001) || NEAR_ZERO(b, 0.0001)
01033                 || NEAR_ZERO(a, 0.0001)) {
01034                 bu_log("eto_tess: r, rd, or rc zero length\n");
01035                 fail = (-2);
01036                 goto failure;
01037         }
01038 
01039         /* Establish tolerances */
01040         if( ttol->rel <= 0.0 || ttol->rel >= 1.0 )  {
01041                 dtol = 0.0;             /* none */
01042         } else {
01043                 /*
01044                  * Convert relative to absolute by scaling smallest of
01045                  * radius and semi-minor axis
01046                  */
01047                 if (tip->eto_r < b)
01048                         dtol = ttol->rel * 2 * tip->eto_r;
01049                 else
01050                         dtol = ttol->rel * 2 * b;
01051         }
01052         if( ttol->abs <= 0.0 )  {
01053                 if( dtol <= 0.0 )  {
01054                         /* No tolerance given, use a default */
01055                         if (tip->eto_r < b)
01056                                 dtol = 2 * 0.10 * tip->eto_r;   /* 10% */
01057                         else
01058                                 dtol = 2 * 0.10 * b;    /* 10% */
01059                 } else {
01060                         /* Use absolute-ized relative tolerance */
01061                 }
01062         } else {
01063                 /* Absolute tolerance was given, pick smaller */
01064                 if( ttol->rel <= 0.0 || dtol > ttol->abs )
01065                         dtol = ttol->abs;
01066         }
01067         /* To ensure normal tolerance, remain below this angle */
01068         if( ttol->norm > 0.0 )
01069                 ntol = ttol->norm;
01070         else
01071                 /* tolerate everything */
01072                 ntol = bn_pi;
01073 
01074         /* (x, y) coords for an ellipse */
01075         ell = rt_mk_ell( &npts, a, b, dtol, ntol );
01076         /* generate coordinate axes */
01077         VMOVE( Nu, tip->eto_N );
01078         VUNITIZE( Nu );                 /* z axis of coord sys */
01079         bn_vec_ortho( Bu, Nu );         /* x axis */
01080         VUNITIZE( Bu );
01081         VCROSS( Au, Nu, Bu );           /* y axis */
01082 
01083         /* number of segments required in eto circles */
01084         nells = rt_num_circular_segments( dtol, tip->eto_r );
01085         theta = bn_twopi / nells;       /* put ellipse every theta rads */
01086         /* get horizontal and vertical components of C and Rd */
01087         cv = VDOT( tip->eto_C, Nu );
01088         ch = sqrt( VDOT( tip->eto_C, tip->eto_C ) - cv * cv );
01089         /* angle between C and Nu */
01090         phi = acos( cv / MAGNITUDE(tip->eto_C) );
01091         dv = tip->eto_rd * sin(phi);
01092         dh = -tip->eto_rd * cos(phi);
01093 
01094         /* make sure ellipse doesn't overlap itself when revolved */
01095         if (ch > tip->eto_r || dh > tip->eto_r) {
01096                 bu_log("eto_tess: revolved ellipse overlaps itself\n");
01097                 fail = (-3);
01098                 goto failure;
01099         }
01100 
01101         /* get memory for nells ellipses */
01102         eto_ells = (fastf_t *)bu_malloc(nells * npts * sizeof(point_t), "ells[]");
01103         norms = (vect_t *)bu_calloc( nells*npts , sizeof( vect_t ) , "rt_eto_tess: norms" );
01104 
01105         /* place each ellipse properly to make eto */
01106         for (i = 0, ang = 0.; i < nells; i++, ang += theta) {
01107                 /* direction of current ellipse */
01108                 VCOMB2( Xu, cos(ang), Bu, sin(ang), Au );
01109                 VUNITIZE( Xu );
01110                 /* vertex of ellipse */
01111                 VJOIN1( Ell_V, tip->eto_V, tip->eto_r, Xu );
01112                 /* coord system for ellipse: x,y directions are Dp,Cp */
01113                 VCOMB2( Cp, ch, Xu, cv, Nu );
01114                 VCOMB2( Dp, dh, Xu, dv, Nu );
01115                 VUNITIZE( Cp );
01116                 VUNITIZE( Dp );
01117                 /* make ellipse */
01118                 for (j = 0; j < npts; j++) {
01119                         VJOIN2( ETO_PTA(i,j),
01120                                 Ell_V, ell[j][X], Dp, ell[j][Y], Cp );
01121                         VBLEND2( ETO_NMA(i,j),
01122                                 a*a*ell[j][X], Dp , b*b*ell[j][Y], Cp );
01123                         VUNITIZE( ETO_NMA(i,j) );
01124                 }
01125         }
01126 
01127         *r = nmg_mrsv( m );     /* Make region, empty shell, vertex */
01128         s = BU_LIST_FIRST(shell, &(*r)->s_hd);
01129 
01130         verts = (struct vertex **)bu_calloc( npts*nells, sizeof(struct vertex *),
01131                 "rt_eto_tess *verts[]" );
01132         faces = (struct faceuse **)bu_calloc( npts*nells, sizeof(struct faceuse *),
01133                 "rt_eto_tess *faces[]" );
01134 
01135         /* Build the topology of the eto */
01136         nfaces = 0;
01137         for( i = 0; i < nells; i++ )  {
01138                 for( j = 0; j < npts; j++ )  {
01139                         vertp[0] = &verts[ ETO_PT(i+0,j+0) ];
01140                         vertp[1] = &verts[ ETO_PT(i+0,j+1) ];
01141                         vertp[2] = &verts[ ETO_PT(i+1,j+1) ];
01142                         vertp[3] = &verts[ ETO_PT(i+1,j+0) ];
01143                         if( (faces[nfaces++] = nmg_cmface( s, vertp, 4 )) == (struct faceuse *)0 )  {
01144                                 bu_log("rt_eto_tess() nmg_cmface failed, i=%d/%d, j=%d/%d\n",
01145                                         i, nells, j, npts );
01146                                 nfaces--;
01147                         }
01148                 }
01149         }
01150 
01151         /* Associate vertex geometry */
01152         for( i = 0; i < nells; i++ )  {
01153                 for( j = 0; j < npts; j++ )  {
01154                         nmg_vertex_gv( verts[ETO_PT(i,j)], ETO_PTA(i,j) );
01155                 }
01156         }
01157 
01158         /* Associate face geometry */
01159         for( i=0; i < nfaces; i++ )  {
01160                 if( nmg_fu_planeeqn( faces[i], tol ) < 0 )
01161                 {
01162                         fail = (-1);
01163                         goto failure;
01164                 }
01165         }
01166 
01167         /* associate vertexuse normals */
01168         for( i=0 ; i<nells ; i++ )
01169         {
01170                 for( j=0 ; j<npts ; j++ )
01171                 {
01172                         struct vertexuse *vu;
01173                         vect_t rev_norm;
01174 
01175                         VREVERSE( rev_norm , ETO_NMA(i,j) );
01176 
01177                         NMG_CK_VERTEX( verts[ETO_PT(i,j)] );
01178 
01179                         for( BU_LIST_FOR( vu , vertexuse , &verts[ETO_PT(i,j)]->vu_hd ) )
01180                         {
01181                                 struct faceuse *fu;
01182 
01183                                 NMG_CK_VERTEXUSE( vu );
01184 
01185                                 fu = nmg_find_fu_of_vu( vu );
01186                                 NMG_CK_FACEUSE( fu );
01187 
01188                                 if( fu->orientation == OT_SAME )
01189                                         nmg_vertexuse_nv( vu , ETO_NMA(i,j) );
01190                                 else if( fu->orientation == OT_OPPOSITE )
01191                                         nmg_vertexuse_nv( vu , rev_norm );
01192                         }
01193                 }
01194         }
01195 
01196         /* Compute "geometry" for region and shell */
01197         nmg_region_a( *r, tol );
01198 
01199 failure:
01200         bu_free( (char *)ell, "rt_mk_ell pts" );
01201         bu_free( (char *)eto_ells, "ells[]" );
01202         bu_free( (char *)verts, "rt_eto_tess *verts[]" );
01203         bu_free( (char *)faces, "rt_eto_tess *faces[]" );
01204         bu_free( (char *)norms, "rt_eto_tess: norms[]" );
01205 
01206         return( fail );
01207 }
01208 
01209 /**
01210  *                      R T _ E T O _ I M P O R T
01211  *
01212  *  Import a eto from the database format to the internal format.
01213  *  Apply modeling transformations at the same time.
01214  */
01215 int
01216 rt_eto_import(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
01217 {
01218         struct rt_eto_internal  *tip;
01219         union record            *rp;
01220 
01221         BU_CK_EXTERNAL( ep );
01222         rp = (union record *)ep->ext_buf;
01223         /* Check record type */
01224         if( rp->u_id != ID_SOLID )  {
01225                 bu_log("rt_eto_import: defective record\n");
01226                 return(-1);
01227         }
01228 
01229         RT_CK_DB_INTERNAL( ip );
01230         ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
01231         ip->idb_type = ID_ETO;
01232         ip->idb_meth = &rt_functab[ID_ETO];
01233         ip->idb_ptr = bu_malloc(sizeof(struct rt_eto_internal), "rt_eto_internal");
01234         tip = (struct rt_eto_internal *)ip->idb_ptr;
01235         tip->eto_magic = RT_ETO_INTERNAL_MAGIC;
01236 
01237         /* Apply modeling transformations */
01238         MAT4X3PNT( tip->eto_V, mat, &rp->s.s_values[0*3] );
01239         MAT4X3VEC( tip->eto_N, mat, &rp->s.s_values[1*3] );
01240         MAT4X3VEC( tip->eto_C, mat, &rp->s.s_values[2*3] );
01241         tip->eto_r  = rp->s.s_values[3*3] / mat[15];
01242         tip->eto_rd = rp->s.s_values[3*3+1] / mat[15];
01243 
01244         if( tip->eto_r < SMALL || tip->eto_rd < SMALL )  {
01245                 bu_log("rt_eto_import:  zero length R or Rd vector\n");
01246                 return(-1);
01247         }
01248 
01249         return(0);              /* OK */
01250 }
01251 
01252 /**
01253  *                      R T _ E T O _ E X P O R T
01254  *
01255  *  The name will be added by the caller.
01256  */
01257 int
01258 rt_eto_export(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
01259 {
01260         struct rt_eto_internal  *tip;
01261         union record            *eto;
01262 
01263         RT_CK_DB_INTERNAL(ip);
01264         if( ip->idb_type != ID_ETO )  return(-1);
01265         tip = (struct rt_eto_internal *)ip->idb_ptr;
01266         RT_ETO_CK_MAGIC(tip);
01267 
01268         BU_CK_EXTERNAL(ep);
01269         ep->ext_nbytes = sizeof(union record);
01270         ep->ext_buf = (genptr_t)bu_calloc( 1, ep->ext_nbytes, "eto external");
01271         eto = (union record *)ep->ext_buf;
01272 
01273         eto->s.s_id = ID_SOLID;
01274         eto->s.s_type = ETO;
01275 
01276         if (MAGNITUDE(tip->eto_C) < RT_LEN_TOL
01277                 || MAGNITUDE(tip->eto_N) < RT_LEN_TOL
01278                 || tip->eto_r < RT_LEN_TOL
01279                 || tip->eto_rd < RT_LEN_TOL) {
01280                 bu_log("rt_eto_export: not all dimensions positive!\n");
01281                 return(-1);
01282         }
01283 
01284         if (tip->eto_rd > MAGNITUDE(tip->eto_C) ) {
01285                 bu_log("rt_eto_export: semi-minor axis cannot be longer than semi-major axis!\n");
01286                 return(-1);
01287         }
01288 
01289         /* Warning:  type conversion */
01290         VSCALE( &eto->s.s_values[0*3], tip->eto_V, local2mm );
01291         VSCALE( &eto->s.s_values[1*3], tip->eto_N, local2mm );
01292         VSCALE( &eto->s.s_values[2*3], tip->eto_C, local2mm );
01293         eto->s.s_values[3*3] = tip->eto_r * local2mm;
01294         eto->s.s_values[3*3+1] = tip->eto_rd * local2mm;
01295 
01296         return(0);
01297 }
01298 
01299 /**
01300  *                      R T _ E T O _ I M P O R T 5
01301  *
01302  *  Import a eto from the database format to the internal format.
01303  *  Apply modeling transformations at the same time.
01304  */
01305 int
01306 rt_eto_import5(struct rt_db_internal *ip, const struct bu_external *ep, register const fastf_t *mat, const struct db_i *dbip)
01307 {
01308         struct rt_eto_internal  *tip;
01309         fastf_t                 vec[11];
01310 
01311         BU_CK_EXTERNAL( ep );
01312 
01313         BU_ASSERT_LONG( ep->ext_nbytes, ==, SIZEOF_NETWORK_DOUBLE * 11 );
01314 
01315         RT_CK_DB_INTERNAL( ip );
01316         ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
01317         ip->idb_type = ID_ETO;
01318         ip->idb_meth = &rt_functab[ID_ETO];
01319         ip->idb_ptr = bu_malloc(sizeof(struct rt_eto_internal), "rt_eto_internal");
01320 
01321         tip = (struct rt_eto_internal *)ip->idb_ptr;
01322         tip->eto_magic = RT_ETO_INTERNAL_MAGIC;
01323 
01324         /* Convert from database (network) to internal (host) format */
01325         ntohd( (unsigned char *)vec, ep->ext_buf, 11 );
01326 
01327         /* Apply modeling transformations */
01328         MAT4X3PNT( tip->eto_V, mat, &vec[0*3] );
01329         MAT4X3VEC( tip->eto_N, mat, &vec[1*3] );
01330         MAT4X3VEC( tip->eto_C, mat, &vec[2*3] );
01331         tip->eto_r  = vec[3*3] / mat[15];
01332         tip->eto_rd = vec[3*3+1] / mat[15];
01333 
01334         if( tip->eto_r < SMALL || tip->eto_rd < SMALL )  {
01335                 bu_log("rt_eto_import:  zero length R or Rd vector\n");
01336                 return(-1);
01337         }
01338 
01339         return(0);              /* OK */
01340 }
01341 
01342 /**
01343  *                      R T _ E T O _ E X P O R T 5
01344  *
01345  *  The name will be added by the caller.
01346  */
01347 int
01348 rt_eto_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
01349 {
01350         struct rt_eto_internal  *tip;
01351         fastf_t                 vec[11];
01352 
01353         RT_CK_DB_INTERNAL(ip);
01354         if( ip->idb_type != ID_ETO )  return(-1);
01355         tip = (struct rt_eto_internal *)ip->idb_ptr;
01356         RT_ETO_CK_MAGIC(tip);
01357 
01358         BU_CK_EXTERNAL(ep);
01359         ep->ext_nbytes = SIZEOF_NETWORK_DOUBLE * 11;
01360         ep->ext_buf = (genptr_t)bu_malloc( ep->ext_nbytes, "eto external");
01361 
01362         if (MAGNITUDE(tip->eto_C) < RT_LEN_TOL
01363                 || MAGNITUDE(tip->eto_N) < RT_LEN_TOL
01364                 || tip->eto_r < RT_LEN_TOL
01365                 || tip->eto_rd < RT_LEN_TOL) {
01366                 bu_log("rt_eto_export: not all dimensions positive!\n");
01367                 return(-1);
01368         }
01369 
01370         if (tip->eto_rd > MAGNITUDE(tip->eto_C) ) {
01371                 bu_log("rt_eto_export: semi-minor axis cannot be longer than semi-major axis!\n");
01372                 return(-1);
01373         }
01374 
01375         /* scale 'em into local buffer */
01376         VSCALE( &vec[0*3], tip->eto_V, local2mm );
01377         VSCALE( &vec[1*3], tip->eto_N, local2mm );
01378         VSCALE( &vec[2*3], tip->eto_C, local2mm );
01379         vec[3*3] = tip->eto_r * local2mm;
01380         vec[3*3+1] = tip->eto_rd * local2mm;
01381 
01382         /* Convert from internal (host) to database (network) format */
01383         htond( ep->ext_buf, (unsigned char *)vec, 11 );
01384 
01385         return(0);
01386 }
01387 
01388 /**
01389  *                      R T _ E T O _ D E S C R I B E
01390  *
01391  *  Make human-readable formatted presentation of this solid.
01392  *  First line describes type of solid.
01393  *  Additional lines are indented one tab, and give parameter values.
01394  */
01395 int
01396 rt_eto_describe(struct bu_vls *str, const struct rt_db_internal *ip, int verbose, double mm2local)
01397 {
01398         register struct rt_eto_internal *tip =
01399                 (struct rt_eto_internal *)ip->idb_ptr;
01400         char                            buf[256];
01401 
01402         RT_ETO_CK_MAGIC(tip);
01403         bu_vls_strcat( str, "Elliptical Torus (ETO)\n");
01404 
01405         sprintf(buf, "\tV (%g, %g, %g)\n",
01406                 INTCLAMP(tip->eto_V[X] * mm2local),
01407                 INTCLAMP(tip->eto_V[Y] * mm2local),
01408                 INTCLAMP(tip->eto_V[Z] * mm2local) );
01409         bu_vls_strcat( str, buf );
01410 
01411         sprintf(buf, "\tN=(%g, %g, %g)\n",
01412                 INTCLAMP(tip->eto_N[X] * mm2local),
01413                 INTCLAMP(tip->eto_N[Y] * mm2local),
01414                 INTCLAMP(tip->eto_N[Z] * mm2local) );
01415         bu_vls_strcat( str, buf );
01416 
01417         sprintf(buf, "\tC=(%g, %g, %g) mag=%g\n",
01418                 INTCLAMP(tip->eto_C[X] * mm2local),
01419                 INTCLAMP(tip->eto_C[Y] * mm2local),
01420                 INTCLAMP(tip->eto_C[Z] * mm2local),
01421                 INTCLAMP(MAGNITUDE(tip->eto_C) * mm2local) );
01422         bu_vls_strcat( str, buf );
01423 
01424         sprintf(buf, "\tr=%g\n", INTCLAMP(tip->eto_r * mm2local));
01425         bu_vls_strcat( str, buf );
01426 
01427         sprintf(buf, "\td=%g\n", INTCLAMP(tip->eto_rd * mm2local));
01428         bu_vls_strcat( str, buf );
01429 
01430         return(0);
01431 }
01432 
01433 /**
01434  *                      R T _ E T O _ I F R E E
01435  *
01436  *  Free the storage associated with the rt_db_internal version of this solid.
01437  */
01438 void
01439 rt_eto_ifree(struct rt_db_internal *ip)
01440 {
01441         register struct rt_eto_internal *tip;
01442 
01443         RT_CK_DB_INTERNAL(ip);
01444         tip = (struct rt_eto_internal *)ip->idb_ptr;
01445         RT_ETO_CK_MAGIC(tip);
01446 
01447         bu_free( (char *)tip, "eto ifree" );
01448         ip->idb_ptr = GENPTR_NULL;      /* sanity */
01449 }
01450 
01451 /*@}*/
01452 /*
01453  * Local Variables:
01454  * mode: C
01455  * tab-width: 8
01456  * c-basic-offset: 4
01457  * indent-tabs-mode: t
01458  * End:
01459  * ex: shiftwidth=4 tabstop=8
01460  */

Generated on Mon Sep 18 01:24:50 2006 for BRL-CAD by  doxygen 1.4.6