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
00039
00040 #include "globals.hh"
00041
00042 #include "G4EllipticalCone.hh"
00043
00044 #include "G4ClippablePolygon.hh"
00045 #include "G4SolidExtentList.hh"
00046 #include "G4VoxelLimits.hh"
00047 #include "G4AffineTransform.hh"
00048 #include "G4GeometryTolerance.hh"
00049
00050 #include "meshdefs.hh"
00051
00052 #include "Randomize.hh"
00053
00054 #include "G4VGraphicsScene.hh"
00055 #include "G4Polyhedron.hh"
00056 #include "G4NURBS.hh"
00057 #include "G4NURBSbox.hh"
00058 #include "G4VisExtent.hh"
00059
00060
00061
00062 using namespace CLHEP;
00063
00065
00066
00067
00068 G4EllipticalCone::G4EllipticalCone(const G4String& pName,
00069 G4double pxSemiAxis,
00070 G4double pySemiAxis,
00071 G4double pzMax,
00072 G4double pzTopCut)
00073 : G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.),
00074 zTopCut(0.)
00075 {
00076
00077 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00078
00079
00080
00081 if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.) || (pzMax <= 0.) )
00082 {
00083 std::ostringstream message;
00084 message << "Invalid semi-axis or height - " << GetName();
00085 G4Exception("G4EllipticalCone::G4EllipticalCone()", "GeomSolids0002",
00086 FatalErrorInArgument, message);
00087 }
00088 if ( pzTopCut <= 0 )
00089 {
00090 std::ostringstream message;
00091 message << "Invalid z-coordinate for cutting plane - " << GetName();
00092 G4Exception("G4EllipticalCone::G4EllipticalCone()", "InvalidSetup",
00093 FatalErrorInArgument, message);
00094 }
00095
00096 SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax );
00097 SetZCut(pzTopCut);
00098 }
00099
00101
00102
00103
00104
00105 G4EllipticalCone::G4EllipticalCone( __void__& a )
00106 : G4VSolid(a), fpPolyhedron(0), kRadTolerance(0.), fCubicVolume(0.),
00107 fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zheight(0.),
00108 semiAxisMax(0.), zTopCut(0.)
00109 {
00110 }
00111
00113
00114
00115
00116 G4EllipticalCone::~G4EllipticalCone()
00117 {
00118 }
00119
00121
00122
00123
00124 G4EllipticalCone::G4EllipticalCone(const G4EllipticalCone& rhs)
00125 : G4VSolid(rhs),
00126 fpPolyhedron(0), kRadTolerance(rhs.kRadTolerance),
00127 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
00128 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis), zheight(rhs.zheight),
00129 semiAxisMax(rhs.semiAxisMax), zTopCut(rhs.zTopCut)
00130 {
00131 }
00132
00134
00135
00136
00137 G4EllipticalCone& G4EllipticalCone::operator = (const G4EllipticalCone& rhs)
00138 {
00139
00140
00141 if (this == &rhs) { return *this; }
00142
00143
00144
00145 G4VSolid::operator=(rhs);
00146
00147
00148
00149 fpPolyhedron = 0; kRadTolerance = rhs.kRadTolerance;
00150 fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
00151 xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
00152 zheight = rhs.zheight; semiAxisMax = rhs.semiAxisMax; zTopCut = rhs.zTopCut;
00153
00154 return *this;
00155 }
00156
00158
00159
00160
00161 G4bool
00162 G4EllipticalCone::CalculateExtent( const EAxis axis,
00163 const G4VoxelLimits &voxelLimit,
00164 const G4AffineTransform &transform,
00165 G4double &min, G4double &max ) const
00166 {
00167 G4SolidExtentList extentList( axis, voxelLimit );
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177 G4int numPhi = kMaxMeshSections;
00178 G4double sigPhi = twopi/numPhi;
00179
00180
00181
00182
00183
00184
00185
00186 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
00187 G4double dxFudgeBot = xSemiAxis*2.*zheight*rFudge,
00188 dyFudgeBot = ySemiAxis*2.*zheight*rFudge;
00189 G4double dxFudgeTop = xSemiAxis*(zheight-zTopCut)*rFudge,
00190 dyFudgeTop = ySemiAxis*(zheight-zTopCut)*rFudge;
00191
00192
00193
00194
00195
00196
00197 G4ClippablePolygon endPoly1, endPoly2, phiPoly;
00198
00199 G4double phi = 0,
00200 cosPhi = std::cos(phi),
00201 sinPhi = std::sin(phi);
00202 G4ThreeVector v0( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut ),
00203 v1( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut ),
00204 w0, w1;
00205 transform.ApplyPointTransform( v0 );
00206 transform.ApplyPointTransform( v1 );
00207 do
00208 {
00209 phi += sigPhi;
00210 if (numPhi == 1) phi = 0;
00211 cosPhi = std::cos(phi),
00212 sinPhi = std::sin(phi);
00213
00214 w0 = G4ThreeVector( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut );
00215 w1 = G4ThreeVector( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut );
00216 transform.ApplyPointTransform( w0 );
00217 transform.ApplyPointTransform( w1 );
00218
00219
00220
00221
00222 endPoly1.AddVertexInOrder( v0 );
00223 endPoly2.AddVertexInOrder( v1 );
00224
00225
00226
00227
00228 phiPoly.ClearAllVertices();
00229
00230 phiPoly.AddVertexInOrder( v0 );
00231 phiPoly.AddVertexInOrder( v1 );
00232 phiPoly.AddVertexInOrder( w1 );
00233 phiPoly.AddVertexInOrder( w0 );
00234
00235 if (phiPoly.PartialClip( voxelLimit, axis ))
00236 {
00237
00238
00239
00240 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
00241
00242 extentList.AddSurface( phiPoly );
00243 }
00244
00245
00246
00247
00248 v0 = w0;
00249 v1 = w1;
00250 } while( --numPhi > 0 );
00251
00252
00253
00254
00255 if (endPoly1.PartialClip( voxelLimit, axis ))
00256 {
00257 static const G4ThreeVector normal(0,0,+1);
00258 endPoly1.SetNormal( transform.TransformAxis(normal) );
00259 extentList.AddSurface( endPoly1 );
00260 }
00261
00262 if (endPoly2.PartialClip( voxelLimit, axis ))
00263 {
00264 static const G4ThreeVector normal(0,0,-1);
00265 endPoly2.SetNormal( transform.TransformAxis(normal) );
00266 extentList.AddSurface( endPoly2 );
00267 }
00268
00269
00270
00271
00272 return extentList.GetExtent( min, max );
00273 }
00274
00276
00277
00278
00279
00280
00281 EInside G4EllipticalCone::Inside(const G4ThreeVector& p) const
00282 {
00283 G4double rad2oo,
00284 rad2oi;
00285
00286 EInside in;
00287
00288 static const G4double halfRadTol = 0.5*kRadTolerance;
00289 static const G4double halfCarTol = 0.5*kCarTolerance;
00290
00291
00292
00293
00294 if ( (p.z() < -zTopCut - halfCarTol)
00295 || (p.z() > zTopCut + halfCarTol ) )
00296 {
00297 return in = kOutside;
00298 }
00299
00300 rad2oo= sqr(p.x()/( xSemiAxis + halfRadTol ))
00301 + sqr(p.y()/( ySemiAxis + halfRadTol ));
00302
00303 if ( rad2oo > sqr( zheight-p.z() ) )
00304 {
00305 return in = kOutside;
00306 }
00307
00308
00309
00310 rad2oi = sqr(p.x()/( xSemiAxis - halfRadTol ))
00311 + sqr(p.y()/( ySemiAxis - halfRadTol ));
00312
00313 if (rad2oi < sqr( zheight-p.z() ) )
00314 {
00315 in = ( ( p.z() < -zTopCut + halfRadTol )
00316 || ( p.z() > zTopCut - halfRadTol ) ) ? kSurface : kInside;
00317 }
00318 else
00319 {
00320 in = kSurface;
00321 }
00322
00323 return in;
00324 }
00325
00327
00328
00329
00330 G4ThreeVector G4EllipticalCone::SurfaceNormal( const G4ThreeVector& p) const
00331 {
00332
00333 G4double rx = sqr(p.x()/xSemiAxis),
00334 ry = sqr(p.y()/ySemiAxis);
00335
00336 G4double rds = std::sqrt(rx + ry);
00337
00338 G4ThreeVector norm;
00339
00340 if( (p.z() < -zTopCut) && ((rx+ry) < sqr(zTopCut + zheight)) )
00341 {
00342 return G4ThreeVector( 0., 0., -1. );
00343 }
00344
00345 if( (p.z() > (zheight > zTopCut ? zheight : zTopCut)) &&
00346 ((rx+ry) < sqr(zheight-zTopCut)) )
00347 {
00348 return G4ThreeVector( 0., 0., 1. );
00349 }
00350
00351 if( p.z() > rds + 2.*zTopCut - zheight )
00352 {
00353 if ( p.z() > zTopCut )
00354 {
00355 if( p.x() == 0. )
00356 {
00357 norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., 1. );
00358 return norm /= norm.mag();
00359 }
00360 if( p.y() == 0. )
00361 {
00362 norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., 1. );
00363 return norm /= norm.mag();
00364 }
00365
00366 G4double k = std::fabs(p.x()/p.y());
00367 G4double c2 = sqr(zheight-zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis));
00368 G4double x = std::sqrt(c2);
00369 G4double y = k*x;
00370
00371 x /= sqr(xSemiAxis);
00372 y /= sqr(ySemiAxis);
00373
00374 norm = G4ThreeVector( p.x() < 0. ? -x : x,
00375 p.y() < 0. ? -y : y,
00376 - ( zheight - zTopCut ) );
00377 norm /= norm.mag();
00378 norm += G4ThreeVector( 0., 0., 1. );
00379 return norm /= norm.mag();
00380 }
00381
00382 return G4ThreeVector( 0., 0., 1. );
00383 }
00384
00385 if( p.z() < rds - 2.*zTopCut - zheight )
00386 {
00387 if( p.x() == 0. )
00388 {
00389 norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., -1. );
00390 return norm /= norm.mag();
00391 }
00392 if( p.y() == 0. )
00393 {
00394 norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., -1. );
00395 return norm /= norm.mag();
00396 }
00397
00398 G4double k = std::fabs(p.x()/p.y());
00399 G4double c2 = sqr(zheight+zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis));
00400 G4double x = std::sqrt(c2);
00401 G4double y = k*x;
00402
00403 x /= sqr(xSemiAxis);
00404 y /= sqr(ySemiAxis);
00405
00406 norm = G4ThreeVector( p.x() < 0. ? -x : x,
00407 p.y() < 0. ? -y : y,
00408 - ( zheight - zTopCut ) );
00409 norm /= norm.mag();
00410 norm += G4ThreeVector( 0., 0., -1. );
00411 return norm /= norm.mag();
00412 }
00413
00414 norm = G4ThreeVector(p.x()/sqr(xSemiAxis), p.y()/sqr(ySemiAxis), rds);
00415
00416 G4double k = std::tan(pi/8.);
00417 G4double c = -zTopCut - k*(zTopCut + zheight);
00418
00419 if( p.z() < -k*rds + c )
00420 return G4ThreeVector (0.,0.,-1.);
00421
00422 return norm /= norm.mag();
00423 }
00424
00426
00427
00428
00429
00430 G4double G4EllipticalCone::DistanceToIn( const G4ThreeVector& p,
00431 const G4ThreeVector& v ) const
00432 {
00433 static const G4double halfTol = 0.5*kCarTolerance;
00434
00435 G4double distMin = kInfinity;
00436
00437
00438
00439 G4double sigz = p.z()+zTopCut;
00440
00441
00442
00443
00444
00445 if (sigz < halfTol)
00446 {
00447
00448
00449
00450
00451 if (v.z() <= 0)
00452 {
00453
00454
00455
00456
00457 if (sigz < 0) return kInfinity;
00458
00459
00460
00461
00462
00463
00464 if ( sqr(p.x()/( xSemiAxis - halfTol ))
00465 + sqr(p.y()/( ySemiAxis - halfTol )) <= sqr( zheight+zTopCut ) )
00466 return kInfinity;
00467
00468 }
00469 else
00470 {
00471
00472
00473
00474 G4double q = -sigz/v.z();
00475
00476
00477
00478
00479 G4double xi = p.x() + q*v.x(),
00480 yi = p.y() + q*v.y();
00481
00482
00483
00484
00485 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight + zTopCut ) )
00486 {
00487
00488
00489
00490 return (sigz < -halfTol) ? q : 0;
00491 }
00492 else if (xi/(xSemiAxis*xSemiAxis)*v.x()
00493 + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
00494 {
00495
00496
00497
00498
00499
00500 }
00501 }
00502 }
00503
00504
00505
00506
00507 sigz = p.z() - zTopCut;
00508
00509 if (sigz > -halfTol)
00510 {
00511 if (v.z() >= 0)
00512 {
00513
00514 if (sigz > 0) return kInfinity;
00515
00516 if ( sqr(p.x()/( xSemiAxis - halfTol ))
00517 + sqr(p.y()/( ySemiAxis - halfTol )) <= sqr( zheight-zTopCut ) )
00518 return kInfinity;
00519
00520 }
00521 else {
00522 G4double q = -sigz/v.z();
00523
00524 G4double xi = p.x() + q*v.x(),
00525 yi = p.y() + q*v.y();
00526
00527 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight - zTopCut ) )
00528 {
00529 return (sigz > -halfTol) ? q : 0;
00530 }
00531 else if (xi/(xSemiAxis*xSemiAxis)*v.x()
00532 + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
00533 {
00534
00535 }
00536 }
00537 }
00538
00539
00540 #if 0
00541
00542
00543
00544 if (p.z() < -zTopCut - 0.5*kCarTolerance)
00545 {
00546 if (v.z() <= 0.0)
00547 return distMin;
00548
00549 G4double lambda = (-zTopCut - p.z())/v.z();
00550
00551 if ( sqr((lambda*v.x()+p.x())/xSemiAxis) +
00552 sqr((lambda*v.y()+p.y())/ySemiAxis) <=
00553 sqr(zTopCut + zheight + 0.5*kRadTolerance) )
00554 {
00555 return distMin = std::fabs(lambda);
00556 }
00557 }
00558
00559 if (p.z() > zTopCut+0.5*kCarTolerance)
00560 {
00561 if (v.z() >= 0.0)
00562 { return distMin; }
00563
00564 G4double lambda = (zTopCut - p.z()) / v.z();
00565
00566 if ( sqr((lambda*v.x() + p.x())/xSemiAxis) +
00567 sqr((lambda*v.y() + p.y())/ySemiAxis) <=
00568 sqr(zheight - zTopCut + 0.5*kRadTolerance) )
00569 {
00570 return distMin = std::fabs(lambda);
00571 }
00572 }
00573
00574 if (p.z() > zTopCut - halfTol
00575 && p.z() < zTopCut + halfTol )
00576 {
00577 if (v.z() > 0.)
00578 { return kInfinity; }
00579
00580 return distMin = 0.;
00581 }
00582
00583 if (p.z() < -zTopCut + halfTol
00584 && p.z() > -zTopCut - halfTol)
00585 {
00586 if (v.z() < 0.)
00587 { return distMin = kInfinity; }
00588
00589 return distMin = 0.;
00590 }
00591
00592 #endif
00593
00594
00595
00596
00597 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
00598 G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) +
00599 v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
00600 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) -
00601 sqr(zheight - p.z());
00602
00603 G4double discr = B*B - 4.*A*C;
00604
00605
00606
00607 if ( discr < -halfTol )
00608 { return distMin; }
00609
00610
00611
00612 if ( (discr >= - halfTol ) && (discr < halfTol ) )
00613 {
00614 return distMin = std::fabs(-B/(2.*A));
00615 }
00616
00617 G4double plus = (-B+std::sqrt(discr))/(2.*A);
00618 G4double minus = (-B-std::sqrt(discr))/(2.*A);
00619
00620
00621
00622 if ( ( std::fabs(plus) < halfTol )||( std::fabs(minus) < halfTol ) )
00623 {
00624 G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
00625 p.y()/(ySemiAxis*ySemiAxis),
00626 -( p.z() - zheight ));
00627 if ( truenorm*v >= 0)
00628 {
00629 return kInfinity;
00630 }
00631 else
00632 {
00633 return 0;
00634 }
00635 }
00636
00637
00638 G4double lambda = 0;
00639
00640 if ( minus > halfTol && minus < distMin )
00641 {
00642 lambda = minus ;
00643
00644 G4ThreeVector pin = p + lambda*v;
00645 if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
00646 {
00647 G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
00648 pin.y()/(ySemiAxis*ySemiAxis),
00649 - ( pin.z() - zheight ));
00650 if ( truenorm*v < 0)
00651 {
00652 distMin = lambda;
00653 }
00654 }
00655 }
00656 if ( plus > halfTol && plus < distMin )
00657 {
00658 lambda = plus ;
00659
00660 G4ThreeVector pin = p + lambda*v;
00661 if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
00662 {
00663 G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
00664 pin.y()/(ySemiAxis*ySemiAxis),
00665 - ( pin.z() - zheight ) );
00666 if ( truenorm*v < 0)
00667 {
00668 distMin = lambda;
00669 }
00670 }
00671 }
00672 if (distMin < halfTol) distMin=0.;
00673 return distMin ;
00674 }
00675
00677
00678
00679
00680
00681 G4double G4EllipticalCone::DistanceToIn(const G4ThreeVector& p) const
00682 {
00683 G4double distR, distR2, distZ, maxDim;
00684 G4double distRad;
00685
00686
00687
00688
00689 if( (p.z() <= -zTopCut) && (sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
00690 <= sqr(zTopCut + zheight + 0.5*kCarTolerance )) )
00691 {
00692
00693 return distZ = std::fabs(zTopCut + p.z());
00694 }
00695
00696 if( (p.z() >= zTopCut) && (sqr(p.x()/xSemiAxis)+sqr(p.y()/ySemiAxis)
00697 <= sqr(zheight - zTopCut + kCarTolerance/2.0 )) )
00698 {
00699 return distZ = std::fabs(p.z() - zTopCut);
00700 }
00701
00702
00703
00704
00705
00706 maxDim = xSemiAxis >= ySemiAxis ? xSemiAxis:ySemiAxis;
00707 distRad = std::sqrt(p.x()*p.x()+p.y()*p.y());
00708
00709 if( p.z() > maxDim*distRad + zTopCut*(1.+maxDim)-sqr(maxDim)*zheight )
00710 {
00711 distR2 = sqr(p.z() - zTopCut) + sqr(distRad - maxDim*(zheight - zTopCut));
00712 return std::sqrt( distR2 );
00713 }
00714
00715 if( distRad > maxDim*( zheight - p.z() ) )
00716 {
00717 if( p.z() > maxDim*distRad - (zTopCut*(1.+maxDim)+sqr(maxDim)*zheight) )
00718 {
00719 G4double zVal = (p.z()-maxDim*(distRad-maxDim*zheight))/(1.+sqr(maxDim));
00720 G4double rVal = maxDim*(zheight - zVal);
00721 return distR = std::sqrt(sqr(p.z() - zVal) + sqr(distRad - rVal));
00722 }
00723 }
00724
00725 if( distRad <= maxDim*(zheight - p.z()) )
00726 {
00727 distR2 = sqr(distRad - maxDim*(zheight + zTopCut)) + sqr(p.z() + zTopCut);
00728 return std::sqrt( distR2 );
00729 }
00730
00731 return distR = 0;
00732 }
00733
00735
00736
00737
00738
00739 G4double G4EllipticalCone::DistanceToOut(const G4ThreeVector& p,
00740 const G4ThreeVector& v,
00741 const G4bool calcNorm,
00742 G4bool *validNorm,
00743 G4ThreeVector *n ) const
00744 {
00745 G4double distMin, lambda;
00746 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
00747
00748 distMin = kInfinity;
00749 surface = kNoSurf;
00750
00751 if (v.z() < 0.0)
00752 {
00753 lambda = (-p.z() - zTopCut)/v.z();
00754
00755 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) +
00756 sqr((p.y() + lambda*v.y())/ySemiAxis)) <
00757 sqr(zheight + zTopCut + 0.5*kCarTolerance) )
00758 {
00759 distMin = std::fabs(lambda);
00760
00761 if (!calcNorm) { return distMin; }
00762 }
00763 distMin = std::fabs(lambda);
00764 surface = kPlaneSurf;
00765 }
00766
00767 if (v.z() > 0.0)
00768 {
00769 lambda = (zTopCut - p.z()) / v.z();
00770
00771 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis)
00772 + sqr((p.y() + lambda*v.y())/ySemiAxis) )
00773 < (sqr(zheight - zTopCut + 0.5*kCarTolerance)) )
00774 {
00775 distMin = std::fabs(lambda);
00776 if (!calcNorm) { return distMin; }
00777 }
00778 distMin = std::fabs(lambda);
00779 surface = kPlaneSurf;
00780 }
00781
00782
00783
00784
00785 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
00786 G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) +
00787 v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
00788 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
00789 - sqr(zheight - p.z());
00790
00791 G4double discr = B*B - 4.*A*C;
00792
00793 if ( discr >= - 0.5*kCarTolerance && discr < 0.5*kCarTolerance )
00794 {
00795 if(!calcNorm) { return distMin = std::fabs(-B/(2.*A)); }
00796 }
00797
00798 else if ( discr > 0.5*kCarTolerance )
00799 {
00800 G4double plus = (-B+std::sqrt(discr))/(2.*A);
00801 G4double minus = (-B-std::sqrt(discr))/(2.*A);
00802
00803 if ( plus > 0.5*kCarTolerance && minus > 0.5*kCarTolerance )
00804 {
00805
00806
00807 lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;
00808 }
00809 else
00810 {
00811
00812
00813
00814 lambda = plus > -0.5*kCarTolerance ? plus : 0;
00815 }
00816
00817 if ( std::fabs(lambda) < distMin )
00818 {
00819 if( std::fabs(lambda) > 0.5*kCarTolerance)
00820 {
00821 distMin = std::fabs(lambda);
00822 surface = kCurvedSurf;
00823 }
00824 else
00825 {
00826 G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
00827 p.y()/(ySemiAxis*ySemiAxis),
00828 -( p.z() - zheight ));
00829 if( truenorm.dot(v) > 0 )
00830 {
00831 distMin = 0.0;
00832 surface = kCurvedSurf;
00833 }
00834 }
00835 }
00836 }
00837
00838
00839
00840 if (calcNorm)
00841 {
00842 if (surface == kNoSurf)
00843 {
00844 *validNorm = false;
00845 }
00846 else
00847 {
00848 *validNorm = true;
00849 switch (surface)
00850 {
00851 case kPlaneSurf:
00852 {
00853 *n = G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
00854 }
00855 break;
00856
00857 case kCurvedSurf:
00858 {
00859 G4ThreeVector pexit = p + distMin*v;
00860 G4ThreeVector truenorm( pexit.x()/(xSemiAxis*xSemiAxis),
00861 pexit.y()/(ySemiAxis*ySemiAxis),
00862 -( pexit.z() - zheight ) );
00863 truenorm /= truenorm.mag();
00864 *n= truenorm;
00865 }
00866 break;
00867
00868 default:
00869 DumpInfo();
00870 std::ostringstream message;
00871 G4int oldprc = message.precision(16);
00872 message << "Undefined side for valid surface normal to solid."
00873 << G4endl
00874 << "Position:" << G4endl
00875 << " p.x() = " << p.x()/mm << " mm" << G4endl
00876 << " p.y() = " << p.y()/mm << " mm" << G4endl
00877 << " p.z() = " << p.z()/mm << " mm" << G4endl
00878 << "Direction:" << G4endl
00879 << " v.x() = " << v.x() << G4endl
00880 << " v.y() = " << v.y() << G4endl
00881 << " v.z() = " << v.z() << G4endl
00882 << "Proposed distance :" << G4endl
00883 << " distMin = " << distMin/mm << " mm";
00884 message.precision(oldprc);
00885 G4Exception("G4EllipticalCone::DistanceToOut(p,v,..)",
00886 "GeomSolids1002", JustWarning, message);
00887 break;
00888 }
00889 }
00890 }
00891
00892 if (distMin<0.5*kCarTolerance) { distMin=0; }
00893
00894 return distMin;
00895 }
00896
00898
00899
00900
00901 G4double G4EllipticalCone::DistanceToOut(const G4ThreeVector& p) const
00902 {
00903 G4double rds,roo,roo1, distR, distZ, distMin=0.;
00904 G4double minAxis = xSemiAxis < ySemiAxis ? xSemiAxis : ySemiAxis;
00905
00906 #ifdef G4SPECSDEBUG
00907 if( Inside(p) == kOutside )
00908 {
00909 DumpInfo();
00910 std::ostringstream message;
00911 G4int oldprc = message.precision(16);
00912 message << "Point p is outside !?" << G4endl
00913 << "Position:" << G4endl
00914 << " p.x() = " << p.x()/mm << " mm" << G4endl
00915 << " p.y() = " << p.y()/mm << " mm" << G4endl
00916 << " p.z() = " << p.z()/mm << " mm";
00917 message.precision(oldprc) ;
00918 G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
00919 JustWarning, message);
00920 }
00921 #endif
00922
00923
00924
00925
00926
00927 if( sqr(p.x()/minAxis)+sqr(p.y()/minAxis) < sqr(zheight - p.z()) )
00928 {
00929 rds = std::sqrt(sqr(p.x()) + sqr(p.y()));
00930 roo = minAxis*(zheight-p.z());
00931 roo1 = minAxis*(zheight-zTopCut);
00932
00933 distZ=zTopCut - std::fabs(p.z()) ;
00934 distR=(roo-rds)/(std::sqrt(1+sqr(minAxis)));
00935
00936 if(rds>roo1)
00937 {
00938 distMin=(zTopCut-p.z())*(roo-rds)/(roo-roo1);
00939 distMin=std::min(distMin,distR);
00940 }
00941 distMin=std::min(distR,distZ);
00942 }
00943
00944 return distMin;
00945 }
00946
00948
00949
00950
00951 G4GeometryType G4EllipticalCone::GetEntityType() const
00952 {
00953 return G4String("G4EllipticalCone");
00954 }
00955
00957
00958
00959
00960 G4VSolid* G4EllipticalCone::Clone() const
00961 {
00962 return new G4EllipticalCone(*this);
00963 }
00964
00966
00967
00968
00969 std::ostream& G4EllipticalCone::StreamInfo( std::ostream& os ) const
00970 {
00971 G4int oldprc = os.precision(16);
00972 os << "-----------------------------------------------------------\n"
00973 << " *** Dump for solid - " << GetName() << " ***\n"
00974 << " ===================================================\n"
00975 << " Solid type: G4EllipticalCone\n"
00976 << " Parameters: \n"
00977
00978 << " semi-axis x: " << xSemiAxis/mm << " mm \n"
00979 << " semi-axis y: " << ySemiAxis/mm << " mm \n"
00980 << " height z: " << zheight/mm << " mm \n"
00981 << " half length in z: " << zTopCut/mm << " mm \n"
00982 << "-----------------------------------------------------------\n";
00983 os.precision(oldprc);
00984
00985 return os;
00986 }
00987
00989
00990
00991
00992
00993
00994 G4ThreeVector G4EllipticalCone::GetPointOnSurface() const
00995 {
00996
00997 G4double phi, sinphi, cosphi, aOne, aTwo, aThree,
00998 chose, zRand, rRand1, rRand2;
00999
01000 G4double rOne = std::sqrt(sqr(xSemiAxis)
01001 + sqr(ySemiAxis))*(zheight - zTopCut);
01002 G4double rTwo = std::sqrt(sqr(xSemiAxis)
01003 + sqr(ySemiAxis))*(zheight + zTopCut);
01004
01005 aOne = pi*(rOne + rTwo)*std::sqrt(sqr(rOne - rTwo)+sqr(2.*zTopCut));
01006 aTwo = pi*xSemiAxis*ySemiAxis*sqr(zheight+zTopCut);
01007 aThree = pi*xSemiAxis*ySemiAxis*sqr(zheight-zTopCut);
01008
01009 phi = RandFlat::shoot(0.,twopi);
01010 cosphi = std::cos(phi);
01011 sinphi = std::sin(phi);
01012
01013 if(zTopCut >= zheight) aThree = 0.;
01014
01015 chose = RandFlat::shoot(0.,aOne+aTwo+aThree);
01016 if((chose>=0.) && (chose<aOne))
01017 {
01018 zRand = RandFlat::shoot(-zTopCut,zTopCut);
01019 return G4ThreeVector(xSemiAxis*(zheight-zRand)*cosphi,
01020 ySemiAxis*(zheight-zRand)*sinphi,zRand);
01021 }
01022 else if((chose>=aOne) && (chose<aOne+aTwo))
01023 {
01024 do
01025 {
01026 rRand1 = RandFlat::shoot(0.,1.) ;
01027 rRand2 = RandFlat::shoot(0.,1.) ;
01028 } while ( rRand2 >= rRand1 ) ;
01029
01030
01031 return G4ThreeVector(rRand1*xSemiAxis*(zheight+zTopCut)*cosphi,
01032 rRand1*ySemiAxis*(zheight+zTopCut)*sinphi, -zTopCut);
01033
01034 }
01035
01036
01037
01038 do
01039 {
01040 rRand1 = RandFlat::shoot(0.,1.) ;
01041 rRand2 = RandFlat::shoot(0.,1.) ;
01042 } while ( rRand2 >= rRand1 ) ;
01043
01044 return G4ThreeVector(rRand1*xSemiAxis*(zheight-zTopCut)*cosphi,
01045 rRand1*ySemiAxis*(zheight-zTopCut)*sinphi, zTopCut);
01046 }
01047
01048
01049
01050
01051
01052 void G4EllipticalCone::DescribeYourselfTo (G4VGraphicsScene& scene) const
01053 {
01054 scene.AddSolid(*this);
01055 }
01056
01057 G4VisExtent G4EllipticalCone::GetExtent() const
01058 {
01059
01060
01061 G4double maxDim;
01062 maxDim = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
01063 maxDim = maxDim > zTopCut ? maxDim : zTopCut;
01064
01065 return G4VisExtent (-maxDim, maxDim,
01066 -maxDim, maxDim,
01067 -maxDim, maxDim);
01068 }
01069
01070 G4NURBS* G4EllipticalCone::CreateNURBS () const
01071 {
01072
01073
01074 return new G4NURBSbox(xSemiAxis, ySemiAxis,zheight);
01075 }
01076
01077 G4Polyhedron* G4EllipticalCone::CreatePolyhedron () const
01078 {
01079 return new G4PolyhedronEllipticalCone(xSemiAxis, ySemiAxis, zheight, zTopCut);
01080 }
01081
01082 G4Polyhedron* G4EllipticalCone::GetPolyhedron () const
01083 {
01084 if ( (!fpPolyhedron)
01085 || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01086 fpPolyhedron->GetNumberOfRotationSteps()) )
01087 {
01088 delete fpPolyhedron;
01089 fpPolyhedron = CreatePolyhedron();
01090 }
01091 return fpPolyhedron;
01092 }