111 :
G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
112 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
113 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.),
114 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.),
194 G4bool IsParallel = false ;
195 G4bool IsConverged = false ;
238 G4bool tmpisvalid = false ;
240 std::vector<Intersection> xbuf ;
254 if ( std::fabs(p.
z()) <=
L )
260 + 2*
fDy2minus1*phi)*(v.
x()*std::cos(phi) + v.
y()*std::sin(phi)))
261 / (2.*
fPhiTwist*(v.
y()*std::cos(phi) - v.
x()*std::sin(phi)));
269 xbuf.push_back(xbuftmp) ;
278 areacode[0], isvalid[0],
279 0, validate, &gp, &gv);
290 c[6] = 120*(-52*phiyz - 120*
fDz*v.
x() + 60*
fdeltaX*v.
z()
294 c[4] = 12*(127*phiyz + 640*
fDz*v.
x() - 320*
fdeltaX*v.
z()
304 G4cout <<
"coef = " << c[0] <<
" "
318 for (
G4int i = 0 ; i<num ; ++i )
323 G4cout <<
"Solution " << i <<
" : " << srd[i] <<
G4endl ;
325 phi = std::fmod(srd[i] , pihalf) ;
326 u = (1/std::cos(phi)*(2*phixz + 4*
fDz*phi*v.
x()
336 xbuf.push_back(xbuftmp) ;
339 G4cout <<
"solution " << i <<
" = " << phi <<
" , " << u <<
G4endl ;
356 for (
size_t k = 0 ; k<xbuf.size() ; ++k )
359 G4cout <<
"Solution " << k <<
" : "
360 <<
"reconstructed phiR = " << xbuf[k].phi
361 <<
", uR = " << xbuf[k].u <<
G4endl ;
367 IsConverged = false ;
369 for (
G4int i = 1 ; i<maxint ; ++i )
372 surfacenormal =
NormAng(phi,u) ;
374 deltaX = ( tmpxx - xxonsurface ).mag() ;
375 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
387 G4cout <<
"Step i = " << i <<
", distance = "
388 << tmpdist <<
", " << deltaX <<
G4endl ;
395 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
398 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
401 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
404 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
418 if (tmpdist >= 0) tmpisvalid =
true;
426 if (tmpdist >= 0) tmpisvalid =
true;
431 G4Exception(
"G4TwistTrapParallelSide::DistanceToSurface()",
433 "Feature NOT implemented !");
444 xbuf[k].distance = tmpdist ;
445 xbuf[k].areacode = tmpareacode ;
446 xbuf[k].isvalid = tmpisvalid ;
462 G4int nxxtmp = xbuf.size() ;
464 if ( nxxtmp<2 || IsParallel )
480 xbuf.push_back(xbuftmp) ;
495 xbuf.push_back(xbuftmp) ;
497 for (
size_t k = nxxtmp ; k<xbuf.size() ; ++k )
500 G4cout <<
"Solution " << k <<
" : "
501 <<
"reconstructed phiR = " << xbuf[k].phi
502 <<
", uR = " << xbuf[k].u <<
G4endl ;
508 IsConverged = false ;
510 for (
G4int i = 1 ; i<maxint ; ++i )
513 surfacenormal =
NormAng(phi,u) ;
515 deltaX = ( tmpxx - xxonsurface ).mag() ;
516 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
527 G4cout <<
"Step i = " << i <<
", distance = "
528 << tmpdist <<
", " << deltaX <<
G4endl ;
535 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
538 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
541 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
544 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
558 if (tmpdist >= 0) tmpisvalid =
true;
566 if (tmpdist >= 0) tmpisvalid =
true;
571 G4Exception(
"G4TwistedBoxSide::DistanceToSurface()",
573 "Feature NOT implemented !");
585 xbuf[k].distance = tmpdist ;
586 xbuf[k].areacode = tmpareacode ;
587 xbuf[k].isvalid = tmpisvalid ;
605 for (
size_t i = 0 ; i<xbuf.size() ; ++i )
607 distance[i] = xbuf[i].distance;
609 areacode[i] = xbuf[i].areacode ;
610 isvalid[i] = xbuf[i].isvalid ;
613 isvalid[i], nxx, validate, &gp, &gv);
615 G4cout <<
"element Nr. " << i
616 <<
", local Intersection = " << xbuf[i].xx
617 <<
", distance = " << xbuf[i].distance
618 <<
", u = " << xbuf[i].u
619 <<
", phi = " << xbuf[i].phi
620 <<
", isvalid = " << xbuf[i].isvalid
626 G4cout <<
"G4TwistTrapParallelSide finished " <<
G4endl ;
627 G4cout << nxx <<
" possible physical solutions found" <<
G4endl ;
628 for (
G4int k= 0 ; k< nxx ; ++k )
630 G4cout <<
"global intersection Point found: " << gxx[k] <<
G4endl ;
684 for (
G4int i = 1 ; i<maxint ; ++i )
687 surfacenormal =
NormAng(phiR,uR) ;
689 deltaX = ( xx - xxonsurface ).mag() ;
692 G4cout <<
"i = " << i <<
", distance = "
693 << distance[0] <<
", " << deltaX <<
G4endl ;
700 if ( deltaX <= ctol ) { break ; }
709 if ( phiR > halfphi ) phiR = halfphi ;
710 if ( phiR < -halfphi ) phiR = -halfphi ;
711 if ( uR > uMax ) uR = uMax ;
712 if ( uR < uMin ) uR = uMin ;
715 distance[0] = ( p - xx ).mag() ;
716 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
721 G4cout <<
"refined solution " << phiR <<
" , " << uR <<
" , " <<
G4endl ;
730 G4cout <<
"intersection Point found: " << gxx[0] <<
G4endl ;
759 G4cout <<
"GetAreaCode: yprime = " << yprime <<
G4endl ;
760 G4cout <<
"Intervall is " << fXAxisMin <<
" to " << fXAxisMax <<
G4endl ;
775 if (yprime < fXAxisMin + ctol)
778 if (yprime <= fXAxisMin - ctol) isoutside =
true;
781 else if (yprime > fXAxisMax - ctol)
784 if (yprime >= fXAxisMax + ctol) isoutside =
true;
795 if (xx.
z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
798 else if (xx.
z() >
fAxisMax[zaxis] - ctol)
804 if (xx.
z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
812 G4int tmpareacode = areacode & (~sInside);
813 areacode = tmpareacode;
825 if (yprime < fXAxisMin )
829 else if (yprime > fXAxisMax)
859 G4Exception(
"G4TwistTrapParallelSide::GetAreaCode()",
861 "Feature NOT implemented !");
918 G4Exception(
"G4TwistTrapParallelSide::SetCorners()",
920 "Method NOT implemented !");
938 direction = direction.
unit();
944 direction = direction.
unit();
950 direction = direction.
unit();
956 direction = direction.
unit();
962 G4Exception(
"G4TwistTrapParallelSide::SetCorners()",
964 "Feature NOT implemented !");
1038 for (
G4int i = 0 ; i<
n ; ++i )
1045 for (
G4int j = 0 ; j<k ; ++j )
1048 u = umax - j*(umax-umin)/(k-1) ;
1051 xyz[nnode][0] = p.
x() ;
1052 xyz[nnode][1] = p.
y() ;
1053 xyz[nnode][2] = p.
z() ;
1055 if ( i<
n-1 && j<k-1 )
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
static constexpr double L
static constexpr double pi
G4bool DistanceSort(const Intersection &a, const Intersection &b)
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
G4GLOB_DLL std::ostream G4cout
void set(double x, double y, double z)
HepRotation inverse() const
HepRotation & rotateZ(double delta)
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
virtual G4double GetBoundaryMin(G4double phi)
virtual G4double GetBoundaryMax(G4double phi)
virtual void SetBoundaries()
G4ThreeVector NormAng(G4double phi, G4double u)
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
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)
virtual ~G4TwistTrapParallelSide()
virtual void SetCorners()
G4ThreeVector ProjectPoint(const G4ThreeVector &p, G4bool isglobal=false)
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
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)
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)
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
static const G4int sAxisMin
static const G4int sC0Max1Max
static const G4int sAxis1
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
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
static const G4double kInfinity
static double normal(HepRandomEngine *eptr)
const char * name(G4int ptype)