41 double pRmin1,
double pRmax1,
42 double pRmin2,
double pRmax2,
44 double pSPhi,
double pDPhi)
45 :
VUSolid(pName.c_str()), fRmin1(pRmin1), fRmin2(pRmin2),
46 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
54 std::ostringstream message;
55 message <<
"Invalid Z half-length for Solid: " <<
GetName() << std::endl
62 if (((pRmin1 >= pRmax1) || (pRmin2 >= pRmax2) || (pRmin1 < 0)) && (pRmin2 < 0))
64 std::ostringstream message;
65 message <<
"Invalid values of radii for Solid: " <<
GetName() << std::endl
66 <<
" pRmin1 = " << pRmin1 <<
", pRmin2 = " << pRmin2
67 <<
", pRmax1 = " << pRmax1 <<
", pRmax2 = " << pRmax2;
71 if ((pRmin1 == 0.0) && (pRmin2 > 0.0))
75 if ((pRmin2 == 0.0) && (pRmin1 > 0.0))
93 :
VUSolid(
""), kRadTolerance(0.), kAngTolerance(0.),
94 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
95 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
96 cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
115 :
VUSolid(rhs), kRadTolerance(rhs.kRadTolerance),
116 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
117 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
118 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
119 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
120 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
121 cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone)
141 VUSolid::operator=(rhs);
178 double distZ, distRMin, distRMax;
180 double pRMin, widRMin;
181 double pRMax, widRMax;
187 UVector3 nR, nr(0., 0., 0.), nPs, nPe;
189 distZ = std::fabs(std::fabs(p.
z) -
fDz);
190 rho = std::sqrt(p.
x * p.
x + p.
y * p.
y);
194 distRMin = std::fabs(pRMin - widRMin) /
secRMin;
198 distRMax = std::fabs(pRMax - widRMax) /
secRMax;
204 pPhi = std::atan2(p.
y, p.
x);
206 if (pPhi <
fSPhi - delta)
215 distSPhi = std::fabs(pPhi -
fSPhi);
235 if (distRMax <= delta)
247 if (distSPhi <= dAngle)
252 if (distEPhi <= dAngle)
275 Warning, 1,
"Point p is not on surface !?");
279 else if (noSurfaces == 1)
285 norm = sumnorm.
Unit();
290 return (
bool) noSurfaces;
303 double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin;
304 double pRMin, widRMin;
305 double pRMax, widRMax;
307 distZ = std::fabs(std::fabs(p.
z) -
fDz);
308 rho = std::sqrt(p.
x * p.
x + p.
y * p.
y);
312 distRMin = std::fabs(pRMin - widRMin) /
secRMin;
316 distRMax = std::fabs(pRMax - widRMax) /
secRMax;
318 if (distRMin < distRMax)
320 if (distZ < distRMin)
333 if (distZ < distRMax)
346 phi = std::atan2(p.
y, p.
x);
359 distSPhi = std::fabs(phi -
fSPhi) * rho;
362 distEPhi = std::fabs(phi -
fSPhi -
fDPhi) * rho;
366 if (distSPhi < distEPhi)
368 if (distSPhi < distMin)
375 if (distEPhi < distMin)
410 "Undefined side for valid surface normal to solid.");
448 double rMaxAv, rMaxOAv;
449 double rMinAv, rMinOAv;
452 double tolORMin, tolORMin2, tolIRMin, tolIRMin2;
453 double tolORMax2, tolIRMax, tolIRMax2;
454 double tolODz, tolIDz;
456 double Dist, sd, xi, yi, zi, ri = 0., risec, rhoi2, cosPsi;
458 double t1, t2, t3, b, c, d;
459 double nt1, nt2, nt3;
467 if (rMinAv > halfRadTolerance)
469 rMinOAv = rMinAv - halfRadTolerance;
476 rMaxOAv = rMaxAv + halfRadTolerance;
480 tolIDz =
fDz - halfCarTolerance;
481 tolODz =
fDz + halfCarTolerance;
483 if (std::fabs(p.
z) >= tolIDz)
487 sd = (std::fabs(p.
z) -
fDz) / std::fabs(v.
z);
496 rhoi2 = xi * xi + yi * yi ;
519 tolORMin2 = tolORMin * tolORMin;
520 tolIRMin2 = tolIRMin * tolIRMin;
529 tolIRMax2 = tolIRMax * tolIRMax;
536 if ((tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2))
580 t1 = 1.0 - v.
z * v.
z;
581 t2 = p.
x * v.
x + p.
y * v.
y;
582 t3 = p.
x * p.
x + p.
y * p.
y;
591 nt3 = t3 - rout * rout;
607 if ((rout < 0) && (nt3 <= 0))
614 sd = c / (-b - std::sqrt(d));
618 sd = -b + std::sqrt(d);
623 if ((b <= 0) && (c >= 0))
625 sd = c / (-b + std::sqrt(d));
631 sd = -b + std::sqrt(d);
644 double fTerm = sd - std::fmod(sd, dRmax);
649 if (std::fabs(zi) <= tolODz)
678 if ((t3 > (rin + halfRadTolerance *
secRMin)*
679 (rin + halfRadTolerance * secRMin))
680 && (nt2 < 0) && (d >= 0) && (std::fabs(p.
z) <= tolIDz))
687 risec = std::sqrt(xi * xi + yi * yi) *
secRMax;
694 if (norm.
Dot(v) <= 0)
702 if (norm.
Dot(v) <= 0)
714 sd = -0.5 * nt3 / nt2;
724 if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
766 nt3 = t3 - rin * rin;
782 sd = c / (-b - std::sqrt(d));
786 sd = -b + std::sqrt(d);
794 double fTerm = sd - std::fmod(sd, dRmax);
799 if (std::fabs(zi) <= tolODz)
810 if (sd > halfRadTolerance)
818 risec = std::sqrt(xi * xi + yi * yi) *
secRMin;
819 norm =
UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
820 if (norm.
Dot(v) <= 0)
829 if (sd > halfRadTolerance)
839 risec = std::sqrt(xi * xi + yi * yi) *
secRMin;
841 if (norm.
Dot(v) <= 0)
866 sd = c / (-b - std::sqrt(d));
870 sd = -b + std::sqrt(d);
877 if ((sd >= 0) && (std::fabs(zi) <= tolODz))
882 double fTerm = sd - std::fmod(sd, dRmax);
893 if (sd > halfRadTolerance)
901 risec = std::sqrt(xi * xi + yi * yi) *
secRMin;
902 norm =
UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
903 if (norm.
Dot(v) <= 0)
912 if (sd > halfRadTolerance)
922 risec = std::sqrt(xi * xi + yi * yi) *
secRMin;
923 norm =
UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
924 if (norm.
Dot(v) <= 0)
936 sd = -b - std::sqrt(d);
940 sd = c / (-b + std::sqrt(d));
945 if ((sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz))
950 double fTerm = sd - std::fmod(sd, dRmax);
961 if (sd > halfRadTolerance)
969 risec = std::sqrt(xi * xi + yi * yi) *
secRMin;
970 norm =
UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
971 if (norm.
Dot(v) <= 0)
980 if (sd > halfRadTolerance)
990 risec = std::sqrt(xi * xi + yi * yi) *
secRMin;
991 norm =
UVector3(-xi / risec, -yi / risec, tanRMin / secRMin);
992 if (norm.
Dot(v) <= 0)
1009 if (std::fabs(p.
z) <= tolODz)
1042 sd = -b - std::sqrt(d);
1046 sd = c / (-b + std::sqrt(d));
1048 zi = p.
z + sd * v.
z;
1055 sd = c / (-b - std::sqrt(d));
1059 sd = -b + std::sqrt(d);
1062 zi = p.
z + sd * v.
z;
1064 if ((sd >= 0) && (std::fabs(zi) <= tolODz))
1069 double fTerm = sd - std::fmod(sd, dRmax);
1074 xi = p.
x + sd * v.
x;
1075 yi = p.
y + sd * v.
y;
1107 sd = c / (-b - std::sqrt(d));
1111 sd = -b + std::sqrt(d);
1113 zi = p.
z + sd * v.
z;
1115 if ((sd >= 0) && (std::fabs(zi) <= tolODz))
1120 double fTerm = sd - std::fmod(sd, dRmax);
1125 xi = p.
x + sd * v.
x;
1126 yi = p.
y + sd * v.
y;
1165 if (Dist < halfCarTolerance)
1176 zi = p.
z + sd * v.
z;
1178 if (std::fabs(zi) <= tolODz)
1180 xi = p.
x + sd * v.
x;
1181 yi = p.
y + sd * v.
y;
1182 rhoi2 = xi * xi + yi * yi;
1186 if ((rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2))
1208 if (Dist < halfCarTolerance)
1219 zi = p.
z + sd * v.
z;
1221 if (std::fabs(zi) <= tolODz)
1223 xi = p.
x + sd * v.
x;
1224 yi = p.
y + sd * v.
y;
1225 rhoi2 = xi * xi + yi * yi;
1229 if ((rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2))
1244 if (snxt < halfCarTolerance)
1261 double safe = 0.0, rho, safeR1, safeR2, safeZ, safePhi;
1262 double pRMin, pRMax;
1265 rho = std::sqrt(p.
x * p.
x + p.
y * p.
y);
1266 safeZ = std::fabs(p.
z) -
fDz;
1267 safeR1 = 0; safeR2 = 0;
1272 safeR1 = (pRMin - rho) /
secRMin;
1275 safeR2 = (rho - pRMax) /
secRMax;
1277 if (safeR1 > safeR2)
1289 safe = (rho - pRMax) /
secRMax;
1301 if ((outside) && (safePhi > safe))
1308 safe = 0.0;
return safe;
1310 if (!aAccurate)
return safe;
1316 safsq += safeR1 * safeR1;
1321 safsq += safeR2 * safeR2;
1326 safsq += safeZ * safeZ;
1329 if (count == 1)
return safe;
1330 return std::sqrt(safsq);
1355 double snxt, srd, sphi, pdist;
1360 double t1, t2, t3, rout, rin, nt1, nt2, nt3;
1361 double b, c, d, sr2, sr3;
1370 double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi;
1371 double zi, ri, deltaRoi2;
1379 if (pdist > halfCarTolerance)
1395 if (pdist > halfCarTolerance)
1397 snxt = -pdist / v.
z;
1402 aNormalVector =
UVector3(0, 0, -1);
1432 t1 = 1.0 - v.
z * v.
z;
1433 t2 = p.
x * v.
x + p.
y * v.
y;
1434 t3 = p.
x * p.
x + p.
y * p.
y;
1439 nt3 = t3 - rout * rout;
1443 deltaRoi2 = snxt * snxt * t1 + 2 * snxt * t2 + t3
1448 deltaRoi2 = snxt * snxt * t1 + 2 * snxt * t2 + t3
1456 if (nt1 && (deltaRoi2 > 0.0))
1469 if (nt3 > -halfRadTolerance && nt2 >= 0)
1471 risec = std::sqrt(t3) *
secRMax;
1481 srd = -b - std::sqrt(d);
1485 srd = c / (-b + std::sqrt(d));
1488 zi = p.
z + srd * v.
z;
1491 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
1499 if ((ri < 0) || (srd < halfRadTolerance))
1506 sr2 = c / (-b - std::sqrt(d));
1510 sr2 = -b + std::sqrt(d);
1512 zi = p.
z + sr2 * v.
z;
1515 if ((ri >= 0) && (sr2 > halfRadTolerance))
1523 if ((-halfRadTolerance <= sr2) && (sr2 <= halfRadTolerance))
1540 risec = std::sqrt(t3) *
secRMax;
1546 else if (nt2 && (deltaRoi2 > 0.0))
1550 risec = std::sqrt(t3) *
secRMax;
1565 if (slentol <= halfCarTolerance)
1576 xi = p.
x + slentol * v.
x;
1577 yi = p.
y + slentol * v.
y;
1578 risec = std::sqrt(xi * xi + yi * yi) *
secRMax;
1581 if (norm.
Dot(v) > 0)
1583 aNormalVector = norm.
Unit();
1605 nt3 = t3 - rin * rin;
1623 risec = std::sqrt(p.
x * p.
x + p.
y * p.
y) *
secRMin;
1632 sr2 = -b - std::sqrt(d);
1636 sr2 = c / (-b + std::sqrt(d));
1638 zi = p.
z + sr2 * v.
z;
1641 if ((ri >= 0.0) && (-halfRadTolerance <= sr2) && (sr2 <= halfRadTolerance))
1649 if ((ri < 0) || (sr2 < halfRadTolerance))
1653 sr3 = c / (-b - std::sqrt(d));
1657 sr3 = -b + std::sqrt(d);
1663 if (sr3 > halfRadTolerance)
1667 zi = p.
z + sr3 * v.
z;
1677 else if (sr3 > -halfRadTolerance)
1685 else if ((sr2 < srd) && (sr2 > halfCarTolerance))
1690 else if (sr2 > -halfCarTolerance)
1697 if (slentol <= halfCarTolerance)
1706 xi = p.
x + slentol * v.
x;
1707 yi = p.
y + slentol * v.
y;
1708 if (sidetol ==
kRMax)
1710 risec = std::sqrt(xi * xi + yi * yi) *
secRMax;
1715 risec = std::sqrt(xi * xi + yi * yi) *
secRMin;
1718 if (norm.
Dot(v) > 0)
1722 aNormalVector = norm.
Unit();
1748 vphi = std::atan2(v.
y, v.
x);
1750 if (vphi <
fSPhi - halfAngTolerance)
1754 else if (vphi >
fSPhi +
fDPhi + halfAngTolerance)
1774 && (pDistE <= halfCarTolerance)))
1776 && (pDistE > halfCarTolerance))))
1781 sphi = pDistS / compS;
1782 if (sphi >= -halfCarTolerance)
1784 xi = p.
x + sphi * v.
x;
1785 yi = p.
y + sphi * v.
y;
1794 if ((
fSPhi - halfAngTolerance <= vphi)
1795 && (
fSPhi +
fDPhi + halfAngTolerance >= vphi))
1807 if (pDistS > -halfCarTolerance)
1825 sphi2 = pDistE / compE;
1829 if ((sphi2 > -halfCarTolerance) && (sphi2 < sphi))
1831 xi = p.
x + sphi2 * v.
x;
1832 yi = p.
y + sphi2 * v.
y;
1841 if (!((
fSPhi - halfAngTolerance <= vphi)
1842 && (
fSPhi +
fDPhi + halfAngTolerance >= vphi)))
1845 if (pDistE <= -halfCarTolerance)
1861 if (pDistE <= -halfCarTolerance)
1883 if ((
fSPhi - halfAngTolerance <= vphi)
1884 && (vphi <=
fSPhi +
fDPhi + halfAngTolerance))
1910 xi = p.
x + snxt * v.
x;
1911 yi = p.
y + snxt * v.
y;
1912 risec = std::sqrt(xi * xi + yi * yi) *
secRMax;
1917 xi = p.
x + snxt * v.
x;
1918 yi = p.
y + snxt * v.
y;
1919 risec = std::sqrt(xi * xi + yi * yi) *
secRMin;
1952 aNormalVector =
UVector3(0, 0, -1);
1958 std::ostringstream message;
1959 int oldprc = message.precision(16);
1960 message <<
"Undefined side for valid surface normal to solid."
1962 <<
"Position:" << std::endl << std::endl
1963 <<
"p.x = " << p.
x <<
" mm" << std::endl
1964 <<
"p.y = " << p.
y <<
" mm" << std::endl
1965 <<
"p.z = " << p.
z <<
" mm" << std::endl << std::endl
1966 <<
"pho at z = " << std::sqrt(p.
x * p.
x + p.
y * p.
y)
1967 <<
" mm" << std::endl << std::endl;
1968 if (p.
x != 0. || p.
y != 0.)
1970 message <<
"point phi = " << std::atan2(p.
y, p.
x) / (
UUtils::kPi / 180.0)
1971 <<
" degree" << std::endl << std::endl;
1973 message <<
"Direction:" << std::endl << std::endl
1974 <<
"v.x = " << v.
x << std::endl
1975 <<
"v.y = " << v.
y << std::endl
1976 <<
"v.z = " << v.
z << std::endl << std::endl
1977 <<
"Proposed distance :" << std::endl << std::endl
1978 <<
"snxt = " << snxt <<
" mm" << std::endl;
1979 message.precision(oldprc);
1983 if (snxt < halfCarTolerance)
1997 double safe = 0.0, rho, safeZ;
2002 int oldprc = cout.precision(16);
2005 cout <<
"Position:" << std::endl << std::endl;
2006 cout <<
"p.x = " << p.
x <<
" mm" << std::endl;
2007 cout <<
"p.y = " << p.
y <<
" mm" << std::endl;
2008 cout <<
"p.z = " << p.
z <<
" mm" << std::endl << std::endl;
2009 cout <<
"pho at z = " << std::sqrt(p.
x * p.
x + p.
y * p.
y)
2010 <<
" mm" << std::endl << std::endl;
2011 if ((p.
x != 0.) || (p.
x != 0.))
2013 cout <<
"point phi = " << std::atan2(p.
y, p.
x) /
degree
2014 <<
" degree" << std::endl << std::endl;
2016 cout.precision(oldprc);
2022 rho = std::sqrt(p.
x * p.
x + p.
y * p.
y);
2025 safeZ =
fDz - std::fabs(p.
z);
2046 return std::string(
"Cons");
2056 return new UCons(*
this);
2065 int oldprc = os.precision(16);
2066 os <<
"-----------------------------------------------------------\n"
2067 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
2068 <<
" ===================================================\n"
2069 <<
" Solid type: UCons\n"
2070 <<
" Parameters: \n"
2071 <<
" inside -fDz radius : " <<
fRmin1 <<
" mm \n"
2072 <<
" outside -fDz radius: " <<
fRmax1 <<
" mm \n"
2073 <<
" inside +fDz radius : " <<
fRmin2 <<
" mm \n"
2074 <<
" outside +fDz radius: " <<
fRmax2 <<
" mm \n"
2075 <<
" half length in Z : " <<
fDz <<
" mm \n"
2076 <<
" starting angle of segment: " <<
fSPhi / (
UUtils::kPi / 180.0) <<
" degrees \n"
2077 <<
" delta angle of segment : " <<
fDPhi / (
UUtils::kPi / 180.0) <<
" degrees \n"
2078 <<
"-----------------------------------------------------------\n";
2079 os.precision(oldprc);
2094 double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
2095 double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
2117 cosu = std::cos(phi);
2118 sinu = std::sin(phi);
2126 chose =
UUtils::Random(0., Aone + Atwo + Athree + Afour + 2.*Afive);
2128 if ((chose >= 0.) && (chose < Aone))
2133 return UVector3(rtwo * cosu * (qtwo - zRand),
2134 rtwo * sinu * (qtwo - zRand), zRand);
2142 else if ((chose >= Aone) && (chose <= Aone + Atwo))
2147 return UVector3(rone * cosu * (qone - zRand),
2148 rone * sinu * (qone - zRand), zRand);
2156 else if ((chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree))
2158 return UVector3(rRand1 * cosu, rRand1 * sinu, -1 *
fDz);
2160 else if ((chose >= Aone + Atwo + Athree)
2161 && (chose < Aone + Atwo + Athree + Afour))
2163 return UVector3(rRand2 * cosu, rRand2 * sinu,
fDz);
2165 else if ((chose >= Aone + Atwo + Athree + Afour)
2166 && (chose < Aone + Atwo + Athree + Afour + Afive))
2172 rRand1 * std::sin(
fSPhi), zRand);
double SafetyToPhi(const UVector3 &p, const double rho, bool &outside) const
double SafetyFromOutside(const UVector3 &p, bool precise) const
static double frTolerance
VUSolid::EnumInside Inside(const UVector3 &p) const
UGeometryType GetEntityType() const
std::ostream & StreamInfo(std::ostream &os) const
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
UCons & operator=(const UCons &rhs)
const std::string & GetName() const
double SafetyFromInsideR(const UVector3 &p, const double rho, bool) const
double SafetyFromInside(const UVector3 &p, bool precise) const
double GetOuterRadiusPlusZ() const
static double Tolerance()
static const double kInfinity
double GetStartPhiAngle() const
bool Normal(const UVector3 &p, UVector3 &n) const
double GetZHalfLength() const
static double faTolerance
double Dot(const UVector3 &) const
double GetInnerRadiusPlusZ() const
double GetRadiusInRing(double rmin, double rmax)
double GetOuterRadiusMinusZ() const
void Extent(UVector3 &aMin, UVector3 &aMax) const
virtual void GetParametersList(int, double *) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void CheckPhiAngles(double sPhi, double dPhi)
double GetDeltaPhiAngle() const
static const double degree
std::string UGeometryType
double GetInnerRadiusMinusZ() const
UVector3 GetPointOnSurface() const
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
UVector3 ApproxSurfaceNormal(const UVector3 &p) const
double Random(double min=0.0, double max=1.0)
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const