50 axis0min, axis1min, axis0max, axis1max),
55 G4Exception(
"G4TwistTubsSide::G4TwistTubsSide()",
"GeomSolids0002",
90 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
223 for (
auto i=0; i<2; ++i)
249 isvalid[0], 0, validate, &gp, &gv);
269 distance[0] = - c / b;
270 xx[0] = p + distance[0]*v;
278 if (distance[0] >= 0) isvalid[0] =
true;
286 if (distance[0] >= 0) isvalid[0] =
true;
295 if (distance[0] >= 0) isvalid[0] =
true;
301 areacode[0], isvalid[0],
302 0, validate, &gp, &gv);
308 isvalid[0], 1, validate, &gp, &gv);
321 isvalid[0], 0, validate, &gp, &gv);
333 G4bool tmpisvalid[2] = {
false,
false};
335 for (
auto i=0; i<2; ++i)
342 if ( b *
D < 0 && std::fabs(bminusD /
D) < protection )
345 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
349 tmpdist[i] = factor * bminusD;
353 tmpxx[i] = p + tmpdist[i]*v;
360 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
369 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
376 if (tmpxx[i].x() > 0)
379 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
388 if (tmpdist[0] <= tmpdist[1])
390 distance[0] = tmpdist[0];
391 distance[1] = tmpdist[1];
396 areacode[0] = tmpareacode[0];
397 areacode[1] = tmpareacode[1];
398 isvalid[0] = tmpisvalid[0];
399 isvalid[1] = tmpisvalid[1];
403 distance[0] = tmpdist[1];
404 distance[1] = tmpdist[0];
409 areacode[0] = tmpareacode[1];
410 areacode[1] = tmpareacode[0];
411 isvalid[0] = tmpisvalid[1];
412 isvalid[1] = tmpisvalid[0];
416 isvalid[0], 2, validate, &gp, &gv);
418 isvalid[1], 2, validate, &gp, &gv);
422 for (
G4int k=0; k<2; ++k)
424 if (!isvalid[k])
continue;
427 * xx[k].z() , xx[k].z());
428 G4double deltaY = (xx[k] - xxonsurface).mag();
435 for (l=0; l<maxcount; ++l)
439 surfacenormal, xx[k]);
440 deltaY = (xx[k] - xxonsurface).mag();
441 if (deltaY > lastdeltaY) { }
445 xxonsurface.
set(xx[k].x(),
446 fKappa * std::fabs(xx[k].x()) * xx[k].z(),
451 std::ostringstream message;
452 message <<
"Exceeded maxloop count!" <<
G4endl
453 <<
" maxloop count " << maxcount;
454 G4Exception(
"G4TwistTubsFlatSide::DistanceToSurface(p,v)",
467 isvalid[0], 0, validate, &gp, &gv);
494 for (
auto i=0; i<2; ++i)
516 for (
auto i=0; i<2; ++i)
521 if ((gp - lastgxx[0]).mag() < halftol
522 || (gp - lastgxx[1]).mag() < halftol)
582 else if (p.
z() <
C.z())
593 else if (p.
z() <
A.z())
603 for (
auto i=0; i<2; ++i)
608 B = x0[i] + ((
A.z() - x0[i].
z()) / d[i].z()) * d[i];
614 D = x0[i] + ((
C.z() - x0[i].
z()) / d[i].z()) * d[i];
684 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol)
686 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
696 if (distToACB * distToCAD > 0 && distToACB < 0)
705 if (distToACB * distToCAD > 0)
709 if (distToACB <= distToCAD)
711 distance[0] = distToACB;
716 distance[0] = distToCAD;
726 distance[0] = distToACB;
731 distance[0] = distToCAD;
771 if (distToanm * distTocmn > 0 && distToanm < 0)
775 "Point p is behind the surfaces.");
779 if (std::fabs(distToanm) <= halftol)
785 else if (std::fabs(distTocmn) <= halftol)
792 if (distToanm <= distTocmn)
850 if (xx.
x() <=
fAxisMin[xaxis] - ctol) isoutside =
true;
853 else if (xx.
x() >
fAxisMax[xaxis] - ctol)
856 if (xx.
x() >=
fAxisMax[xaxis] + ctol) isoutside =
true;
867 if (xx.
z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
870 else if (xx.
z() >
fAxisMax[zaxis] - ctol)
876 if (xx.
z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
884 G4int tmpareacode = areacode & (~sInside);
885 areacode = tmpareacode;
932 "Feature NOT implemented !");
955 x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
956 y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
961 x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
962 y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
967 x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
968 y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
973 x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
974 y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
981 std::ostringstream message;
982 message <<
"Feature NOT implemented !" <<
G4endl
984 <<
" fAxis[1] = " <<
fAxis[1];
997 "Method NOT implemented !");
1013 direction = direction.
unit();
1019 direction = direction.
unit();
1025 direction = direction.
unit();
1031 direction = direction.
unit();
1038 std::ostringstream message;
1039 message <<
"Feature NOT implemented !" <<
G4endl
1041 <<
" fAxis[1] = " <<
fAxis[1];
1063 for (
G4int i = 0 ; i<
n ; ++i )
1067 for (
G4int j = 0 ; j<k ; ++j )
1076 x = xmin + j*(xmax-xmin)/(k-1) ;
1080 x = xmax - j*(xmax-xmin)/(k-1) ;
1085 xyz[nnode][0] = p.
x() ;
1086 xyz[nnode][1] = p.
y() ;
1087 xyz[nnode][2] = p.
z() ;
1089 if ( i<
n-1 && j<k-1 )
1098 * (
GetNode(i+1,j+1,k,
n,iside)+1) ;
G4double C(G4double temp)
G4double B(G4double temperature)
G4double D(G4double temp)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Hep3Vector cross(const Hep3Vector &) const
void set(double x, double y, double z)
HepRotation & rotateZ(double delta)
virtual G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &A, const G4ThreeVector &B, const G4ThreeVector &C, const G4ThreeVector &D, const G4int parity, G4ThreeVector &xx, G4ThreeVector &n)
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
virtual G4double GetBoundaryMax(G4double phi)
virtual void SetCorners()
virtual ~G4TwistTubsSide()
virtual G4double GetBoundaryMin(G4double phi)
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
virtual void SetBoundaries()
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
G4TwistTubsSide(const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const G4double kappa, const EAxis axis0=kXAxis, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
G4int GetAreacode(G4int i) const
G4double GetDistance(G4int i) const
G4bool IsValid(G4int i) const
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=nullptr)
G4ThreeVector GetXX(G4int i) const
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=nullptr)
virtual G4int AmIOnLeftSide(const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
static const G4int sC0Min1Min
static const G4int sC0Min1Max
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sOutside
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
static const G4int sAxisMax
static const G4int sAxis0
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
static const G4int sAxisMin
static const G4int sC0Max1Max
static const G4int sAxis1
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sBoundary
static const G4int sAxisZ
G4bool IsOutside(G4int areacode) const
virtual G4double DistanceToBoundary(G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
static const G4int sCorner
static const G4int sC0Max1Min
static const G4int sInside
CurrentStatus fCurStatWithV
static const G4int sAxisX
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4SurfCurNormal fCurrentNormal
virtual void GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
static const G4double kInfinity
static double normal(HepRandomEngine *eptr)
const char * name(G4int ptype)