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 #include "globals.hh"
00037
00038 #include "G4Paraboloid.hh"
00039
00040 #include "G4VoxelLimits.hh"
00041 #include "G4AffineTransform.hh"
00042
00043 #include "meshdefs.hh"
00044
00045 #include "Randomize.hh"
00046
00047 #include "G4VGraphicsScene.hh"
00048 #include "G4Polyhedron.hh"
00049 #include "G4NURBS.hh"
00050 #include "G4NURBSbox.hh"
00051 #include "G4VisExtent.hh"
00052
00053 using namespace CLHEP;
00054
00056
00057
00058
00059 G4Paraboloid::G4Paraboloid(const G4String& pName,
00060 G4double pDz,
00061 G4double pR1,
00062 G4double pR2)
00063 : G4VSolid(pName), fpPolyhedron(0), fSurfaceArea(0.), fCubicVolume(0.)
00064
00065 {
00066 if( (pDz <= 0.) || (pR2 <= pR1) || (pR1 < 0.) )
00067 {
00068 std::ostringstream message;
00069 message << "Invalid dimensions. Negative Input Values or R1>=R2 - "
00070 << GetName();
00071 G4Exception("G4Paraboloid::G4Paraboloid()", "GeomSolids0002",
00072 FatalErrorInArgument, message,
00073 "Z half-length must be larger than zero or R1>=R2.");
00074 }
00075
00076 r1 = pR1;
00077 r2 = pR2;
00078 dz = pDz;
00079
00080
00081
00082
00083
00084
00085 k1 = (r2 * r2 - r1 * r1) / 2 / dz;
00086 k2 = (r2 * r2 + r1 * r1) / 2;
00087 }
00088
00090
00091
00092
00093
00094 G4Paraboloid::G4Paraboloid( __void__& a )
00095 : G4VSolid(a), fpPolyhedron(0), fSurfaceArea(0.), fCubicVolume(0.),
00096 dz(0.), r1(0.), r2(0.), k1(0.), k2(0.)
00097 {
00098 }
00099
00101
00102
00103
00104 G4Paraboloid::~G4Paraboloid()
00105 {
00106 }
00107
00109
00110
00111
00112 G4Paraboloid::G4Paraboloid(const G4Paraboloid& rhs)
00113 : G4VSolid(rhs), fpPolyhedron(0),
00114 fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume),
00115 dz(rhs.dz), r1(rhs.r1), r2(rhs.r2), k1(rhs.k1), k2(rhs.k2)
00116 {
00117 }
00118
00119
00121
00122
00123
00124 G4Paraboloid& G4Paraboloid::operator = (const G4Paraboloid& rhs)
00125 {
00126
00127
00128 if (this == &rhs) { return *this; }
00129
00130
00131
00132 G4VSolid::operator=(rhs);
00133
00134
00135
00136 fpPolyhedron = 0;
00137 fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
00138 dz = rhs.dz; r1 = rhs.r1; r2 = rhs.r2; k1 = rhs.k1; k2 = rhs.k2;
00139
00140 return *this;
00141 }
00142
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00157
00158
00159
00160 G4bool
00161 G4Paraboloid::CalculateExtent(const EAxis pAxis,
00162 const G4VoxelLimits& pVoxelLimit,
00163 const G4AffineTransform& pTransform,
00164 G4double& pMin, G4double& pMax) const
00165 {
00166 G4double xMin = -r2 + pTransform.NetTranslation().x(),
00167 xMax = r2 + pTransform.NetTranslation().x(),
00168 yMin = -r2 + pTransform.NetTranslation().y(),
00169 yMax = r2 + pTransform.NetTranslation().y(),
00170 zMin = -dz + pTransform.NetTranslation().z(),
00171 zMax = dz + pTransform.NetTranslation().z();
00172
00173 if(!pTransform.IsRotated()
00174 || pTransform.NetRotation()(G4ThreeVector(0, 0, 1)) == G4ThreeVector(0, 0, 1))
00175 {
00176 if(pVoxelLimit.IsXLimited())
00177 {
00178 if(pVoxelLimit.GetMaxXExtent() < xMin - 0.5 * kCarTolerance
00179 || pVoxelLimit.GetMinXExtent() > xMax + 0.5 * kCarTolerance)
00180 {
00181 return false;
00182 }
00183 else
00184 {
00185 if(pVoxelLimit.GetMinXExtent() > xMin)
00186 {
00187 xMin = pVoxelLimit.GetMinXExtent();
00188 }
00189 if(pVoxelLimit.GetMaxXExtent() < xMax)
00190 {
00191 xMax = pVoxelLimit.GetMaxXExtent();
00192 }
00193 }
00194 }
00195 if(pVoxelLimit.IsYLimited())
00196 {
00197 if(pVoxelLimit.GetMaxYExtent() < yMin - 0.5 * kCarTolerance
00198 || pVoxelLimit.GetMinYExtent() > yMax + 0.5 * kCarTolerance)
00199 {
00200 return false;
00201 }
00202 else
00203 {
00204 if(pVoxelLimit.GetMinYExtent() > yMin)
00205 {
00206 yMin = pVoxelLimit.GetMinYExtent();
00207 }
00208 if(pVoxelLimit.GetMaxYExtent() < yMax)
00209 {
00210 yMax = pVoxelLimit.GetMaxYExtent();
00211 }
00212 }
00213 }
00214 if(pVoxelLimit.IsZLimited())
00215 {
00216 if(pVoxelLimit.GetMaxZExtent() < zMin - 0.5 * kCarTolerance
00217 || pVoxelLimit.GetMinZExtent() > zMax + 0.5 * kCarTolerance)
00218 {
00219 return false;
00220 }
00221 else
00222 {
00223 if(pVoxelLimit.GetMinZExtent() > zMin)
00224 {
00225 zMin = pVoxelLimit.GetMinZExtent();
00226 }
00227 if(pVoxelLimit.GetMaxZExtent() < zMax)
00228 {
00229 zMax = pVoxelLimit.GetMaxZExtent();
00230 }
00231 }
00232 }
00233 switch(pAxis)
00234 {
00235 case kXAxis:
00236 pMin = xMin;
00237 pMax = xMax;
00238 break;
00239 case kYAxis:
00240 pMin = yMin;
00241 pMax = yMax;
00242 break;
00243 case kZAxis:
00244 pMin = zMin;
00245 pMax = zMax;
00246 break;
00247 default:
00248 pMin = 0;
00249 pMax = 0;
00250 return false;
00251 }
00252 }
00253 else
00254 {
00255 G4bool existsAfterClip=true;
00256
00257
00258
00259 G4int noPolygonVertices=0;
00260 G4ThreeVectorList* vertices
00261 = CreateRotatedVertices(pTransform,noPolygonVertices);
00262
00263 if(pAxis == kXAxis || pAxis == kYAxis || pAxis == kZAxis)
00264 {
00265
00266 pMin = kInfinity;
00267 pMax = -kInfinity;
00268
00269 for(G4ThreeVectorList::iterator it = vertices->begin();
00270 it < vertices->end(); it++)
00271 {
00272 if(pMin > (*it)[pAxis]) pMin = (*it)[pAxis];
00273 if((*it)[pAxis] < pVoxelLimit.GetMinExtent(pAxis))
00274 {
00275 pMin = pVoxelLimit.GetMinExtent(pAxis);
00276 }
00277 if(pMax < (*it)[pAxis])
00278 {
00279 pMax = (*it)[pAxis];
00280 }
00281 if((*it)[pAxis] > pVoxelLimit.GetMaxExtent(pAxis))
00282 {
00283 pMax = pVoxelLimit.GetMaxExtent(pAxis);
00284 }
00285 }
00286
00287 if(pMin > pVoxelLimit.GetMaxExtent(pAxis)
00288 || pMax < pVoxelLimit.GetMinExtent(pAxis)) { existsAfterClip = false; }
00289 }
00290 else
00291 {
00292 pMin = 0;
00293 pMax = 0;
00294 existsAfterClip = false;
00295 }
00296 delete vertices;
00297 return existsAfterClip;
00298 }
00299 return true;
00300 }
00301
00303
00304
00305
00306 EInside G4Paraboloid::Inside(const G4ThreeVector& p) const
00307 {
00308
00309
00310 if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance) { return kOutside; }
00311
00312 G4double rho2 = p.perp2(),
00313 rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance),
00314 A = rho2 - ((k1 *p.z() + k2) + 0.25 * kCarTolerance * kCarTolerance);
00315
00316 if(A < 0 && sqr(A) > rhoSurfTimesTol2)
00317 {
00318
00319
00320
00321 if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
00322 {
00323
00324
00325
00326 return kSurface;
00327 }
00328 else
00329 {
00330 return kInside;
00331 }
00332 }
00333 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
00334 {
00335
00336
00337 return kSurface;
00338 }
00339 else
00340 {
00341 return kOutside;
00342 }
00343 }
00344
00346
00347
00348 G4ThreeVector G4Paraboloid::SurfaceNormal( const G4ThreeVector& p) const
00349 {
00350 G4ThreeVector n(0, 0, 0);
00351 if(std::fabs(p.z()) > dz + 0.5 * kCarTolerance)
00352 {
00353
00354
00355 n = G4ThreeVector(0, 0, p.z()/std::fabs(p.z()));
00356 }
00357 else if(std::fabs(p.z()) > dz - 0.5 * kCarTolerance)
00358 {
00359
00360
00361
00362 if(p.z() < 0)
00363 {
00364 if(p.perp2() > sqr(r1 + 0.5 * kCarTolerance))
00365 {
00366 n = G4ThreeVector(p.x(), p.y(), -k1 / 2).unit();
00367 }
00368 else if(r1 < 0.5 * kCarTolerance
00369 || p.perp2() > sqr(r1 - 0.5 * kCarTolerance))
00370 {
00371 n = G4ThreeVector(p.x(), p.y(), 0.).unit()
00372 + G4ThreeVector(0., 0., -1.).unit();
00373 n = n.unit();
00374 }
00375 else
00376 {
00377 n = G4ThreeVector(0., 0., -1.);
00378 }
00379 }
00380 else
00381 {
00382 if(p.perp2() > sqr(r2 + 0.5 * kCarTolerance))
00383 {
00384 n = G4ThreeVector(p.x(), p.y(), 0.).unit();
00385 }
00386 else if(r2 < 0.5 * kCarTolerance
00387 || p.perp2() > sqr(r2 - 0.5 * kCarTolerance))
00388 {
00389 n = G4ThreeVector(p.x(), p.y(), 0.).unit()
00390 + G4ThreeVector(0., 0., 1.).unit();
00391 n = n.unit();
00392 }
00393 else
00394 {
00395 n = G4ThreeVector(0., 0., 1.);
00396 }
00397 }
00398 }
00399 else
00400 {
00401 G4double rho2 = p.perp2();
00402 G4double rhoSurfTimesTol2 = (k1 * p.z() + k2) * sqr(kCarTolerance);
00403 G4double A = rho2 - ((k1 *p.z() + k2)
00404 + 0.25 * kCarTolerance * kCarTolerance);
00405
00406 if(A < 0 && sqr(A) > rhoSurfTimesTol2)
00407 {
00408
00409
00410
00411 if(p.mag2() != 0) { n = p.unit(); }
00412 }
00413 else if(A <= 0 || sqr(A) < rhoSurfTimesTol2)
00414 {
00415
00416
00417 n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
00418 }
00419 else
00420 {
00421 n = G4ThreeVector(p.x(), p.y(), - k1 / 2).unit();
00422 }
00423 }
00424
00425 if(n.mag2() == 0)
00426 {
00427 std::ostringstream message;
00428 message << "No normal defined for this point p." << G4endl
00429 << " p = " << 1 / mm * p << " mm";
00430 G4Exception("G4Paraboloid::SurfaceNormal(p)", "GeomSolids1002",
00431 JustWarning, message);
00432 }
00433 return n;
00434 }
00435
00437
00438
00439
00440
00441
00442 G4double G4Paraboloid::DistanceToIn( const G4ThreeVector& p,
00443 const G4ThreeVector& v ) const
00444 {
00445 G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
00446 G4double tol2 = kCarTolerance*kCarTolerance;
00447 G4double tolh = 0.5*kCarTolerance;
00448
00449 if(r2 && p.z() > - tolh + dz)
00450 {
00451
00452
00453 if(v.z() < 0)
00454 {
00455 G4double intersection = (dz - p.z()) / v.z();
00456 if(sqr(p.x() + v.x()*intersection)
00457 + sqr(p.y() + v.y()*intersection) < sqr(r2 + 0.5 * kCarTolerance))
00458 {
00459 if(p.z() < tolh + dz)
00460 { return 0; }
00461 else
00462 { return intersection; }
00463 }
00464 }
00465 else
00466 {
00467 return kInfinity;
00468 }
00469 }
00470 else if(r1 && p.z() < tolh - dz)
00471 {
00472
00473
00474 if(v.z() > 0)
00475 {
00476 G4double intersection = (-dz - p.z()) / v.z();
00477 if(sqr(p.x() + v.x()*intersection)
00478 + sqr(p.y() + v.y()*intersection) < sqr(r1 + 0.5 * kCarTolerance))
00479 {
00480 if(p.z() > -tolh - dz)
00481 {
00482 return 0;
00483 }
00484 else
00485 {
00486 return intersection;
00487 }
00488 }
00489 }
00490 else
00491 {
00492 return kInfinity;
00493 }
00494 }
00495
00496 G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y(),
00497 vRho2 = v.perp2(), intersection,
00498 B = (k1 * p.z() + k2 - rho2) * vRho2;
00499
00500 if ( ( (rho2 > paraRho2) && (sqr(rho2-paraRho2-0.25*tol2) > tol2*paraRho2) )
00501 || (p.z() < - dz+kCarTolerance)
00502 || (p.z() > dz-kCarTolerance) )
00503 {
00504
00505
00506 if(vRho2<tol2)
00507 {
00508 intersection = ((rho2 - k2)/k1 - p.z())/v.z();
00509 if(intersection < 0) { return kInfinity; }
00510 else if(std::fabs(p.z() + v.z() * intersection) <= dz)
00511 {
00512 return intersection;
00513 }
00514 else
00515 {
00516 return kInfinity;
00517 }
00518 }
00519 else if(A*A + B < 0)
00520 {
00521 return kInfinity;
00522 }
00523 else
00524 {
00525 intersection = (A - std::sqrt(B + sqr(A))) / vRho2;
00526 if(intersection < 0)
00527 {
00528 return kInfinity;
00529 }
00530 else if(std::fabs(p.z() + intersection * v.z()) < dz + tolh)
00531 {
00532 return intersection;
00533 }
00534 else
00535 {
00536 return kInfinity;
00537 }
00538 }
00539 }
00540 else if(sqr(rho2 - paraRho2 - .25 * tol2) <= tol2 * paraRho2)
00541 {
00542
00543
00544 G4ThreeVector normal(p.x(), p.y(), -k1/2);
00545 if(normal.dot(v) <= 0)
00546 { return 0; }
00547 else
00548 { return kInfinity; }
00549 }
00550 else
00551 {
00552 std::ostringstream message;
00553 if(Inside(p) == kInside)
00554 {
00555 message << "Point p is inside! - " << GetName() << G4endl;
00556 }
00557 else
00558 {
00559 message << "Likely a problem in this function, for solid: " << GetName()
00560 << G4endl;
00561 }
00562 message << " p = " << p * (1/mm) << " mm" << G4endl
00563 << " v = " << v * (1/mm) << " mm";
00564 G4Exception("G4Paraboloid::DistanceToIn(p,v)", "GeomSolids1002",
00565 JustWarning, message);
00566 return 0;
00567 }
00568 }
00569
00571
00572
00573
00574
00575 G4double G4Paraboloid::DistanceToIn(const G4ThreeVector& p) const
00576 {
00577 G4double safz = -dz+std::fabs(p.z());
00578 if(safz<0) { safz=0; }
00579 G4double safr = kInfinity;
00580
00581 G4double rho = p.x()*p.x()+p.y()*p.y();
00582 G4double paraRho = (p.z()-k2)/k1;
00583 G4double sqrho = std::sqrt(rho);
00584
00585 if(paraRho<0)
00586 {
00587 safr=sqrho-r2;
00588 if(safr>safz) { safz=safr; }
00589 return safz;
00590 }
00591
00592 G4double sqprho = std::sqrt(paraRho);
00593 G4double dRho = sqrho-sqprho;
00594 if(dRho<0) { return safz; }
00595
00596 G4double talf = -2.*k1*sqprho;
00597 G4double tmp = 1+talf*talf;
00598 if(tmp<0.) { return safz; }
00599
00600 G4double salf = talf/std::sqrt(tmp);
00601 safr = std::fabs(dRho*salf);
00602 if(safr>safz) { safz=safr; }
00603
00604 return safz;
00605 }
00606
00608
00609
00610
00611 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p,
00612 const G4ThreeVector& v,
00613 const G4bool calcNorm,
00614 G4bool *validNorm,
00615 G4ThreeVector *n ) const
00616 {
00617 G4double rho2 = p.perp2(), paraRho2 = std::fabs(k1 * p.z() + k2);
00618 G4double vRho2 = v.perp2(), intersection;
00619 G4double tol2 = kCarTolerance*kCarTolerance;
00620 G4double tolh = 0.5*kCarTolerance;
00621
00622 if(calcNorm) { *validNorm = false; }
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632 G4double A = k1 / 2 * v.z() - p.x() * v.x() - p.y() * v.y();
00633
00634
00635
00636 G4double B = (-rho2 + paraRho2) * vRho2;
00637
00638 if ( rho2 < paraRho2 && sqr(rho2 - paraRho2 - 0.25 * tol2) > tol2 * paraRho2
00639 && std::fabs(p.z()) < dz - kCarTolerance)
00640 {
00641
00642
00643 if(v.z() > 0)
00644 {
00645
00646
00647
00648
00649
00650 intersection = (dz - p.z()) / v.z();
00651 G4ThreeVector ip = p + intersection * v;
00652
00653 if(ip.perp2() < sqr(r2 + kCarTolerance))
00654 {
00655 if(calcNorm)
00656 {
00657 *n = G4ThreeVector(0, 0, 1);
00658 if(r2 < tolh || ip.perp2() > sqr(r2 - tolh))
00659 {
00660 *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
00661 *n = n->unit();
00662 }
00663 *validNorm = true;
00664 }
00665 return intersection;
00666 }
00667 }
00668 else if(v.z() < 0)
00669 {
00670
00671
00672
00673
00674
00675 intersection = (-dz - p.z()) / v.z();
00676 G4ThreeVector ip = p + intersection * v;
00677
00678 if(ip.perp2() < sqr(r1 + tolh))
00679 {
00680 if(calcNorm)
00681 {
00682 *n = G4ThreeVector(0, 0, -1);
00683 if(r1 < tolh || ip.perp2() > sqr(r1 - tolh))
00684 {
00685 *n += G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
00686 *n = n->unit();
00687 }
00688 *validNorm = true;
00689 }
00690 return intersection;
00691 }
00692 }
00693
00694
00695
00696 if(vRho2 == 0)
00697 {
00698 intersection = ((rho2 - k2)/k1 - p.z())/v.z();
00699 if(calcNorm)
00700 {
00701 G4ThreeVector intersectionP = p + v * intersection;
00702 *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
00703 *n = n->unit();
00704
00705 *validNorm = true;
00706 }
00707 return intersection;
00708 }
00709 else if( ((A <= 0) && (B >= sqr(A) * (sqr(vRho2) - 1))) || (A >= 0))
00710 {
00711
00712
00713
00714
00715 A = A/vRho2;
00716 B = (k1 * p.z() + k2 - rho2)/vRho2;
00717 intersection = B/(-A + std::sqrt(B + sqr(A)));
00718 if(calcNorm)
00719 {
00720 G4ThreeVector intersectionP = p + v * intersection;
00721 *n = G4ThreeVector(intersectionP.x(), intersectionP.y(), -k1/2);
00722 *n = n->unit();
00723 *validNorm = true;
00724 }
00725 return intersection;
00726 }
00727 std::ostringstream message;
00728 message << "There is no intersection between given line and solid!"
00729 << G4endl
00730 << " p = " << p << G4endl
00731 << " v = " << v;
00732 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
00733 JustWarning, message);
00734
00735 return kInfinity;
00736 }
00737 else if ( (rho2 < paraRho2 + kCarTolerance
00738 || sqr(rho2 - paraRho2 - 0.25 * tol2) < tol2 * paraRho2 )
00739 && std::fabs(p.z()) < dz + tolh)
00740 {
00741
00742
00743 G4ThreeVector normal = G4ThreeVector (p.x(), p.y(), -k1/2);
00744
00745 if(std::fabs(p.z()) > dz - tolh)
00746 {
00747
00748
00749 if( ((v.z() > 0) && (p.z() > 0)) || ((v.z() < 0) && (p.z() < 0)) )
00750 {
00751 if(calcNorm)
00752 {
00753 *validNorm = true;
00754 if(p.z() > 0)
00755 { *n = G4ThreeVector(0, 0, 1); }
00756 else
00757 { *n = G4ThreeVector(0, 0, -1); }
00758 }
00759 return 0;
00760 }
00761
00762 if(v.z() == 0)
00763 {
00764
00765
00766
00767
00768 G4double r = (p.z() > 0)? r2 : r1;
00769 G4double pDotV = p.dot(v);
00770 A = vRho2 * ( sqr(r) - sqr(p.x()) - sqr(p.y()));
00771 intersection = (-pDotV + std::sqrt(A + sqr(pDotV))) / vRho2;
00772
00773 if(calcNorm)
00774 {
00775 *validNorm = true;
00776
00777 *n = (G4ThreeVector(0, 0, p.z()/std::fabs(p.z()))
00778 + G4ThreeVector(p.x() + v.x() * intersection, p.y() + v.y()
00779 * intersection, -k1/2).unit()).unit();
00780 }
00781 return intersection;
00782 }
00783 }
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805 if(v.z() > 0)
00806 {
00807
00808
00809 intersection = (dz - p.z()) / v.z();
00810 G4ThreeVector ip = p + intersection * v;
00811
00812 if(ip.perp2() < sqr(r2 - tolh))
00813 {
00814 if(calcNorm)
00815 {
00816 *validNorm = true;
00817 *n = G4ThreeVector(0, 0, 1);
00818 }
00819 return intersection;
00820 }
00821 else if(ip.perp2() < sqr(r2 + tolh))
00822 {
00823 if(calcNorm)
00824 {
00825 *validNorm = true;
00826 *n = G4ThreeVector(0, 0, 1)
00827 + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
00828 *n = n->unit();
00829 }
00830 return intersection;
00831 }
00832 }
00833 if( v.z() < 0)
00834 {
00835
00836
00837 intersection = (-dz - p.z()) / v.z();
00838 G4ThreeVector ip = p + intersection * v;
00839
00840 if(ip.perp2() < sqr(r1 - tolh))
00841 {
00842 if(calcNorm)
00843 {
00844 *validNorm = true;
00845 *n = G4ThreeVector(0, 0, -1);
00846 }
00847 return intersection;
00848 }
00849 else if(ip.perp2() < sqr(r1 + tolh))
00850 {
00851 if(calcNorm)
00852 {
00853 *validNorm = true;
00854 *n = G4ThreeVector(0, 0, -1)
00855 + G4ThreeVector(ip.x(), ip.y(), - k1 / 2).unit();
00856 *n = n->unit();
00857 }
00858 return intersection;
00859 }
00860 }
00861
00862
00863
00864 if(std::fabs(vRho2) > tol2)
00865 {
00866 A = A/vRho2;
00867 B = (k1 * p.z() + k2 - rho2);
00868 if(std::fabs(B)>kCarTolerance)
00869 {
00870 B = (B)/vRho2;
00871 intersection = B/(-A + std::sqrt(B + sqr(A)));
00872 }
00873 else
00874 {
00875 if(normal.dot(v) >= 0)
00876 {
00877 if(calcNorm)
00878 {
00879 *validNorm = true;
00880 *n = normal.unit();
00881 }
00882 return 0;
00883 }
00884 intersection = 2.*A;
00885 }
00886 }
00887 else
00888 {
00889 intersection = ((rho2 - k2) / k1 - p.z()) / v.z();
00890 }
00891
00892 if(calcNorm)
00893 {
00894 *validNorm = true;
00895 *n = G4ThreeVector(p.x() + intersection * v.x(), p.y()
00896 + intersection * v.y(), - k1 / 2);
00897 *n = n->unit();
00898 }
00899 return intersection;
00900 }
00901 else
00902 {
00903 #ifdef G4SPECSDEBUG
00904 if(kOutside == Inside(p))
00905 {
00906 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
00907 JustWarning, "Point p is outside!");
00908 }
00909 else
00910 G4Exception("G4Paraboloid::DistanceToOut(p,v,...)", "GeomSolids1002",
00911 JustWarning, "There's an error in this functions code.");
00912 #endif
00913 return kInfinity;
00914 }
00915 return 0;
00916 }
00917
00919
00920
00921
00922 G4double G4Paraboloid::DistanceToOut(const G4ThreeVector& p) const
00923 {
00924 G4double safe=0.0,rho,safeR,safeZ ;
00925 G4double tanRMax,secRMax,pRMax ;
00926
00927 #ifdef G4SPECSDEBUG
00928 if( Inside(p) == kOutside )
00929 {
00930 G4cout << G4endl ;
00931 DumpInfo();
00932 std::ostringstream message;
00933 G4int oldprc = message.precision(16);
00934 message << "Point p is outside !?" << G4endl
00935 << "Position:" << G4endl
00936 << " p.x() = " << p.x()/mm << " mm" << G4endl
00937 << " p.y() = " << p.y()/mm << " mm" << G4endl
00938 << " p.z() = " << p.z()/mm << " mm";
00939 message.precision(oldprc) ;
00940 G4Exception("G4Paraboloid::DistanceToOut(p)", "GeomSolids1002",
00941 JustWarning, message);
00942 }
00943 #endif
00944
00945 rho = p.perp();
00946 safeZ = dz - std::fabs(p.z()) ;
00947
00948 tanRMax = (r2 - r1)*0.5/dz ;
00949 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
00950 pRMax = tanRMax*p.z() + (r1+r2)*0.5 ;
00951 safeR = (pRMax - rho)/secRMax ;
00952
00953 if (safeZ < safeR) { safe = safeZ; }
00954 else { safe = safeR; }
00955 if ( safe < 0.5 * kCarTolerance ) { safe = 0; }
00956 return safe ;
00957 }
00958
00960
00961
00962
00963 G4GeometryType G4Paraboloid::GetEntityType() const
00964 {
00965 return G4String("G4Paraboloid");
00966 }
00967
00969
00970
00971
00972 G4VSolid* G4Paraboloid::Clone() const
00973 {
00974 return new G4Paraboloid(*this);
00975 }
00976
00978
00979
00980
00981 std::ostream& G4Paraboloid::StreamInfo( std::ostream& os ) const
00982 {
00983 G4int oldprc = os.precision(16);
00984 os << "-----------------------------------------------------------\n"
00985 << " *** Dump for solid - " << GetName() << " ***\n"
00986 << " ===================================================\n"
00987 << " Solid type: G4Paraboloid\n"
00988 << " Parameters: \n"
00989 << " z half-axis: " << dz/mm << " mm \n"
00990 << " radius at -dz: " << r1/mm << " mm \n"
00991 << " radius at dz: " << r2/mm << " mm \n"
00992 << "-----------------------------------------------------------\n";
00993 os.precision(oldprc);
00994
00995 return os;
00996 }
00997
00999
01000
01001
01002 G4ThreeVector G4Paraboloid::GetPointOnSurface() const
01003 {
01004 G4double A = (fSurfaceArea == 0)? CalculateSurfaceArea(): fSurfaceArea;
01005 G4double z = RandFlat::shoot(0.,1.);
01006 G4double phi = RandFlat::shoot(0., twopi);
01007 if(pi*(sqr(r1) + sqr(r2))/A >= z)
01008 {
01009 G4double rho;
01010 if(pi * sqr(r1) / A > z)
01011 {
01012 rho = r1 * std::sqrt(RandFlat::shoot(0., 1.));
01013 return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), -dz);
01014 }
01015 else
01016 {
01017 rho = r2 * std::sqrt(RandFlat::shoot(0., 1));
01018 return G4ThreeVector(rho * std::cos(phi), rho * std::sin(phi), dz);
01019 }
01020 }
01021 else
01022 {
01023 z = RandFlat::shoot(0., 1.)*2*dz - dz;
01024 return G4ThreeVector(std::sqrt(z*k1 + k2)*std::cos(phi),
01025 std::sqrt(z*k1 + k2)*std::sin(phi), z);
01026 }
01027 }
01028
01029 G4ThreeVectorList*
01030 G4Paraboloid::CreateRotatedVertices(const G4AffineTransform& pTransform,
01031 G4int& noPolygonVertices) const
01032 {
01033 G4ThreeVectorList *vertices;
01034 G4ThreeVector vertex;
01035 G4double meshAnglePhi, cosMeshAnglePhiPer2,
01036 crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi,
01037 sRho, dRho, rho, lastRho = 0., swapRho;
01038 G4double rx, ry, rz, k3, k4, zm;
01039 G4int crossSectionPhi, noPhiCrossSections, noRhoSections;
01040
01041
01042
01043 noPhiCrossSections = G4int(twopi/kMeshAngleDefault)+1;
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054 meshAnglePhi=twopi/(noPhiCrossSections-1);
01055
01056 sAnglePhi = -meshAnglePhi*0.5*0;
01057 cosMeshAnglePhiPer2 = std::cos(meshAnglePhi / 2.);
01058
01059 noRhoSections = G4int(pi/2/kMeshAngleDefault) + 1;
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072 sRho = r1;
01073 dRho = (r2 - r1) / double(noRhoSections - 1);
01074
01075 vertices=new G4ThreeVectorList();
01076
01077 if (vertices)
01078 {
01079 for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
01080 crossSectionPhi++)
01081 {
01082 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
01083 coscrossAnglePhi=std::cos(crossAnglePhi);
01084 sincrossAnglePhi=std::sin(crossAnglePhi);
01085 lastRho = 0;
01086 for (int iRho=0; iRho < noRhoSections;
01087 iRho++)
01088 {
01089
01090
01091 if(iRho == noRhoSections - 1)
01092 {
01093 rho = r2;
01094 }
01095 else
01096 {
01097 rho = iRho * dRho + sRho;
01098
01099
01100
01101
01102 k3 = k1 / (2*rho + dRho);
01103 k4 = rho - k3 * (sqr(rho) - k2) / k1;
01104 zm = (sqr(k1 / (2 * k3)) - k2) / k1;
01105 rho += std::sqrt(k1 * zm + k2) - zm * k3 - k4;
01106 }
01107
01108 rho += (1 / cosMeshAnglePhiPer2 - 1) * (iRho * dRho + sRho);
01109
01110 if(rho < lastRho)
01111 {
01112 swapRho = lastRho;
01113 lastRho = rho + dRho;
01114 rho = swapRho;
01115 }
01116 else
01117 {
01118 lastRho = rho + dRho;
01119 }
01120
01121 rx = coscrossAnglePhi*rho;
01122 ry = sincrossAnglePhi*rho;
01123 rz = (sqr(iRho * dRho + sRho) - k2) / k1;
01124 vertex = G4ThreeVector(rx,ry,rz);
01125 vertices->push_back(pTransform.TransformPoint(vertex));
01126 }
01127 }
01128 noPolygonVertices = noRhoSections ;
01129 }
01130 else
01131 {
01132 DumpInfo();
01133 G4Exception("G4Paraboloid::CreateRotatedVertices()",
01134 "GeomSolids0003", FatalException,
01135 "Error in allocation of vertices. Out of memory !");
01136 }
01137 return vertices;
01138 }
01139
01141
01142
01143
01144 void G4Paraboloid::DescribeYourselfTo (G4VGraphicsScene& scene) const
01145 {
01146 scene.AddSolid(*this);
01147 }
01148
01149 G4NURBS* G4Paraboloid::CreateNURBS () const
01150 {
01151
01152
01153 return new G4NURBSbox(r1, r1, dz);
01154 }
01155
01156 G4Polyhedron* G4Paraboloid::CreatePolyhedron () const
01157 {
01158 return new G4PolyhedronParaboloid(r1, r2, dz, 0., twopi);
01159 }
01160
01161
01162 G4Polyhedron* G4Paraboloid::GetPolyhedron () const
01163 {
01164 if (!fpPolyhedron ||
01165 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01166 fpPolyhedron->GetNumberOfRotationSteps())
01167 {
01168 delete fpPolyhedron;
01169 fpPolyhedron = CreatePolyhedron();
01170 }
01171 return fpPolyhedron;
01172 }