82 if ( ! (fDx1 == fDx2 && fDx3 == fDx4 ) ) {
83 std::ostringstream message;
84 message <<
"TwistedTrapBoxSide is not used as a the side of a box: "
87 G4Exception(
"G4TwistBoxSide::G4TwistBoxSide()",
"GeomSolids0002",
97 fTAlph = std::tan(fAlph) ;
103 fDx4plus2 = fDx4 + fDx2 ;
104 fDx4minus2 = fDx4 - fDx2 ;
105 fDx3plus1 = fDx3 + fDx1 ;
106 fDx3minus1 = fDx3 - fDx1 ;
107 fDy2plus1 = fDy2 + fDy1 ;
108 fDy2minus1 = fDy2 - fDy1 ;
110 fa1md1 = 2*fDx2 - 2*fDx1 ;
111 fa2md2 = 2*fDx4 - 2*fDx3 ;
114 fPhiTwist = PhiTwist ;
115 fAngleSide = AngleSide ;
117 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
118 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
135 :
G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
136 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
137 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.),
138 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.),
178 GetPhiUAtX(xx,phi,u) ;
212 G4bool IsParallel = false ;
213 G4bool IsConverged = false ;
255 G4bool tmpisvalid = false ;
257 std::vector<Intersection> xbuf ;
264 G4double phixz = fPhiTwist * ( p.
x() * v.
z() - p.
z() * v.
x() ) ;
265 G4double phiyz = fPhiTwist * ( p.
y() * v.
z() - p.
z() * v.
y() ) ;
271 if ( (std::fabs(p.
z()) <= L) && fPhiTwist ) {
273 phi = p.
z() * fPhiTwist /
L ;
275 q = (2.* fPhiTwist*((v.
x() - fTAlph*v.
y())*std::cos(phi)
276 + (fTAlph*v.
x() + v.
y())*std::sin(phi)));
279 u = (2*(-(fdeltaY*phi*v.
x()) + fPhiTwist*p.
y()*v.
x()
280 + fdeltaX*phi*v.
y() - fPhiTwist*p.
x()*v.
y())
281 + (fDx4plus2*fPhiTwist + 2*fDx4minus2*phi)
282 * (v.
y()*std::cos(phi) - v.
x()*std::sin(phi))) / q;
291 xbuf.push_back(xbuftmp) ;
302 areacode[0], isvalid[0],
303 0, validate, &gp, &gv);
319 c[7] = -14400*(-2*phixz + 2*fTAlph*phiyz + fDx4plus2*fPhiTwist*v.
z()) ;
320 c[6] = 28800*(phiyz + 2*fDz*v.
x() - (fdeltaX + fDx4minus2)*v.
z() + fTAlph*(phixz - 2*fDz*v.
y() + fdeltaY*v.
z())) ;
321 c[5] = -1200*(10*phixz - 48*fDz*v.
y() + 24*fdeltaY*v.
z() + fDx4plus2*fPhiTwist*v.
z() - 2*fTAlph*(5*phiyz + 24*fDz*v.
x() - 12*fdeltaX*v.
z())) ;
322 c[4] = -2400*(phiyz + 10*fDz*v.
x() - 5*fdeltaX*v.
z() + fDx4minus2*v.
z() + fTAlph*(phixz - 10*fDz*v.
y() + 5*fdeltaY*v.
z())) ;
323 c[3] = 24*(2*phixz - 200*fDz*v.
y() + 100*fdeltaY*v.
z() - fDx4plus2*fPhiTwist*v.
z() - 2*fTAlph*(phiyz + 100*fDz*v.
x() - 50*fdeltaX*v.
z())) ;
324 c[2] = -16*(7*fTAlph* phixz + 7*phiyz - 6*fDz*v.
x() + 6*fDz*fTAlph*v.
y() + 3*(fdeltaX + fDx4minus2 - fdeltaY*fTAlph)*v.
z()) ;
325 c[1] = 4*(9*phixz - 9*fTAlph*phiyz - 56*fDz*fTAlph*v.
x() - 56*fDz*v.
y() + 28*(fdeltaY + fdeltaX*fTAlph)*v.
z()) ;
326 c[0] = 36*(2* fDz*(v.
x() - fTAlph*v.
y()) - fdeltaX*v.
z() + fdeltaY*fTAlph*v.
z()) ;
330 G4cout <<
"coef = " << c[0] <<
" "
344 for (
G4int i = 0 ; i<num ; i++ ) {
345 if ( (si[i]==0.0) && fPhiTwist ) {
347 G4cout <<
"Solution " << i <<
" : " << srd[i] <<
G4endl ;
349 phi = std::fmod(srd[i] , pihalf) ;
351 u = (2*phiyz + 4*fDz*phi*v.
y() - 2*fdeltaY*phi*v.
z()
352 - fDx4plus2*fPhiTwist*v.
z()*std::sin(phi)
353 - 2*fDx4minus2*phi*v.
z()*std::sin(phi))
354 /(2*fPhiTwist*v.
z()*std::cos(phi)
355 + 2*fPhiTwist*fTAlph*v.
z()*std::sin(phi)) ;
363 xbuf.push_back(xbuftmp) ;
366 G4cout <<
"solution " << i <<
" = " << phi <<
" , " << u <<
G4endl ;
385 for (
size_t k = 0 ; k<xbuf.size() ; k++ ) {
388 G4cout <<
"Solution " << k <<
" : "
389 <<
"reconstructed phiR = " << xbuf[k].phi
390 <<
", uR = " << xbuf[k].u <<
G4endl ;
396 IsConverged = false ;
398 for (
G4int i = 1 ; i<maxint ; i++ ) {
400 xxonsurface = SurfacePoint(phi,u) ;
401 surfacenormal = NormAng(phi,u) ;
403 deltaX = ( tmpxx - xxonsurface ).mag() ;
404 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
405 if ( theta < 0.001 ) {
414 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist <<
", " << deltaX <<
G4endl ;
418 GetPhiUAtX(tmpxx, phi, u) ;
421 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
424 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
430 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
433 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
443 tmpareacode = GetAreaCode(tmpxx);
445 if (tmpdist >= 0) tmpisvalid =
true;
448 tmpareacode = GetAreaCode(tmpxx,
false);
450 if (tmpdist >= 0) tmpisvalid =
true;
455 "Feature NOT implemented !");
467 xbuf[k].distance = tmpdist ;
468 xbuf[k].areacode = tmpareacode ;
469 xbuf[k].isvalid = tmpisvalid ;
478 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
484 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() ,
EqualIntersection ) , xbuf.end() ) ;
489 G4int nxxtmp = xbuf.size() ;
491 if ( nxxtmp<2 || IsParallel ) {
509 xbuf.push_back(xbuftmp) ;
525 xbuf.push_back(xbuftmp) ;
527 for (
size_t k = nxxtmp ; k<xbuf.size() ; k++ ) {
530 G4cout <<
"Solution " << k <<
" : "
531 <<
"reconstructed phiR = " << xbuf[k].phi
532 <<
", uR = " << xbuf[k].u <<
G4endl ;
538 IsConverged = false ;
540 for (
G4int i = 1 ; i<maxint ; i++ ) {
542 xxonsurface = SurfacePoint(phi,u) ;
543 surfacenormal = NormAng(phi,u) ;
545 deltaX = ( tmpxx - xxonsurface ).mag() ;
546 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
547 if ( theta < 0.001 ) {
555 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist <<
", " << deltaX <<
G4endl ;
559 GetPhiUAtX(tmpxx, phi, u) ;
562 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
565 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
571 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
574 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
584 tmpareacode = GetAreaCode(tmpxx);
586 if (tmpdist >= 0) tmpisvalid =
true;
589 tmpareacode = GetAreaCode(tmpxx,
false);
591 if (tmpdist >= 0) tmpisvalid =
true;
594 G4Exception(
"G4TwistedBoxSide::DistanceToSurface()",
596 "Feature NOT implemented !");
608 xbuf[k].distance = tmpdist ;
609 xbuf[k].areacode = tmpareacode ;
610 xbuf[k].isvalid = tmpisvalid ;
623 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() ,
EqualIntersection ) , xbuf.end() ) ;
626 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
632 for (
size_t i = 0 ; i<xbuf.size() ; i++ ) {
634 distance[i] = xbuf[i].distance;
636 areacode[i] = xbuf[i].areacode ;
637 isvalid[i] = xbuf[i].isvalid ;
640 isvalid[i], nxx, validate, &gp, &gv);
643 G4cout <<
"element Nr. " << i
644 <<
", local Intersection = " << xbuf[i].xx
645 <<
", distance = " << xbuf[i].distance
646 <<
", u = " << xbuf[i].u
647 <<
", phi = " << xbuf[i].phi
648 <<
", isvalid = " << xbuf[i].isvalid
657 G4cout << nxx <<
" possible physical solutions found" <<
G4endl ;
658 for (
G4int k= 0 ; k< nxx ; k++ ) {
659 G4cout <<
"global intersection Point found: " << gxx[k] <<
G4endl ;
715 for (
G4int i = 1 ; i<maxint ; i++ ) {
717 xxonsurface = SurfacePoint(phiR,uR) ;
718 surfacenormal = NormAng(phiR,uR) ;
720 deltaX = ( xx - xxonsurface ).mag() ;
723 G4cout <<
"i = " << i <<
", distance = " << distance[0] <<
", " << deltaX <<
G4endl ;
728 GetPhiUAtX(xx, phiR, uR) ;
730 if ( deltaX <= ctol ) { break ; }
737 G4double uMax = GetBoundaryMax(phiR) ;
739 if ( phiR > halfphi ) phiR = halfphi ;
740 if ( phiR < -halfphi ) phiR = -halfphi ;
741 if ( uR > uMax ) uR = uMax ;
742 if ( uR < -uMax ) uR = -uMax ;
744 xxonsurface = SurfacePoint(phiR,uR) ;
745 distance[0] = ( p - xx ).mag() ;
746 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
751 G4cout <<
"refined solution " << phiR <<
" , " << uR <<
" , " <<
G4endl ;
760 G4cout <<
"intersection Point found: " << gxx[0] <<
G4endl ;
785 GetPhiUAtX(xx, phi,yprime ) ;
787 G4double fYAxisMax = GetBoundaryMax(phi) ;
792 G4cout <<
"GetAreaCode: yprime = " << yprime <<
G4endl ;
793 G4cout <<
"Intervall is " << fYAxisMin <<
" to " << fYAxisMax <<
G4endl ;
808 if (yprime < fYAxisMin + ctol) {
810 if (yprime <= fYAxisMin - ctol) isoutside =
true;
812 }
else if (yprime > fYAxisMax - ctol) {
814 if (yprime >= fYAxisMax + ctol) isoutside =
true;
824 if (xx.
z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
826 }
else if (xx.
z() >
fAxisMax[zaxis] - ctol) {
829 if (areacode & sBoundary) areacode |=
sCorner;
831 if (xx.
z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
839 areacode = tmpareacode;
840 }
else if ((areacode & sBoundary) !=
sBoundary) {
848 if (yprime < fYAxisMin ) {
850 }
else if (yprime > fYAxisMax) {
858 if (areacode & sBoundary) areacode |=
sCorner;
863 if (areacode & sBoundary) areacode |=
sCorner;
867 if ((areacode & sBoundary) !=
sBoundary) {
875 "Feature NOT implemented !");
883 void G4TwistBoxSide::SetCorners()
894 x = -fdeltaX/2. + (fDx2 - fDy1*fTAlph)*std::cos(fPhiTwist/2.) - fDy1*std::sin(fPhiTwist/2.) ;
895 y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/2.) + (-fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
902 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
903 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
910 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
911 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
918 x = fdeltaX/2. + (fDx4 - fDy2*fTAlph)*std::cos(fPhiTwist/2.) + fDy2*std::sin(fPhiTwist/2.) ;
919 y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2.) + (fDx4 - fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
928 "Method NOT implemented !");
935 void G4TwistBoxSide::SetBoundaries()
946 direction = direction.
unit();
952 direction = direction.
unit();
958 direction = direction.
unit();
964 direction = direction.
unit();
972 "Feature NOT implemented !");
984 phi = p.
z()/(2*fDz)*fPhiTwist ;
986 u = -(fTAlph*(fDx4plus2*fPhiTwist + 2*fDx4minus2*phi) + 2*(fdeltaY*phi + fdeltaX*fTAlph*phi - fPhiTwist*(fTAlph*p.
x() + p.
y()))* std::cos(phi) + 2*(-(fdeltaX*phi) + fdeltaY*fTAlph*phi + fPhiTwist*(p.
x() - fTAlph*p.
y()))*std::sin(phi))/(2.*(fPhiTwist + fPhiTwist*fTAlph*fTAlph)) ;
1006 GetPhiUAtX( tmpp, phi, u ) ;
1033 for ( i = 0 ; i<
n ; i++ ) {
1035 z = -fDz+i*(2.*fDz)/(n-1) ;
1036 phi = z*fPhiTwist/(2*fDz) ;
1037 b = GetValueB(phi) ;
1039 for ( j = 0 ; j<k ; j++ ) {
1041 nnode =
GetNode(i,j,k,n,iside) ;
1042 u = -b/2 +j*b/(k-1) ;
1043 p = SurfacePoint(phi,u,
true) ;
1045 xyz[nnode][0] = p.
x() ;
1046 xyz[nnode][1] = p.
y() ;
1047 xyz[nnode][2] = p.
z() ;
1049 if ( i<n-1 && j<k-1 ) {
1051 nface =
GetFace(i,j,k,n,iside) ;
static const G4int sAxisZ
void set(double x, double y, double z)
G4int GetAreacode(G4int i) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
static const G4double kInfinity
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4SurfCurNormal fCurrentNormal
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sOutside
G4TwistBoxSide(const G4String &name, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph, G4double AngleSide)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
HepRotation inverse() const
static double normal(HepRandomEngine *eptr)
static const G4int sC0Min1Min
G4GLOB_DLL std::ostream G4cout
G4bool IsOutside(G4int areacode) const
G4double GetDistance(G4int i) const
static const G4int sAxis1
static const G4int sC0Max1Max
static const G4int sBoundary
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
G4bool IsValid(G4int i) const
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
static const G4int sAxis0
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
static const G4int sAxisMin
static const G4int sAxisY
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
static const G4int sInside
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static const G4int sCorner
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sAxisMax
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
G4ThreeVector GetXX(G4int i) const
virtual G4String GetName() const
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
HepRotation & rotateZ(double delta)
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
virtual ~G4TwistBoxSide()
static constexpr double pi
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
static constexpr double L
G4bool DistanceSort(const Intersection &a, const Intersection &b)
static const G4int sC0Max1Min
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
CurrentStatus fCurStatWithV
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)