108 fRot.rotateZ( AngleSide ) ;
123 :
G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
124 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
125 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.),
126 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.),
200 G4bool IsParallel = false ;
201 G4bool IsConverged = false ;
243 G4bool tmpisvalid = false ;
245 std::vector<Intersection> xbuf ;
259 if ( std::fabs(p.z()) <= L ) {
265 + 2*
fDy2minus1*phi)*(v.x()*std::cos(phi) + v.y()*std::sin(phi)))
266 / (2.*
fPhiTwist*(v.y()*std::cos(phi) - v.x()*std::sin(phi)));
274 xbuf.push_back(xbuftmp) ;
285 areacode[0], isvalid[0],
286 0, validate, &gp, &gv);
308 c[2] = -72*phiyz + 404*(-2*
fDz*v.x() +
fdeltaX*v.z()) ;
309 c[1] = 12*(phixz - 12*
fDz*v.y() + 6*
fdeltaY*v.z()) ;
314 G4cout <<
"coef = " << c[0] <<
" "
329 for (
G4int i = 0 ; i<num ; i++ ) {
332 G4cout <<
"Solution " << i <<
" : " << srd[i] <<
G4endl ;
334 phi = std::fmod(srd[i] , pihalf) ;
344 xbuf.push_back(xbuftmp) ;
347 G4cout <<
"solution " << i <<
" = " << phi <<
" , " << u <<
G4endl ;
366 for (
size_t k = 0 ; k<xbuf.size() ; k++ ) {
369 G4cout <<
"Solution " << k <<
" : "
370 <<
"reconstructed phiR = " << xbuf[k].phi
371 <<
", uR = " << xbuf[k].u <<
G4endl ;
377 IsConverged = false ;
379 for (
G4int i = 1 ; i<maxint ; i++ ) {
382 surfacenormal =
NormAng(phi,u) ;
384 deltaX = ( tmpxx - xxonsurface ).mag() ;
385 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
386 if ( theta < 0.001 ) {
395 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist <<
", " << deltaX <<
G4endl ;
402 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
405 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
411 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
415 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
427 if (tmpdist >= 0) tmpisvalid =
true;
432 if (tmpdist >= 0) tmpisvalid =
true;
435 G4Exception(
"G4TwistTrapParallelSide::DistanceToSurface()",
437 "Feature NOT implemented !");
449 xbuf[k].distance = tmpdist ;
450 xbuf[k].areacode = tmpareacode ;
451 xbuf[k].isvalid = tmpisvalid ;
460 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
466 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() ,
EqualIntersection ) , xbuf.end() ) ;
471 G4int nxxtmp = xbuf.size() ;
473 if ( nxxtmp<2 || IsParallel ) {
489 xbuf.push_back(xbuftmp) ;
505 xbuf.push_back(xbuftmp) ;
507 for (
size_t k = nxxtmp ; k<xbuf.size() ; k++ ) {
510 G4cout <<
"Solution " << k <<
" : "
511 <<
"reconstructed phiR = " << xbuf[k].phi
512 <<
", uR = " << xbuf[k].u <<
G4endl ;
518 IsConverged = false ;
520 for (
G4int i = 1 ; i<maxint ; i++ ) {
523 surfacenormal =
NormAng(phi,u) ;
525 deltaX = ( tmpxx - xxonsurface ).mag() ;
526 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
527 if ( theta < 0.001 ) {
535 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist <<
", " << deltaX <<
G4endl ;
542 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
545 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
551 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
555 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
567 if (tmpdist >= 0) tmpisvalid =
true;
572 if (tmpdist >= 0) tmpisvalid =
true;
575 G4Exception(
"G4TwistedBoxSide::DistanceToSurface()",
577 "Feature NOT implemented !");
589 xbuf[k].distance = tmpdist ;
590 xbuf[k].areacode = tmpareacode ;
591 xbuf[k].isvalid = tmpisvalid ;
604 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() ,
EqualIntersection ) , xbuf.end() ) ;
607 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
613 for (
size_t i = 0 ; i<xbuf.size() ; i++ ) {
615 distance[i] = xbuf[i].distance;
617 areacode[i] = xbuf[i].areacode ;
618 isvalid[i] = xbuf[i].isvalid ;
621 isvalid[i], nxx, validate, &gp, &gv);
624 G4cout <<
"element Nr. " << i
625 <<
", local Intersection = " << xbuf[i].xx
626 <<
", distance = " << xbuf[i].distance
627 <<
", u = " << xbuf[i].u
628 <<
", phi = " << xbuf[i].phi
629 <<
", isvalid = " << xbuf[i].isvalid
637 G4cout <<
"G4TwistTrapParallelSide finished " <<
G4endl ;
638 G4cout << nxx <<
" possible physical solutions found" <<
G4endl ;
639 for (
G4int k= 0 ; k< nxx ; k++ ) {
640 G4cout <<
"global intersection Point found: " << gxx[k] <<
G4endl ;
697 for (
G4int i = 1 ; i<maxint ; i++ ) {
700 surfacenormal =
NormAng(phiR,uR) ;
702 deltaX = ( xx - xxonsurface ).mag() ;
705 G4cout <<
"i = " << i <<
", distance = " << distance[0] <<
", " << deltaX <<
G4endl ;
712 if ( deltaX <= ctol ) { break ; }
722 if ( phiR > halfphi ) phiR = halfphi ;
723 if ( phiR < -halfphi ) phiR = -halfphi ;
724 if ( uR > uMax ) uR = uMax ;
725 if ( uR < uMin ) uR = uMin ;
728 distance[0] = ( p - xx ).mag() ;
729 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
734 G4cout <<
"refined solution " << phiR <<
" , " << uR <<
" , " <<
G4endl ;
743 G4cout <<
"intersection Point found: " << gxx[0] <<
G4endl ;
773 G4cout <<
"GetAreaCode: yprime = " << yprime <<
G4endl ;
774 G4cout <<
"Intervall is " << fXAxisMin <<
" to " << fXAxisMax <<
G4endl ;
789 if (yprime < fXAxisMin + ctol) {
791 if (yprime <= fXAxisMin - ctol) isoutside =
true;
793 }
else if (yprime > fXAxisMax - ctol) {
795 if (yprime >= fXAxisMax + ctol) isoutside =
true;
800 if (xx.z() <
fAxisMin[zaxis] + ctol) {
805 if (xx.z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
807 }
else if (xx.z() >
fAxisMax[zaxis] - ctol) {
812 if (xx.z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
820 areacode = tmpareacode;
829 if (yprime < fXAxisMin ) {
831 }
else if (yprime > fXAxisMax) {
842 }
else if (xx.z() >
fAxisMax[zaxis]) {
854 G4Exception(
"G4TwistTrapParallelSide::GetAreaCode()",
856 "Feature NOT implemented !");
905 G4Exception(
"G4TwistTrapParallelSide::SetCorners()",
907 "Method NOT implemented !");
925 direction = direction.unit();
931 direction = direction.unit();
937 direction = direction.unit();
943 direction = direction.unit();
949 G4Exception(
"G4TwistTrapParallelSide::SetCorners()",
951 "Feature NOT implemented !");
1021 for ( i = 0 ; i<
n ; i++ ) {
1023 z = -
fDz+i*(2.*
fDz)/(n-1) ;
1028 for ( j = 0 ; j<k ; j++ ) {
1030 nnode =
GetNode(i,j,k,n,iside) ;
1031 u = umax - j*(umax-umin)/(k-1) ;
1034 xyz[nnode][0] = p.x() ;
1035 xyz[nnode][1] = p.y() ;
1036 xyz[nnode][2] = p.z() ;
1038 if ( i<n-1 && j<k-1 ) {
1040 nface =
GetFace(i,j,k,n,iside) ;
static const G4int sAxisZ
G4int GetAreacode(G4int i) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4SurfCurNormal fCurrentNormal
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sOutside
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
virtual ~G4TwistTrapParallelSide()
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
virtual G4double GetBoundaryMin(G4double phi)
static const G4int sAxisX
static double normal(HepRandomEngine *eptr)
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
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
virtual G4double GetBoundaryMax(G4double phi)
G4TwistTrapParallelSide(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)
static const G4int sInside
G4ThreeVector NormAng(G4double phi, G4double u)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4ThreeVector ProjectPoint(const G4ThreeVector &p, G4bool isglobal=false)
static const G4int sCorner
static const G4double factor
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
virtual void SetCorners()
static const G4int sAxisMax
const G4double x[NPOINTSGL]
G4ThreeVector GetXX(G4int i) const
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
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)
virtual void SetBoundaries()
CurrentStatus fCurStatWithV
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)