110 fRot.rotateZ( AngleSide ) ;
125 :
G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
126 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
127 fAngleSide(0.), fDx4plus2(0.), fDx4minus2(0.), fDx3plus1(0.), fDx3minus1(0.),
128 fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), fa2md2(0.), fdeltaX(0.),
211 G4bool IsParallel = false ;
212 G4bool IsConverged = false ;
255 G4bool tmpisvalid = false ;
257 std::vector<Intersection> xbuf ;
272 if ( std::fabs(p.z()) <= L )
279 *(v.y()*std::cos(phi) - v.x()*std::sin(phi))))
282 + 4*
fDy1*(
fTAlph*v.x() + v.y()))*std::sin(phi));
289 xbuf.push_back(xbuftmp) ;
298 areacode[0], isvalid[0],
299 0, validate, &gp, &gv);
324 + 2*
fDy1*(2*phiyz + 20*
fDz*v.x()
329 +
fDy1*(4*phixz - 400*
fDz*v.y()
338 + 4*
fDy1*(9*phixz - 9*fTAlph*phiyz - 56*
fDz*fTAlph*v.x()
346 G4cout <<
"coef = " << c[0] <<
" "
359 for (
G4int i = 0 ; i<num ; i++ )
364 G4cout <<
"Solution " << i <<
" : " << srd[i] <<
G4endl ;
366 phi = std::fmod(srd[i] , pihalf) ;
378 xbuf.push_back(xbuftmp) ;
381 G4cout <<
"solution " << i <<
" = " << phi <<
" , " << u <<
G4endl ;
396 for (
size_t k = 0 ; k<xbuf.size() ; k++ )
399 G4cout <<
"Solution " << k <<
" : "
400 <<
"reconstructed phiR = " << xbuf[k].phi
401 <<
", uR = " << xbuf[k].u <<
G4endl ;
407 IsConverged = false ;
409 for (
G4int i = 1 ; i<maxint ; i++ )
412 surfacenormal =
NormAng(phi,u) ;
415 deltaX = ( tmpxx - xxonsurface ).mag() ;
416 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
428 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist
429 <<
", " << deltaX <<
G4endl ;
437 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
440 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
444 if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0 ; }
447 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
461 if (tmpdist >= 0) tmpisvalid =
true;
469 if (tmpdist >= 0) { tmpisvalid =
true; }
474 G4Exception(
"G4TwistTrapAlphaSide::DistanceToSurface()",
476 "Feature NOT implemented !");
488 xbuf[k].distance = tmpdist ;
489 xbuf[k].areacode = tmpareacode ;
490 xbuf[k].isvalid = tmpisvalid ;
497 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
509 G4int nxxtmp = xbuf.size() ;
511 if ( nxxtmp<2 || IsParallel )
527 xbuf.push_back(xbuftmp) ;
542 xbuf.push_back(xbuftmp) ;
544 for (
size_t k = nxxtmp ; k<xbuf.size() ; k++ )
548 G4cout <<
"Solution " << k <<
" : "
549 <<
"reconstructed phiR = " << xbuf[k].phi
550 <<
", uR = " << xbuf[k].u <<
G4endl ;
556 IsConverged = false ;
558 for (
G4int i = 1 ; i<maxint ; i++ )
561 surfacenormal =
NormAng(phi,u) ;
563 deltaX = ( tmpxx - xxonsurface ).mag();
564 theta = std::fabs(std::acos(v*surfacenormal) - pihalf);
575 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist
576 <<
", " << deltaX << G4endl
577 <<
"X = " << tmpxx <<
G4endl ;
584 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
587 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
591 if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0; }
594 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
608 if (tmpdist >= 0) { tmpisvalid =
true; }
616 if (tmpdist >= 0) { tmpisvalid =
true; }
621 G4Exception(
"G4TwistedBoxSide::DistanceToSurface()",
623 "Feature NOT implemented !");
635 xbuf[k].distance = tmpdist ;
636 xbuf[k].areacode = tmpareacode ;
637 xbuf[k].isvalid = tmpisvalid ;
650 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
656 for (
size_t i = 0 ; i<xbuf.size() ; i++ )
658 distance[i] = xbuf[i].distance;
660 areacode[i] = xbuf[i].areacode ;
661 isvalid[i] = xbuf[i].isvalid ;
664 isvalid[i], nxx, validate, &gp, &gv);
666 G4cout <<
"element Nr. " << i
667 <<
", local Intersection = " << xbuf[i].xx
668 <<
", distance = " << xbuf[i].distance
669 <<
", u = " << xbuf[i].u
670 <<
", phi = " << xbuf[i].phi
671 <<
", isvalid = " << xbuf[i].isvalid
679 G4cout << nxx <<
" possible physical solutions found" <<
G4endl ;
680 for (
G4int k= 0 ; k< nxx ; k++ )
682 G4cout <<
"global intersection Point found: " << gxx[k] <<
G4endl ;
738 for (
G4int i = 1 ; i<20 ; i++ )
741 surfacenormal =
NormAng(phiR,uR) ;
743 deltaX = ( xx - xxonsurface ).mag() ;
746 G4cout <<
"i = " << i <<
", distance = " << distance[0]
747 <<
", " << deltaX <<
G4endl
748 <<
"X = " << xx <<
G4endl ;
755 if ( deltaX <= ctol ) { break ; }
762 if ( phiR > halfphi ) { phiR = halfphi ; }
763 if ( phiR < -halfphi ) { phiR = -halfphi ; }
764 if ( uR > uMax ) { uR = uMax ; }
765 if ( uR < -uMax ) { uR = -uMax ; }
768 distance[0] = ( p - xx ).mag() ;
769 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
774 G4cout <<
"refined solution " << phiR <<
" , " << uR <<
" , " <<
G4endl ;
783 G4cout <<
"intersection Point found: " << gxx[0] <<
G4endl ;
813 G4cout <<
"GetAreaCode: yprime = " << yprime <<
G4endl ;
814 G4cout <<
"Intervall is " << fYAxisMin <<
" to " << fYAxisMax <<
G4endl ;
829 if (yprime < fYAxisMin + ctol)
832 if (yprime <= fYAxisMin - ctol) { isoutside =
true; }
835 else if (yprime > fYAxisMax - ctol)
838 if (yprime >= fYAxisMax + ctol) { isoutside =
true; }
843 if (xx.z() <
fAxisMin[zaxis] + ctol)
852 if (xx.z() <=
fAxisMin[zaxis] - ctol) { isoutside =
true; }
854 else if (xx.z() >
fAxisMax[zaxis] - ctol)
862 if (xx.z() >=
fAxisMax[zaxis] + ctol) { isoutside =
true; }
871 areacode = tmpareacode;
883 if (yprime < fYAxisMin )
887 else if (yprime > fYAxisMax)
922 "Feature NOT implemented !");
992 "Method NOT implemented !");
1010 direction = direction.unit();
1016 direction = direction.unit();
1022 direction = direction.unit();
1028 direction = direction.unit();
1037 "Feature NOT implemented !");
1060 -
fPhiTwist*(fTAlph*p.x() + p.y())))*std::cos(phi)
1065 /
fDy1 - 4*std::sin(phi)))
1066 *(std::fabs(((
fa1md1 + 4*
fDy1*fTAlph)*std::cos(phi))
1067 /
fDy1 - 4*std::sin(phi)))
1068 + (std::fabs(4*std::cos(phi)
1070 * (std::fabs(4*std::cos(phi)
1130 for (
G4int i = 0 ; i<
n ; i++ )
1132 z = -
fDz+i*(2.*
fDz)/(n-1) ;
1136 for (
G4int j = 0 ; j<k ; j++ )
1138 nnode =
GetNode(i,j,k,n,iside) ;
1139 u = -b/2 +j*b/(k-1) ;
1142 xyz[nnode][0] = p.x() ;
1143 xyz[nnode][1] = p.y() ;
1144 xyz[nnode][2] = p.z() ;
1146 if ( i<n-1 && j<k-1 )
1148 nface =
GetFace(i,j,k,n,iside) ;
1150 * (
GetNode(i ,j ,k,n,iside)+1) ;
1152 * (
GetNode(i ,j+1,k,n,iside)+1) ;
1154 * (
GetNode(i+1,j+1,k,n,iside)+1) ;
1156 * (
GetNode(i+1,j ,k,n,iside)+1) ;
static const G4int sAxisZ
G4int GetAreacode(G4int i) const
G4ThreeVector ProjectPoint(const G4ThreeVector &p, G4bool isglobal=false)
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
virtual G4double GetBoundaryMax(G4double phi)
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
G4SurfCurNormal fCurrentNormal
virtual void SetCorners()
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
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
virtual ~G4TwistTrapAlphaSide()
G4ThreeVector NormAng(G4double phi, G4double u)
static double normal(HepRandomEngine *eptr)
static const G4int sC0Min1Min
G4TwistTrapAlphaSide(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)
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
G4double GetValueB(G4double phi)
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
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
static const G4int sAxisY
static const G4int sInside
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
static const G4int sCorner
static const G4double factor
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sAxisMax
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
virtual void SetBoundaries()
G4ThreeVector GetXX(G4int i) const
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
virtual G4double GetBoundaryMin(G4double phi)
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
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)