50 axis0min, axis1min, axis0max, axis1max),
55 G4Exception(
"G4TwistTubsSide::G4TwistTubsSide()",
"GeomSolids0002",
90 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
221 for (
auto i=0; i<2; ++i)
223 distance[i] = kInfinity;
226 gxx[i].
set(kInfinity, kInfinity, kInfinity);
247 isvalid[0], 0, validate, &gp, &gv);
256 G4double b = fKappa * (v.
x() * p.
z() + v.
z() * p.
x()) - v.
y();
267 distance[0] = - c / b;
268 xx[0] = p + distance[0]*v;
273 areacode[0] = GetAreaCode(xx[0]);
276 if (distance[0] >= 0) isvalid[0] =
true;
281 areacode[0] = GetAreaCode(xx[0],
false);
284 if (distance[0] >= 0) isvalid[0] =
true;
293 if (distance[0] >= 0) isvalid[0] =
true;
297 distance[0] = kInfinity;
299 areacode[0], isvalid[0],
300 0, validate, &gp, &gv);
306 isvalid[0], 1, validate, &gp, &gv);
319 isvalid[0], 0, validate, &gp, &gv);
328 G4double tmpdist[2] = {kInfinity, kInfinity};
331 G4bool tmpisvalid[2] = {
false,
false};
333 for (
auto i=0; i<2; ++i)
340 if ( b *
D < 0 && std::fabs(bminusD /
D) < protection )
343 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
347 tmpdist[i] = factor * bminusD;
351 tmpxx[i] = p + tmpdist[i]*v;
355 tmpareacode[i] = GetAreaCode(tmpxx[i]);
358 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
364 tmpareacode[i] = GetAreaCode(tmpxx[i],
false);
367 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
374 if (tmpxx[i].x() > 0)
377 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
380 tmpdist[i] = kInfinity;
386 if (tmpdist[0] <= tmpdist[1])
388 distance[0] = tmpdist[0];
389 distance[1] = tmpdist[1];
394 areacode[0] = tmpareacode[0];
395 areacode[1] = tmpareacode[1];
396 isvalid[0] = tmpisvalid[0];
397 isvalid[1] = tmpisvalid[1];
401 distance[0] = tmpdist[1];
402 distance[1] = tmpdist[0];
407 areacode[0] = tmpareacode[1];
408 areacode[1] = tmpareacode[0];
409 isvalid[0] = tmpisvalid[1];
410 isvalid[1] = tmpisvalid[0];
414 isvalid[0], 2, validate, &gp, &gv);
416 isvalid[1], 2, validate, &gp, &gv);
420 for (
G4int k=0; k<2; ++k)
422 if (!isvalid[k])
continue;
424 G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
425 * xx[k].z() , xx[k].z());
426 G4double deltaY = (xx[k] - xxonsurface).mag();
433 for (l=0; l<maxcount; ++l)
437 surfacenormal, xx[k]);
438 deltaY = (xx[k] - xxonsurface).mag();
439 if (deltaY > lastdeltaY) { }
443 xxonsurface.
set(xx[k].x(),
444 fKappa * std::fabs(xx[k].x()) * xx[k].z(),
449 std::ostringstream message;
450 message <<
"Exceeded maxloop count!" <<
G4endl
451 <<
" maxloop count " << maxcount;
452 G4Exception(
"G4TwistTubsFlatSide::DistanceToSurface(p,v)",
465 isvalid[0], 0, validate, &gp, &gv);
492 for (
auto i=0; i<2; ++i)
494 distance[i] = kInfinity;
496 gxx[i].
set(kInfinity, kInfinity, kInfinity);
504 G4int parity = (fKappa >= 0 ? 1 : -1);
514 for (
auto i=0; i<2; ++i)
519 if ((gp - lastgxx[0]).mag() < halftol
520 || (gp - lastgxx[1]).mag() < halftol)
580 else if (p.
z() <
C.z())
591 else if (p.
z() <
A.z())
601 for (
auto i=0; i<2; ++i)
606 B = x0[i] + ((
A.z() - x0[i].z()) / d[i].z()) * d[i];
612 D = x0[i] + ((
C.z() - x0[i].z()) / d[i].z()) * d[i];
621 G4double test = (
A.z() -
C.z()) * parity * pside;
677 xxacb, nacb) * parity;
679 xxcad, ncad) * parity;
682 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol)
684 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
694 if (distToACB * distToCAD > 0 && distToACB < 0)
699 distance[0] = DistanceToPlane(p,
A,
B,
C,
D, parity, xx, normal);
703 if (distToACB * distToCAD > 0)
707 if (distToACB <= distToCAD)
709 distance[0] = distToACB;
714 distance[0] = distToCAD;
724 distance[0] = distToACB;
729 distance[0] = distToCAD;
764 xxanm, nanm) * parity;
766 xxcmn, ncmn) * parity;
769 if (distToanm * distTocmn > 0 && distToanm < 0)
773 "Point p is behind the surfaces.");
777 if (std::fabs(distToanm) <= halftol)
783 else if (std::fabs(distTocmn) <= halftol)
790 if (distToanm <= distTocmn)
802 return DistanceToPlane(p,
A,
M,
N,
D, parity, xx, n);
817 return DistanceToPlane(p, C,
N,
M, B, parity, xx, n);
848 if (xx.
x() <=
fAxisMin[xaxis] - ctol) isoutside =
true;
851 else if (xx.
x() >
fAxisMax[xaxis] - ctol)
854 if (xx.
x() >=
fAxisMax[xaxis] + ctol) isoutside =
true;
865 if (xx.
z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
868 else if (xx.
z() >
fAxisMax[zaxis] - ctol)
874 if (xx.
z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
882 G4int tmpareacode = areacode & (~sInside);
883 areacode = tmpareacode;
930 "Feature NOT implemented !");
938void G4TwistTubsSide::SetCorners(
G4double endInnerRad[2],
953 x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
954 y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
959 x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
960 y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
965 x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
966 y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
971 x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
972 y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
979 std::ostringstream message;
980 message <<
"Feature NOT implemented !" <<
G4endl
982 <<
" fAxis[1] = " <<
fAxis[1];
991void G4TwistTubsSide::SetCorners()
995 "Method NOT implemented !");
1001void G4TwistTubsSide::SetBoundaries()
1011 direction = direction.
unit();
1017 direction = direction.
unit();
1023 direction = direction.
unit();
1029 direction = direction.
unit();
1036 std::ostringstream message;
1037 message <<
"Feature NOT implemented !" <<
G4endl
1039 <<
" fAxis[1] = " <<
fAxis[1];
1061 for (
G4int i = 0 ; i<n ; ++i )
1065 for (
G4int j = 0 ; j<k ; ++j )
1067 nnode =
GetNode(i,j,k,n,iside) ;
1074 x = xmin + j*(xmax-xmin)/(k-1) ;
1078 x = xmax - j*(xmax-xmin)/(k-1) ;
1083 xyz[nnode][0] = p.
x() ;
1084 xyz[nnode][1] = p.
y() ;
1085 xyz[nnode][2] = p.
z() ;
1087 if ( i<n-1 && j<k-1 )
1089 nface =
GetFace(i,j,k,n,iside) ;
1092 * (
GetNode(i ,j ,k,n,iside)+1) ;
1094 * (
GetNode(i+1,j ,k,n,iside)+1) ;
1096 * (
GetNode(i+1,j+1,k,n,iside)+1) ;
1098 * (
GetNode(i ,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)
G4double GetBoundaryMax(G4double phi) override
G4double GetBoundaryMin(G4double phi) override
G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false) override
~G4TwistTubsSide() override
void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside) override
G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false) override
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 DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol) override
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