86 fTAlph = std::tan(fAlph) ;
92 fDx4plus2 = fDx4 + fDx2 ;
93 fDx4minus2 = fDx4 - fDx2 ;
94 fDx3plus1 = fDx3 + fDx1 ;
95 fDx3minus1 = fDx3 - fDx1 ;
96 fDy2plus1 = fDy2 + fDy1 ;
97 fDy2minus1 = fDy2 - fDy1 ;
99 fa1md1 = 2*fDx2 - 2*fDx1 ;
100 fa2md2 = 2*fDx4 - 2*fDx3 ;
102 fPhiTwist = PhiTwist ;
103 fAngleSide = AngleSide ;
105 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
106 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
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.),
166 GetPhiUAtX(xx,phi,u) ;
200 G4bool IsParallel = false ;
201 G4bool IsConverged = false ;
221 distance[i] = kInfinity;
224 gxx[i].
set(kInfinity, kInfinity, kInfinity);
243 G4bool tmpisvalid = false ;
245 std::vector<Intersection> xbuf ;
252 G4double phixz = fPhiTwist * ( p.
x() * v.
z() - p.
z() * v.
x() ) ;
253 G4double phiyz = fPhiTwist * ( p.
y() * v.
z() - p.
z() * v.
y() ) ;
259 if ( std::fabs(p.
z()) <= L ) {
261 phi = p.
z() * fPhiTwist / L ;
263 u = (2*(fdeltaY*phi*v.
x() - fPhiTwist*p.
y()*v.
x() - fdeltaX*phi*v.
y()
264 + fPhiTwist*p.
x()*v.
y()) + (fDy2plus1*fPhiTwist
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) ;
280 distance[0] = kInfinity;
281 gxx[0].
set(kInfinity,kInfinity,kInfinity);
285 areacode[0], isvalid[0],
286 0, validate, &gp, &gv);
302 c[8] = -3600*(-2*phiyz + fDy2plus1*fPhiTwist*v.
z()) ;
303 c[7] = -7200*(phixz - 2*fDz*v.
y() + (fdeltaY + fDy2minus1)*v.
z()) ;
304 c[6] = 120*(-52*phiyz - 120*fDz*v.
x() + 60*fdeltaX*v.
z() + 11*fDy2plus1*fPhiTwist*v.
z()) ;
305 c[5] = 240*(16*phixz - 52*fDz*v.
y() + 26*fdeltaY*v.
z() + 11*fDy2minus1*v.
z()) ;
306 c[4] = 12*(127*phiyz + 640*fDz*v.
x() - 320*fdeltaX*v.
z() + 4*fDy2plus1*fPhiTwist*v.
z()) ;
307 c[3] = -404*phixz + 3048*fDz*v.
y() - 1524*fdeltaY*v.
z() + 96*fDy2minus1*v.
z() ;
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()) ;
310 c[0] = 24*fDz*v.
x() - 12*fdeltaX*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) ;
336 u = (1/std::cos(phi)*(2*phixz + 4*fDz*phi*v.
x() - 2*fdeltaX*phi*v.
z() + (fDy2plus1*fPhiTwist + 2*fDy2minus1*phi)*v.
z()* std::sin(phi)))/(2.*fPhiTwist*v.
z()) ;
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++ ) {
381 xxonsurface = SurfacePoint(phi,u) ;
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 ;
399 GetPhiUAtX(tmpxx, phi, u) ;
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 ;
425 tmpareacode = GetAreaCode(tmpxx);
427 if (tmpdist >= 0) tmpisvalid =
true;
430 tmpareacode = GetAreaCode(tmpxx,
false);
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++ ) {
522 xxonsurface = SurfacePoint(phi,u) ;
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 ;
539 GetPhiUAtX(tmpxx, phi, u) ;
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 ;
565 tmpareacode = GetAreaCode(tmpxx);
567 if (tmpdist >= 0) tmpisvalid =
true;
570 tmpareacode = GetAreaCode(tmpxx,
false);
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 ;
678 distance[i] = kInfinity;
680 gxx[i].
set(kInfinity, kInfinity, kInfinity);
697 for (
G4int i = 1 ; i<maxint ; i++ ) {
699 xxonsurface = SurfacePoint(phiR,uR) ;
700 surfacenormal = NormAng(phiR,uR) ;
702 deltaX = ( xx - xxonsurface ).mag() ;
705 G4cout <<
"i = " << i <<
", distance = " << distance[0] <<
", " << deltaX <<
G4endl ;
710 GetPhiUAtX(xx, phiR, uR) ;
712 if ( deltaX <= ctol ) { break ; }
719 G4double uMax = GetBoundaryMax(phiR) ;
720 G4double uMin = GetBoundaryMin(phiR) ;
722 if ( phiR > halfphi ) phiR = halfphi ;
723 if ( phiR < -halfphi ) phiR = -halfphi ;
724 if ( uR > uMax ) uR = uMax ;
725 if ( uR < uMin ) uR = uMin ;
727 xxonsurface = SurfacePoint(phiR,uR) ;
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 ;
766 GetPhiUAtX(xx, phi,yprime ) ;
768 G4double fXAxisMax = GetBoundaryMax(phi) ;
769 G4double fXAxisMin = GetBoundaryMin(phi) ;
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;
805 if (xx.
z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
807 }
else if (xx.
z() >
fAxisMax[zaxis] - ctol) {
810 if (areacode & sBoundary) areacode |=
sCorner;
812 if (xx.
z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
820 areacode = tmpareacode;
821 }
else if ((areacode & sBoundary) !=
sBoundary) {
829 if (yprime < fXAxisMin ) {
831 }
else if (yprime > fXAxisMax) {
839 if (areacode & sBoundary) areacode |=
sCorner;
844 if (areacode & sBoundary) areacode |=
sCorner;
848 if ((areacode & sBoundary) !=
sBoundary) {
854 G4Exception(
"G4TwistTrapParallelSide::GetAreaCode()",
856 "Feature NOT implemented !");
864 void G4TwistTrapParallelSide::SetCorners()
875 x = -fdeltaX/2. + (-fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
876 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) + (fDx2 - fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
883 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
884 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
890 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
891 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
897 x = fdeltaX/2. + (-fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
898 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (-fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
905 G4Exception(
"G4TwistTrapParallelSide::SetCorners()",
907 "Method NOT implemented !");
914 void G4TwistTrapParallelSide::SetBoundaries()
925 direction = direction.
unit();
931 direction = direction.
unit();
937 direction = direction.
unit();
943 direction = direction.
unit();
949 G4Exception(
"G4TwistTrapParallelSide::SetCorners()",
951 "Feature NOT implemented !");
967 phi = p.
z()/(2*fDz)*fPhiTwist ;
969 u = ((-(fdeltaX*phi) + fPhiTwist*p.
x())* std::cos(phi) + (-(fdeltaY*phi) + fPhiTwist*p.
y())*std::sin(phi))/fPhiTwist ;
990 GetPhiUAtX( tmpp, phi, u ) ;
1021 for ( i = 0 ; i<
n ; i++ ) {
1023 z = -fDz+i*(2.*fDz)/(n-1) ;
1024 phi = z*fPhiTwist/(2*fDz) ;
1025 umin = GetBoundaryMin(phi) ;
1026 umax = GetBoundaryMax(phi) ;
1028 for ( j = 0 ; j<k ; j++ ) {
1030 nnode =
GetNode(i,j,k,n,iside) ;
1031 u = umax - j*(umax-umin)/(k-1) ;
1032 p = SurfacePoint(phi,u,
true) ;
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
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
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 ~G4TwistTrapParallelSide()
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
static const G4int sAxisX
HepRotation inverse() const
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
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
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
G4ThreeVector GetXX(G4int i) 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)
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)