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.),
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] <<
" "
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 ;
452 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
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 ;
599 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
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 ;
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;
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 )
1040 z = -
fDz+i*(2.*
fDz)/(n-1) ;
1045 for (
G4int j = 0 ; j<
k ; ++j )
1047 nnode =
GetNode(i,j,k,n,iside) ;
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 )
1057 nface =
GetFace(i,j,k,n,iside) ;
1059 * (
GetNode(i ,j ,k,n,iside)+1) ;
1061 * (
GetNode(i ,j+1,k,n,iside)+1) ;
1063 * (
GetNode(i+1,j+1,k,n,iside)+1) ;
1065 * (
GetNode(i+1,j ,k,n,iside)+1) ;