00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
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
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145 struct eto_specific {
00146 vect_t eto_V;
00147 vect_t eto_N;
00148 vect_t eto_C;
00149 fastf_t eto_r;
00150 fastf_t eto_rc;
00151 fastf_t eto_rd;
00152 mat_t eto_R;
00153 mat_t eto_invR;
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
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
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;
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
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 );
00211 bn_vec_ortho( Bu, Nu );
00212 VUNITIZE( Bu );
00213 VCROSS( Au, Nu, Bu );
00214 VMOVE( Cu, tip->eto_C );
00215 VUNITIZE( Cu );
00216
00217
00218 cv = VDOT( eto->eto_C, Nu );
00219 ch = sqrt( VDOT( eto->eto_C, eto->eto_C ) - cv * cv );
00220
00221 phi = acos( cv / eto->eto_rc );
00222 dh = eto->eto_rd * cos(phi);
00223
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 ) );
00231 eto->eu = sqrt( 1.0 - eto->ev * eto->ev );
00232 eto->fu = -eto->ev;
00233 eto->fv = eto->eu;
00234
00235
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
00246
00247
00248
00249
00250
00251
00252
00253
00254 VSET( P, 1.0, 0, 0 );
00255 VCROSS( w1, Nu, 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
00263 VSET( P, 0, 1.0, 0 );
00264 VCROSS( w1, Nu, 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
00272 VSET( P, 0, 0, 1.0 );
00273 VCROSS( w1, Nu, 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);
00281 }
00282
00283
00284
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
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
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;
00342 LOCAL vect_t pprime;
00343 LOCAL vect_t work;
00344 LOCAL bn_poly_t C;
00345 LOCAL bn_complex_t val[4];
00346 LOCAL double k[4];
00347 register int i;
00348 LOCAL int j;
00349 LOCAL vect_t cor_pprime;
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
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
00361
00362
00363
00364
00365
00366
00367 cor_proj = VDOT( pprime, dprime );
00368 VSCALE( cor_pprime, dprime, cor_proj );
00369 VSUB2( cor_pprime, pprime, cor_pprime );
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
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
00425
00426
00427
00428
00429
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
00440
00441 C.dgr=4;
00442 C.cf[4] = (A4*A4 - A7*A7*C1);
00443 C.cf[3] = (2*A4*A5 - A7*A7*C2 - 2*A7*A8*C1);
00444 C.cf[2] = (2*A4*A6 + A5*A5 - A7*A7*C3 - 2*A7*A8*C2 - A8*A8*C1);
00445 C.cf[1] = (2*A5*A6 - 2*A7*A8*C3 - A8*A8*C2);
00446 C.cf[0] = term;
00447
00448
00449
00450
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);
00467 }
00468
00469
00470
00471
00472
00473
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
00481 for( j = 0; j < i; ++j )
00482 k[j] -= cor_proj;
00483
00484
00485 switch( i ) {
00486 case 0:
00487 return(0);
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);
00493
00494 case 2:
00495 {
00496
00497 FAST fastf_t u;
00498 if( (u=k[0]) < k[1] ) {
00499
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
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
00516 k[n] = k[n+1];
00517 k[n+1] = u;
00518 }
00519 }
00520 }
00521 }
00522 break;
00523 }
00524
00525
00526
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
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);
00538
00539
00540
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);
00549 }
00550
00551 #define SEG_MISS(SEG) (SEG).seg_stp=(struct soltab *) 0;
00552
00553
00554
00555
00556
00557 void
00558 rt_eto_vshot(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
00559
00560
00561
00562
00563
00564 {
00565 }
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
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
00621
00622
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
00639 VSET( Ru, hitp->hit_vpriv[X], hitp->hit_vpriv[Y], 0. );
00640 VUNITIZE( Ru );
00641 VSCALE( Radius, Ru, eto->eto_r );
00642
00643
00644 cv = VDOT( eto->eto_C, Nu );
00645 ch = sqrt( VDOT( eto->eto_C, eto->eto_C ) - cv * cv );
00646
00647 phi = acos( cv / MAGNITUDE(eto->eto_C) );
00648 dv = eto->eto_rd * sin(phi);
00649 dh = -eto->eto_rd * cos(phi);
00650
00651
00652 VCOMB2( Cp, ch, Ru, cv, Nu );
00653 VCOMB2( Dp, dh, Ru, dv, Nu );
00654 VUNITIZE( Cp );
00655 VUNITIZE( Dp );
00656
00657
00658 VSUB2( Hit_Ell, hitp->hit_vpriv, Radius );
00659 xp = VDOT( Hit_Ell, Dp );
00660
00661
00662
00663
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
00671 k_circ = -1. / MAGNITUDE(Radius);
00672
00673 if (fabs(k_ell) < fabs(k_circ)) {
00674
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
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
00699 VSET( Ru, hitp->hit_vpriv[X], hitp->hit_vpriv[Y], 0. );
00700 VUNITIZE( Ru );
00701 VSCALE( Radius, Ru, eto->eto_r );
00702
00703 VSUB2( Hit_Ell, hitp->hit_vpriv, Radius );
00704
00705
00706
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);
00712
00713
00714 theta_v = atan2(hitp->hit_vpriv[Y], hitp->hit_vpriv[X]);
00715
00716
00717 if (theta_u < 0.)
00718 theta_u += bn_twopi;
00719 if (theta_v < 0.)
00720 theta_v += bn_twopi;
00721
00722
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
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
00748
00749
00750
00751
00752
00753
00754
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;
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;
00765 point_t Ell_V;
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
00784 if( ttol->rel <= 0.0 || ttol->rel >= 1.0 ) {
00785 dtol = 0.0;
00786 } else {
00787
00788
00789
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
00799 if (tip->eto_r < b)
00800 dtol = 2 * 0.10 * tip->eto_r;
00801 else
00802 dtol = 2 * 0.10 * b;
00803 } else {
00804
00805 }
00806 } else {
00807
00808 if( ttol->rel <= 0.0 || dtol > ttol->abs )
00809 dtol = ttol->abs;
00810 }
00811
00812 if( ttol->norm > 0.0 )
00813 ntol = ttol->norm;
00814 else
00815
00816 ntol = bn_pi;
00817
00818
00819 ell = rt_mk_ell( &npts, a, b, dtol, ntol );
00820
00821 VMOVE( Nu, tip->eto_N );
00822 VUNITIZE( Nu );
00823 bn_vec_ortho( Bu, Nu );
00824 VUNITIZE( Bu );
00825 VCROSS( Au, Nu, Bu );
00826
00827
00828 nells = rt_num_circular_segments( dtol, tip->eto_r );
00829 theta = bn_twopi / nells;
00830
00831 cv = VDOT( tip->eto_C, Nu );
00832 ch = sqrt( VDOT( tip->eto_C, tip->eto_C ) - cv * cv );
00833
00834 phi = acos( cv / MAGNITUDE(tip->eto_C) );
00835 dv = tip->eto_rd * sin(phi);
00836 dh = -tip->eto_rd * cos(phi);
00837
00838
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
00845 eto_ells = (fastf_t *)bu_malloc(nells * npts * sizeof(point_t), "ells[]");
00846
00847
00848 for (i = 0, ang = 0.; i < nells; i++, ang += theta) {
00849
00850 VCOMB2( Xu, cos(ang), Bu, sin(ang), Au );
00851 VUNITIZE( Xu );
00852
00853 VJOIN1( Ell_V, tip->eto_V, tip->eto_r, Xu );
00854
00855 VCOMB2( Cp, ch, Xu, cv, Nu );
00856 VCOMB2( Dp, dh, Xu, dv, Nu );
00857 VUNITIZE( Cp );
00858 VUNITIZE( Dp );
00859
00860
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
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
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
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
00892
00893
00894
00895
00896
00897
00898
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
00910 VMOVE( p0, pts->p );
00911 VMOVE( p1, pts->next->p );
00912
00913 m = ( p1[X] - p0[X] ) / ( p1[Y] - p0[Y] );
00914 intr = p0[X] - m * p0[Y];
00915
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
00920 dist = fabs( m * mpt[Y] - mpt[X] + intr ) / sqrt( m * m + 1 );
00921
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
00931 if ( dist > dtol || theta0 > ntol || theta1 > ntol ) {
00932
00933 new = rt_ptalloc();
00934 VMOVE( new->p, mpt );
00935 new->next = pts->next;
00936 pts->next = new;
00937
00938 n = 1;
00939
00940 n += rt_ell4( pts, a, b, dtol, ntol );
00941
00942 n += rt_ell4( new, a, b, dtol, ntol );
00943 } else
00944 n = 0;
00945 return( n );
00946 }
00947
00948
00949
00950
00951
00952
00953
00954
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
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
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
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
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
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;
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;
01014 point_t Ell_V;
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;
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
01040 if( ttol->rel <= 0.0 || ttol->rel >= 1.0 ) {
01041 dtol = 0.0;
01042 } else {
01043
01044
01045
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
01055 if (tip->eto_r < b)
01056 dtol = 2 * 0.10 * tip->eto_r;
01057 else
01058 dtol = 2 * 0.10 * b;
01059 } else {
01060
01061 }
01062 } else {
01063
01064 if( ttol->rel <= 0.0 || dtol > ttol->abs )
01065 dtol = ttol->abs;
01066 }
01067
01068 if( ttol->norm > 0.0 )
01069 ntol = ttol->norm;
01070 else
01071
01072 ntol = bn_pi;
01073
01074
01075 ell = rt_mk_ell( &npts, a, b, dtol, ntol );
01076
01077 VMOVE( Nu, tip->eto_N );
01078 VUNITIZE( Nu );
01079 bn_vec_ortho( Bu, Nu );
01080 VUNITIZE( Bu );
01081 VCROSS( Au, Nu, Bu );
01082
01083
01084 nells = rt_num_circular_segments( dtol, tip->eto_r );
01085 theta = bn_twopi / nells;
01086
01087 cv = VDOT( tip->eto_C, Nu );
01088 ch = sqrt( VDOT( tip->eto_C, tip->eto_C ) - cv * cv );
01089
01090 phi = acos( cv / MAGNITUDE(tip->eto_C) );
01091 dv = tip->eto_rd * sin(phi);
01092 dh = -tip->eto_rd * cos(phi);
01093
01094
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
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
01106 for (i = 0, ang = 0.; i < nells; i++, ang += theta) {
01107
01108 VCOMB2( Xu, cos(ang), Bu, sin(ang), Au );
01109 VUNITIZE( Xu );
01110
01111 VJOIN1( Ell_V, tip->eto_V, tip->eto_r, Xu );
01112
01113 VCOMB2( Cp, ch, Xu, cv, Nu );
01114 VCOMB2( Dp, dh, Xu, dv, Nu );
01115 VUNITIZE( Cp );
01116 VUNITIZE( Dp );
01117
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 );
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
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
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
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
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
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
01211
01212
01213
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
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
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);
01250 }
01251
01252
01253
01254
01255
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
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
01301
01302
01303
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
01325 ntohd( (unsigned char *)vec, ep->ext_buf, 11 );
01326
01327
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);
01340 }
01341
01342
01343
01344
01345
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
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
01383 htond( ep->ext_buf, (unsigned char *)vec, 11 );
01384
01385 return(0);
01386 }
01387
01388
01389
01390
01391
01392
01393
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
01435
01436
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;
01449 }
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460