50 axis0min, axis1min, axis0max, axis1max),
55 G4Exception(
"G4TwistTubsSide::G4TwistTubsSide()",
"GeomSolids0002",
90 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
223 for (
auto i=0; i<2; ++i)
249 isvalid[0], 0, validate, &gp, &gv);
269 distance[0] = - c /
b;
270 xx[0] = p + distance[0]*
v;
278 if (distance[0] >= 0) isvalid[0] =
true;
286 if (distance[0] >= 0) isvalid[0] =
true;
295 if (distance[0] >= 0) isvalid[0] =
true;
301 areacode[0], isvalid[0],
302 0, validate, &gp, &gv);
308 isvalid[0], 1, validate, &gp, &gv);
321 isvalid[0], 0, validate, &gp, &gv);
333 G4bool tmpisvalid[2] = {
false,
false};
335 for (
auto i=0; i<2; ++i)
342 if ( b * D < 0 && std::fabs(bminusD / D) < protection )
345 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
349 tmpdist[i] = factor * bminusD;
353 tmpxx[i] = p + tmpdist[i]*
v;
360 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
369 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
376 if (tmpxx[i].
x() > 0)
379 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
388 if (tmpdist[0] <= tmpdist[1])
390 distance[0] = tmpdist[0];
391 distance[1] = tmpdist[1];
396 areacode[0] = tmpareacode[0];
397 areacode[1] = tmpareacode[1];
398 isvalid[0] = tmpisvalid[0];
399 isvalid[1] = tmpisvalid[1];
403 distance[0] = tmpdist[1];
404 distance[1] = tmpdist[0];
409 areacode[0] = tmpareacode[1];
410 areacode[1] = tmpareacode[0];
411 isvalid[0] = tmpisvalid[1];
412 isvalid[1] = tmpisvalid[0];
416 isvalid[0], 2, validate, &gp, &gv);
418 isvalid[1], 2, validate, &gp, &gv);
424 if (!isvalid[
k])
continue;
427 * xx[k].
z() , xx[k].
z());
428 G4double deltaY = (xx[
k] - xxonsurface).mag();
435 for (l=0; l<maxcount; ++l)
439 surfacenormal, xx[k]);
440 deltaY = (xx[
k] - xxonsurface).mag();
441 if (deltaY > lastdeltaY) { }
445 xxonsurface.
set(xx[k].
x(),
446 fKappa * std::fabs(xx[k].
x()) * xx[k].
z(),
452 message <<
"Exceeded maxloop count!" <<
G4endl
453 <<
" maxloop count " << maxcount;
454 G4Exception(
"G4TwistTubsFlatSide::DistanceToSurface(p,v)",
467 isvalid[0], 0, validate, &gp, &gv);
494 for (
auto i=0; i<2; ++i)
516 for (
auto i=0; i<2; ++i)
521 if ((gp - lastgxx[0]).mag() < halftol
522 || (gp - lastgxx[1]).mag() < halftol)
582 else if (p.
z() < C.
z())
593 else if (p.
z() < A.
z())
603 for (
auto i=0; i<2; ++i)
608 B = x0[i] + ((A.
z() - x0[i].
z()) / d[i].
z()) * d[i];
614 D = x0[i] + ((C.
z() - x0[i].
z()) / d[i].
z()) * d[i];
684 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol)
686 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
696 if (distToACB * distToCAD > 0 && distToACB < 0)
705 if (distToACB * distToCAD > 0)
709 if (distToACB <= distToCAD)
711 distance[0] = distToACB;
716 distance[0] = distToCAD;
726 distance[0] = distToACB;
731 distance[0] = distToCAD;
771 if (distToanm * distTocmn > 0 && distToanm < 0)
775 "Point p is behind the surfaces.");
779 if (std::fabs(distToanm) <= halftol)
785 else if (std::fabs(distTocmn) <= halftol)
792 if (distToanm <= distTocmn)
867 if (xx.
z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
870 else if (xx.
z() >
fAxisMax[zaxis] - ctol)
876 if (xx.
z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
885 areacode = tmpareacode;
932 "Feature NOT implemented !");
955 x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
956 y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
961 x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
962 y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
967 x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
968 y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
973 x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
974 y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
982 message <<
"Feature NOT implemented !" <<
G4endl
984 <<
" fAxis[1] = " <<
fAxis[1];
997 "Method NOT implemented !");
1013 direction = direction.
unit();
1019 direction = direction.
unit();
1025 direction = direction.
unit();
1031 direction = direction.
unit();
1039 message <<
"Feature NOT implemented !" <<
G4endl
1041 <<
" fAxis[1] = " <<
fAxis[1];
1063 for (
G4int i = 0 ; i<
n ; ++i )
1067 for (
G4int j = 0 ; j<
k ; ++j )
1069 nnode =
GetNode(i,j,k,n,iside) ;
1076 x = xmin + j*(xmax-
xmin)/(k-1) ;
1080 x = xmax - j*(xmax-
xmin)/(k-1) ;
1085 xyz[nnode][0] = p.
x() ;
1086 xyz[nnode][1] = p.
y() ;
1087 xyz[nnode][2] = p.
z() ;
1089 if ( i<n-1 && j<k-1 )
1091 nface =
GetFace(i,j,k,n,iside) ;
1094 * (
GetNode(i ,j ,k,n,iside)+1) ;
1096 * (
GetNode(i+1,j ,k,n,iside)+1) ;
1098 * (
GetNode(i+1,j+1,k,n,iside)+1) ;
1100 * (
GetNode(i ,j+1,k,n,iside)+1) ;