Gabriele Cosmo
937 lines
29 KiB
C++
937 lines
29 KiB
C++
//
|
|
// ********************************************************************
|
|
// * License and Disclaimer *
|
|
// * *
|
|
// * The Geant4 software is copyright of the Copyright Holders of *
|
|
// * the Geant4 Collaboration. It is provided under the terms and *
|
|
// * conditions of the Geant4 Software License, included in the file *
|
|
// * LICENSE and available at http://cern.ch/geant4/license . These *
|
|
// * include a list of copyright holders. *
|
|
// * *
|
|
// * Neither the authors of this software system, nor their employing *
|
|
// * institutes,nor the agencies providing financial support for this *
|
|
// * work make any representation or warranty, express or implied, *
|
|
// * regarding this software system or assume any liability for its *
|
|
// * use. Please see the license in the file LICENSE and URL above *
|
|
// * for the full disclaimer and the limitation of liability. *
|
|
// * *
|
|
// * This code implementation is the result of the scientific and *
|
|
// * technical work of the GEANT4 collaboration. *
|
|
// * By using, copying, modifying or distributing the software (or *
|
|
// * any work based on the software) you agree to acknowledge its *
|
|
// * use in resulting scientific publications, and indicate your *
|
|
// * acceptance of all terms of the Geant4 Software license. *
|
|
// ********************************************************************
|
|
//
|
|
// class G4GeomTools implementation
|
|
//
|
|
// 10.10.2016, E.Tcherniaev: initial version.
|
|
// --------------------------------------------------------------------
|
|
|
|
#include "G4GeomTools.hh"
|
|
|
|
#include "geomdefs.hh"
|
|
#include "G4SystemOfUnits.hh"
|
|
#include "G4GeometryTolerance.hh"
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate area of a triangle in 2D
|
|
|
|
G4double G4GeomTools::TriangleArea(G4double Ax, G4double Ay,
|
|
G4double Bx, G4double By,
|
|
G4double Cx, G4double Cy)
|
|
{
|
|
return ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax))*0.5;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate area of a triangle in 2D
|
|
|
|
G4double G4GeomTools::TriangleArea(const G4TwoVector& A,
|
|
const G4TwoVector& B,
|
|
const G4TwoVector& C)
|
|
{
|
|
G4double Ax = A.x(), Ay = A.y();
|
|
return ((B.x()-Ax)*(C.y()-Ay) - (B.y()-Ay)*(C.x()-Ax))*0.5;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate area of a quadrilateral in 2D
|
|
|
|
G4double G4GeomTools::QuadArea(const G4TwoVector& A,
|
|
const G4TwoVector& B,
|
|
const G4TwoVector& C,
|
|
const G4TwoVector& D)
|
|
{
|
|
return ((C.x()-A.x())*(D.y()-B.y()) - (C.y()-A.y())*(D.x()-B.x()))*0.5;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate area of a polygon in 2D
|
|
|
|
G4double G4GeomTools::PolygonArea(const G4TwoVectorList& p)
|
|
{
|
|
auto n = (G4int)p.size();
|
|
if (n < 3) { return 0.0; // degenerate polygon
|
|
}
|
|
G4double area = p[n-1].x()*p[0].y() - p[0].x()*p[n-1].y();
|
|
for(G4int i=1; i<n; ++i)
|
|
{
|
|
area += p[i-1].x()*p[i].y() - p[i].x()*p[i-1].y();
|
|
}
|
|
return area*0.5;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Point inside 2D triangle
|
|
|
|
G4bool G4GeomTools::PointInTriangle(G4double Ax, G4double Ay,
|
|
G4double Bx, G4double By,
|
|
G4double Cx, G4double Cy,
|
|
G4double Px, G4double Py)
|
|
|
|
{
|
|
if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
|
|
{
|
|
if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.) { return false; }
|
|
if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.) { return false; }
|
|
if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.) { return false; }
|
|
}
|
|
else
|
|
{
|
|
if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.) { return false; }
|
|
if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.) { return false; }
|
|
if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.) { return false; }
|
|
}
|
|
return true;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Point inside 2D triangle
|
|
|
|
G4bool G4GeomTools::PointInTriangle(const G4TwoVector& A,
|
|
const G4TwoVector& B,
|
|
const G4TwoVector& C,
|
|
const G4TwoVector& P)
|
|
{
|
|
G4double Ax = A.x(), Ay = A.y();
|
|
G4double Bx = B.x(), By = B.y();
|
|
G4double Cx = C.x(), Cy = C.y();
|
|
G4double Px = P.x(), Py = P.y();
|
|
if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
|
|
{
|
|
if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.) { return false; }
|
|
if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.) { return false; }
|
|
if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.) { return false; }
|
|
}
|
|
else
|
|
{
|
|
if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.) { return false; }
|
|
if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.) { return false; }
|
|
if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.) { return false; }
|
|
}
|
|
return true;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Point inside 2D polygon
|
|
|
|
G4bool G4GeomTools::PointInPolygon(const G4TwoVector& p,
|
|
const G4TwoVectorList& v)
|
|
{
|
|
auto Nv = (G4int)v.size();
|
|
G4bool in = false;
|
|
for (G4int i = 0, k = Nv - 1; i < Nv; k = i++)
|
|
{
|
|
if ((v[i].y() > p.y()) != (v[k].y() > p.y()))
|
|
{
|
|
G4double ctg = (v[k].x()-v[i].x())/(v[k].y()-v[i].y());
|
|
in ^= static_cast<int>(p.x() < (p.y()-v[i].y())*ctg + v[i].x());
|
|
}
|
|
}
|
|
return in;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Detemine whether 2D polygon is convex or not
|
|
|
|
G4bool G4GeomTools::IsConvex(const G4TwoVectorList& polygon)
|
|
{
|
|
static const G4double kCarTolerance =
|
|
G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
|
|
|
|
G4bool gotNegative = false;
|
|
G4bool gotPositive = false;
|
|
auto n = (G4int)polygon.size();
|
|
if (n <= 0) { return false; }
|
|
for (G4int icur=0; icur<n; ++icur)
|
|
{
|
|
G4int iprev = (icur == 0) ? n-1 : icur-1;
|
|
G4int inext = (icur == n-1) ? 0 : icur+1;
|
|
G4TwoVector e1 = polygon[icur] - polygon[iprev];
|
|
G4TwoVector e2 = polygon[inext] - polygon[icur];
|
|
G4double cross = e1.x()*e2.y() - e1.y()*e2.x();
|
|
if (std::abs(cross) < kCarTolerance) { return false; }
|
|
if (cross < 0) { gotNegative = true; }
|
|
if (cross > 0) { gotPositive = true; }
|
|
if (gotNegative && gotPositive) { return false; }
|
|
}
|
|
return true;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Triangulate simple polygon
|
|
|
|
G4bool G4GeomTools::TriangulatePolygon(const G4TwoVectorList& polygon,
|
|
G4TwoVectorList& result)
|
|
{
|
|
result.resize(0);
|
|
std::vector<G4int> triangles;
|
|
G4bool reply = TriangulatePolygon(polygon,triangles);
|
|
|
|
auto n = (G4int)triangles.size();
|
|
for (G4int i=0; i<n; ++i) { result.push_back(polygon[triangles[i]]); }
|
|
return reply;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Triangulation of a simple polygon by "ear clipping"
|
|
|
|
G4bool G4GeomTools::TriangulatePolygon(const G4TwoVectorList& polygon,
|
|
std::vector<G4int>& result)
|
|
{
|
|
result.resize(0);
|
|
|
|
// allocate and initialize list of Vertices in polygon
|
|
//
|
|
auto n = (G4int)polygon.size();
|
|
if (n < 3) { return false; }
|
|
|
|
// we want a counter-clockwise polygon in V
|
|
//
|
|
G4double area = G4GeomTools::PolygonArea(polygon);
|
|
auto V = new G4int[n];
|
|
if (area > 0.) {
|
|
for (G4int i=0; i<n; ++i) { V[i] = i; }
|
|
} else {
|
|
for (G4int i=0; i<n; ++i) { V[i] = (n-1)-i; }
|
|
}
|
|
|
|
// Triangulation: remove nv-2 Vertices, creating 1 triangle every time
|
|
//
|
|
G4int nv = n;
|
|
G4int count = 2*nv; // error detection counter
|
|
for(G4int b=nv-1; nv>2; )
|
|
{
|
|
// ERROR: if we loop, it is probably a non-simple polygon
|
|
if ((count--) <= 0)
|
|
{
|
|
delete [] V;
|
|
if (area < 0.) { std::reverse(result.begin(),result.end()); }
|
|
return false;
|
|
}
|
|
|
|
// three consecutive vertices in current polygon, <a,b,c>
|
|
G4int a = (b < nv) ? b : 0; // previous
|
|
b = (a+1 < nv) ? a+1 : 0; // current
|
|
G4int c = (b+1 < nv) ? b+1 : 0; // next
|
|
|
|
if (CheckSnip(polygon, a,b,c, nv,V))
|
|
{
|
|
// output Triangle
|
|
result.push_back(V[a]);
|
|
result.push_back(V[b]);
|
|
result.push_back(V[c]);
|
|
|
|
// remove vertex b from remaining polygon
|
|
nv--;
|
|
for(G4int i=b; i<nv; ++i) { V[i] = V[i+1]; }
|
|
|
|
count = 2*nv; // resest error detection counter
|
|
}
|
|
}
|
|
delete [] V;
|
|
if (area < 0.) { std::reverse(result.begin(),result.end()); }
|
|
return true;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Helper function for "ear clipping" polygon triangulation.
|
|
// Check for a valid snip
|
|
|
|
G4bool G4GeomTools::CheckSnip(const G4TwoVectorList& contour,
|
|
G4int a, G4int b, G4int c,
|
|
G4int n, const G4int* V)
|
|
{
|
|
static const G4double kCarTolerance =
|
|
G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
|
|
|
|
// check orientation of Triangle
|
|
G4double Ax = contour[V[a]].x(), Ay = contour[V[a]].y();
|
|
G4double Bx = contour[V[b]].x(), By = contour[V[b]].y();
|
|
G4double Cx = contour[V[c]].x(), Cy = contour[V[c]].y();
|
|
if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCarTolerance) { return false; }
|
|
|
|
// check that there is no point inside Triangle
|
|
G4double xmin = std::min(std::min(Ax,Bx),Cx);
|
|
G4double xmax = std::max(std::max(Ax,Bx),Cx);
|
|
G4double ymin = std::min(std::min(Ay,By),Cy);
|
|
G4double ymax = std::max(std::max(Ay,By),Cy);
|
|
for (G4int i=0; i<n; ++i)
|
|
{
|
|
if((i == a) || (i == b) || (i == c)) { continue; }
|
|
G4double Px = contour[V[i]].x();
|
|
if (Px < xmin || Px > xmax) { continue; }
|
|
G4double Py = contour[V[i]].y();
|
|
if (Py < ymin || Py > ymax) { continue; }
|
|
if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) { return false; }
|
|
}
|
|
return true;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Remove collinear and coincident points from 2D polygon
|
|
|
|
void G4GeomTools::RemoveRedundantVertices(G4TwoVectorList& polygon,
|
|
std::vector<G4int>& iout,
|
|
G4double tolerance)
|
|
{
|
|
iout.resize(0);
|
|
// set tolerance squared
|
|
G4double delta = sqr(tolerance);
|
|
// set special value to mark vertices for removal
|
|
G4double removeIt = kInfinity;
|
|
|
|
auto nv = (G4int)polygon.size();
|
|
|
|
// Main loop: check every three consecutive points, if the points
|
|
// are collinear then mark middle point for removal
|
|
//
|
|
G4int icur = 0, iprev = 0, inext = 0, nout = 0;
|
|
for (G4int i=0; i<nv; ++i)
|
|
{
|
|
icur = i; // index of current point
|
|
|
|
for (G4int k=1; k<nv+1; ++k) // set index of previous point
|
|
{
|
|
iprev = icur - k;
|
|
if (iprev < 0) { iprev += nv; }
|
|
if (polygon[iprev].x() != removeIt) { break; }
|
|
}
|
|
|
|
for (G4int k=1; k<nv+1; ++k) // set index of next point
|
|
{
|
|
inext = icur + k;
|
|
if (inext >= nv) { inext -= nv; }
|
|
if (polygon[inext].x() != removeIt) { break; }
|
|
}
|
|
|
|
if (iprev == inext) { break; } // degenerate polygon, stop
|
|
|
|
// Calculate parameters of triangle (iprev->icur->inext),
|
|
// if triangle is too small or too narrow then mark current
|
|
// point for removal
|
|
G4TwoVector e1 = polygon[iprev] - polygon[icur];
|
|
G4TwoVector e2 = polygon[inext] - polygon[icur];
|
|
|
|
// Check length of edges, then check height of the triangle
|
|
G4double leng1 = e1.mag2();
|
|
G4double leng2 = e2.mag2();
|
|
G4double leng3 = (e2-e1).mag2();
|
|
if (leng1 <= delta || leng2 <= delta || leng3 <= delta)
|
|
{
|
|
polygon[icur].setX(removeIt); ++nout;
|
|
}
|
|
else
|
|
{
|
|
G4double lmax = std::max(std::max(leng1,leng2),leng3);
|
|
G4double area = std::abs(e1.x()*e2.y()-e1.y()*e2.x())*0.5;
|
|
if (area/std::sqrt(lmax) <= std::abs(tolerance))
|
|
{
|
|
polygon[icur].setX(removeIt); ++nout;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Remove marked points
|
|
//
|
|
icur = 0;
|
|
if (nv - nout < 3) // degenerate polygon, remove all points
|
|
{
|
|
for (G4int i=0; i<nv; ++i) { iout.push_back(i); }
|
|
polygon.resize(0);
|
|
nv = 0;
|
|
}
|
|
for (G4int i=0; i<nv; ++i) // move points, if required
|
|
{
|
|
if (polygon[i].x() != removeIt) {
|
|
polygon[icur++] = polygon[i];
|
|
} else {
|
|
iout.push_back(i); }
|
|
}
|
|
if (icur < nv) { polygon.resize(icur); }
|
|
return;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Find bounding rectangle of a disk sector
|
|
|
|
G4bool G4GeomTools::DiskExtent(G4double rmin, G4double rmax,
|
|
G4double startPhi, G4double delPhi,
|
|
G4TwoVector& pmin, G4TwoVector& pmax)
|
|
{
|
|
static const G4double kCarTolerance =
|
|
G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
|
|
|
|
// check parameters
|
|
//
|
|
pmin.set(0,0);
|
|
pmax.set(0,0);
|
|
if (rmin < 0) { return false; }
|
|
if (rmax <= rmin + kCarTolerance) { return false; }
|
|
if (delPhi <= 0 + kCarTolerance) { return false; }
|
|
|
|
// calculate extent
|
|
//
|
|
pmin.set(-rmax,-rmax);
|
|
pmax.set( rmax, rmax);
|
|
if (delPhi >= CLHEP::twopi) { return true; }
|
|
|
|
DiskExtent(rmin,rmax,
|
|
std::sin(startPhi),std::cos(startPhi),
|
|
std::sin(startPhi+delPhi),std::cos(startPhi+delPhi),
|
|
pmin,pmax);
|
|
return true;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Find bounding rectangle of a disk sector, fast version.
|
|
// No check of parameters !!!
|
|
|
|
void G4GeomTools::DiskExtent(G4double rmin, G4double rmax,
|
|
G4double sinStart, G4double cosStart,
|
|
G4double sinEnd, G4double cosEnd,
|
|
G4TwoVector& pmin, G4TwoVector& pmax)
|
|
{
|
|
static const G4double kCarTolerance =
|
|
G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
|
|
|
|
// check if 360 degrees
|
|
//
|
|
pmin.set(-rmax,-rmax);
|
|
pmax.set( rmax, rmax);
|
|
|
|
if (std::abs(sinEnd-sinStart) < kCarTolerance &&
|
|
std::abs(cosEnd-cosStart) < kCarTolerance) { return; }
|
|
|
|
// get start and end quadrants
|
|
//
|
|
// 1 | 0
|
|
// ---+---
|
|
// 3 | 2
|
|
//
|
|
G4int icase = (cosEnd < 0) ? 1 : 0;
|
|
if (sinEnd < 0) { icase += 2; }
|
|
if (cosStart < 0) { icase += 4; }
|
|
if (sinStart < 0) { icase += 8; }
|
|
|
|
switch (icase)
|
|
{
|
|
// start quadrant 0
|
|
case 0: // start->end : 0->0
|
|
if (sinEnd < sinStart) { break; }
|
|
pmin.set(rmin*cosEnd,rmin*sinStart);
|
|
pmax.set(rmax*cosStart,rmax*sinEnd );
|
|
break;
|
|
case 1: // start->end : 0->1
|
|
pmin.set(rmax*cosEnd,std::min(rmin*sinStart,rmin*sinEnd));
|
|
pmax.set(rmax*cosStart,rmax );
|
|
break;
|
|
case 2: // start->end : 0->2
|
|
pmin.set(-rmax,-rmax);
|
|
pmax.set(std::max(rmax*cosStart,rmax*cosEnd),rmax);
|
|
break;
|
|
case 3: // start->end : 0->3
|
|
pmin.set(-rmax,rmax*sinEnd);
|
|
pmax.set(rmax*cosStart,rmax);
|
|
break;
|
|
// start quadrant 1
|
|
case 4: // start->end : 1->0
|
|
pmin.set(-rmax,-rmax);
|
|
pmax.set(rmax,std::max(rmax*sinStart,rmax*sinEnd));
|
|
break;
|
|
case 5: // start->end : 1->1
|
|
if (sinEnd > sinStart) { break; }
|
|
pmin.set(rmax*cosEnd,rmin*sinEnd );
|
|
pmax.set(rmin*cosStart,rmax*sinStart);
|
|
break;
|
|
case 6: // start->end : 1->2
|
|
pmin.set(-rmax,-rmax);
|
|
pmax.set(rmax*cosEnd,rmax*sinStart);
|
|
break;
|
|
case 7: // start->end : 1->3
|
|
pmin.set(-rmax,rmax*sinEnd);
|
|
pmax.set(std::max(rmin*cosStart,rmin*cosEnd),rmax*sinStart);
|
|
break;
|
|
// start quadrant 2
|
|
case 8: // start->end : 2->0
|
|
pmin.set(std::min(rmin*cosStart,rmin*cosEnd),rmax*sinStart);
|
|
pmax.set(rmax,rmax*sinEnd);
|
|
break;
|
|
case 9: // start->end : 2->1
|
|
pmin.set(rmax*cosEnd,rmax*sinStart);
|
|
pmax.set(rmax,rmax);
|
|
break;
|
|
case 10: // start->end : 2->2
|
|
if (sinEnd < sinStart) { break; }
|
|
pmin.set(rmin*cosStart,rmax*sinStart);
|
|
pmax.set(rmax*cosEnd,rmin*sinEnd );
|
|
break;
|
|
case 11: // start->end : 2->3
|
|
pmin.set(-rmax,std::min(rmax*sinStart,rmax*sinEnd));
|
|
pmax.set(rmax,rmax);
|
|
break;
|
|
// start quadrant 3
|
|
case 12: // start->end : 3->0
|
|
pmin.set(rmax*cosStart,-rmax);
|
|
pmax.set(rmax,rmax*sinEnd);
|
|
break;
|
|
case 13: // start->end : 3->1
|
|
pmin.set(std::min(rmax*cosStart,rmax*cosEnd),-rmax);
|
|
pmax.set(rmax,rmax);
|
|
break;
|
|
case 14: // start->end : 3->2
|
|
pmin.set(rmax*cosStart,-rmax);
|
|
pmax.set(rmax*cosEnd,std::max(rmin*sinStart,rmin*sinEnd));
|
|
break;
|
|
case 15: // start->end : 3->3
|
|
if (sinEnd > sinStart) { break; }
|
|
pmin.set(rmax*cosStart,rmax*sinEnd);
|
|
pmax.set(rmin*cosEnd,rmin*sinStart);
|
|
break;
|
|
}
|
|
return;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Compute the circumference (perimeter) of an ellipse
|
|
|
|
G4double G4GeomTools::EllipsePerimeter(G4double pA, G4double pB)
|
|
{
|
|
G4double x = std::abs(pA);
|
|
G4double y = std::abs(pB);
|
|
G4double a = std::max(x,y);
|
|
G4double b = std::min(x,y);
|
|
G4double e = std::sqrt((1. - b/a)*(1. + b/a));
|
|
return 4. * a * comp_ellint_2(e);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Compute the lateral surface area of an elliptic cone
|
|
|
|
G4double G4GeomTools::EllipticConeLateralArea(G4double pA,
|
|
G4double pB,
|
|
G4double pH)
|
|
{
|
|
G4double x = std::abs(pA);
|
|
G4double y = std::abs(pB);
|
|
G4double h = std::abs(pH);
|
|
G4double a = std::max(x,y);
|
|
G4double b = std::min(x,y);
|
|
G4double e = std::sqrt((1. - b/a)*(1. + b/a)) / std::hypot(1.,b/h);
|
|
return 2. * a * std::hypot(b,h) * comp_ellint_2(e);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Compute Elliptical Integral of the Second Kind
|
|
//
|
|
// The algorithm is based upon Carlson B.C., "Computation of real
|
|
// or complex elliptic integrals", Numerical Algorithms,
|
|
// Volume 10, Issue 1, 1995 (see equations 2.36 - 2.39)
|
|
//
|
|
// The code was adopted from C code at:
|
|
// http://paulbourke.net/geometry/ellipsecirc/
|
|
|
|
G4double G4GeomTools::comp_ellint_2(G4double e)
|
|
{
|
|
const G4double eps = 1. / 134217728.; // 1 / 2^27
|
|
|
|
G4double a = 1.;
|
|
G4double b = std::sqrt((1. - e)*(1. + e));
|
|
if (b == 1.) { return CLHEP::halfpi; }
|
|
if (b == 0.) { return 1.; }
|
|
|
|
G4double x = 1.;
|
|
G4double y = b;
|
|
G4double S = 0.;
|
|
G4double M = 1.;
|
|
while (x - y > eps*y)
|
|
{
|
|
G4double tmp = (x + y) * 0.5;
|
|
y = std::sqrt(x*y);
|
|
x = tmp;
|
|
M += M;
|
|
S += M * (x - y)*(x - y);
|
|
}
|
|
return 0.5 * CLHEP::halfpi * ((a + b)*(a + b) - S) / (x + y);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calcuate area of a triangle in 3D
|
|
|
|
G4ThreeVector G4GeomTools::TriangleAreaNormal(const G4ThreeVector& A,
|
|
const G4ThreeVector& B,
|
|
const G4ThreeVector& C)
|
|
{
|
|
return ((B-A).cross(C-A))*0.5;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calcuate area of a quadrilateral in 3D
|
|
|
|
G4ThreeVector G4GeomTools::QuadAreaNormal(const G4ThreeVector& A,
|
|
const G4ThreeVector& B,
|
|
const G4ThreeVector& C,
|
|
const G4ThreeVector& D)
|
|
{
|
|
return ((C-A).cross(D-B))*0.5;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate area of a polygon in 3D
|
|
|
|
G4ThreeVector G4GeomTools::PolygonAreaNormal(const G4ThreeVectorList& p)
|
|
{
|
|
auto n = (G4int)p.size();
|
|
if (n < 3) { return {0,0,0}; } // degerate polygon
|
|
G4ThreeVector normal = p[n-1].cross(p[0]);
|
|
for(G4int i=1; i<n; ++i)
|
|
{
|
|
normal += p[i-1].cross(p[i]);
|
|
}
|
|
return normal*0.5;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate distance between point P and line segment AB in 3D
|
|
|
|
G4double G4GeomTools::DistancePointSegment(const G4ThreeVector& P,
|
|
const G4ThreeVector& A,
|
|
const G4ThreeVector& B)
|
|
{
|
|
G4ThreeVector AP = P - A;
|
|
G4ThreeVector AB = B - A;
|
|
|
|
G4double u = AP.dot(AB);
|
|
if (u <= 0) { return AP.mag(); } // closest point is A
|
|
|
|
G4double len2 = AB.mag2();
|
|
if (u >= len2) { return (B-P).mag(); } // closest point is B
|
|
|
|
return ((u/len2)*AB - AP).mag(); // distance to line
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Find closest point on line segment in 3D
|
|
|
|
G4ThreeVector
|
|
G4GeomTools::ClosestPointOnSegment(const G4ThreeVector& P,
|
|
const G4ThreeVector& A,
|
|
const G4ThreeVector& B)
|
|
{
|
|
G4ThreeVector AP = P - A;
|
|
G4ThreeVector AB = B - A;
|
|
|
|
G4double u = AP.dot(AB);
|
|
if (u <= 0) { return A; } // closest point is A
|
|
|
|
G4double len2 = AB.mag2();
|
|
if (u >= len2) { return B; } // closest point is B
|
|
|
|
G4double t = u/len2;
|
|
return A + t*AB; // closest point on segment
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Find closest point on triangle in 3D.
|
|
//
|
|
// The implementation is based on the algorithm published in
|
|
// "Geometric Tools for Computer Graphics", Philip J Scheider and
|
|
// David H Eberly, Elsevier Science (USA), 2003.
|
|
//
|
|
// The algorithm is also available at:
|
|
// http://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
|
|
|
|
G4ThreeVector
|
|
G4GeomTools::ClosestPointOnTriangle(const G4ThreeVector& P,
|
|
const G4ThreeVector& A,
|
|
const G4ThreeVector& B,
|
|
const G4ThreeVector& C)
|
|
{
|
|
G4ThreeVector diff = A - P;
|
|
G4ThreeVector edge0 = B - A;
|
|
G4ThreeVector edge1 = C - A;
|
|
|
|
G4double a = edge0.mag2();
|
|
G4double b = edge0.dot(edge1);
|
|
G4double c = edge1.mag2();
|
|
G4double d = diff.dot(edge0);
|
|
G4double e = diff.dot(edge1);
|
|
|
|
G4double det = a*c - b*b;
|
|
G4double t0 = b*e - c*d;
|
|
G4double t1 = b*d - a*e;
|
|
|
|
/*
|
|
^ t1
|
|
\ 2 |
|
|
\ |
|
|
\ | regions
|
|
\|
|
|
C
|
|
|\
|
|
3 | \ 1
|
|
| \
|
|
| 0 \
|
|
| \
|
|
---- A --- B ----> t0
|
|
| \
|
|
4 | 5 \ 6
|
|
| \
|
|
*/
|
|
|
|
G4int region = -1;
|
|
if (t0+t1 <= det) {
|
|
region = (t0 < 0) ? ((t1 < 0) ? 4 : 3) : ((t1 < 0) ? 5 : 0);
|
|
} else {
|
|
region = (t0 < 0) ? 2 : ((t1 < 0) ? 6 : 1);
|
|
}
|
|
|
|
switch (region)
|
|
{
|
|
case 0: // interior of triangle
|
|
{
|
|
G4double invDet = 1./det;
|
|
return A + (t0*invDet)*edge0 + (t1*invDet)*edge1;
|
|
}
|
|
case 1: // edge BC
|
|
{
|
|
G4double numer = c + e - b - d;
|
|
if (numer <= 0) { return C; }
|
|
G4double denom = a - 2*b + c;
|
|
return (numer >= denom) ? B : C + (numer/denom)*(edge0-edge1);
|
|
}
|
|
case 2: // edge AC or BC
|
|
{
|
|
G4double tmp0 = b + d;
|
|
G4double tmp1 = c + e;
|
|
if (tmp1 > tmp0)
|
|
{
|
|
G4double numer = tmp1 - tmp0;
|
|
G4double denom = a - 2*b + c;
|
|
return (numer >= denom) ? B : C + (numer/denom)*(edge0-edge1);
|
|
}
|
|
// same: (e >= 0) ? A : ((-e >= c) ? C : A + (-e/c)*edge1)
|
|
return (tmp1 <= 0) ? C : (( e >= 0) ? A : A + (-e/c)*edge1);
|
|
}
|
|
case 3: // edge AC
|
|
return (e >= 0) ? A : ((-e >= c) ? C : A + (-e/c)*edge1);
|
|
|
|
case 4: // edge AB or AC
|
|
if (d < 0) { return (-d >= a) ? B : A + (-d/a)*edge0; }
|
|
return (e >= 0) ? A : ((-e >= c) ? C : A + (-e/c)*edge1);
|
|
|
|
case 5: // edge AB
|
|
return (d >= 0) ? A : ((-d >= a) ? B : A + (-d/a)*edge0);
|
|
|
|
case 6: // edge AB or BC
|
|
{
|
|
G4double tmp0 = b + e;
|
|
G4double tmp1 = a + d;
|
|
if (tmp1 > tmp0)
|
|
{
|
|
G4double numer = tmp1 - tmp0;
|
|
G4double denom = a - 2*b + c;
|
|
return (numer >= denom) ? C : B + (numer/denom)*(edge1-edge0);
|
|
}
|
|
// same: (d >= 0) ? A : ((-d >= a) ? B : A + (-d/a)*edge0)
|
|
return (tmp1 <= 0) ? B : (( d >= 0) ? A : A + (-d/a)*edge0);
|
|
}
|
|
default: // impossible case
|
|
return {kInfinity,kInfinity,kInfinity};
|
|
}
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate bounding box of a spherical sector
|
|
|
|
G4bool
|
|
G4GeomTools::SphereExtent(G4double rmin, G4double rmax,
|
|
G4double startTheta, G4double delTheta,
|
|
G4double startPhi, G4double delPhi,
|
|
G4ThreeVector& pmin, G4ThreeVector& pmax)
|
|
{
|
|
static const G4double kCarTolerance =
|
|
G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
|
|
|
|
// check parameters
|
|
//
|
|
pmin.set(0,0,0);
|
|
pmax.set(0,0,0);
|
|
if (rmin < 0) { return false; }
|
|
if (rmax <= rmin + kCarTolerance) { return false; }
|
|
if (delTheta <= 0 + kCarTolerance) { return false; }
|
|
if (delPhi <= 0 + kCarTolerance) { return false; }
|
|
|
|
G4double stheta = startTheta;
|
|
G4double dtheta = delTheta;
|
|
if (stheta < 0 && stheta > CLHEP::pi) { return false; }
|
|
if (stheta + dtheta > CLHEP::pi) { dtheta = CLHEP::pi - stheta; }
|
|
if (dtheta <= 0 + kCarTolerance) { return false; }
|
|
|
|
// calculate extent
|
|
//
|
|
pmin.set(-rmax,-rmax,-rmax);
|
|
pmax.set( rmax, rmax, rmax);
|
|
if (dtheta >= CLHEP::pi && delPhi >= CLHEP::twopi) { return true; }
|
|
|
|
G4double etheta = stheta + dtheta;
|
|
G4double sinStart = std::sin(stheta);
|
|
G4double cosStart = std::cos(stheta);
|
|
G4double sinEnd = std::sin(etheta);
|
|
G4double cosEnd = std::cos(etheta);
|
|
|
|
G4double rhomin = rmin*std::min(sinStart,sinEnd);
|
|
G4double rhomax = rmax;
|
|
if (stheta > CLHEP::halfpi) { rhomax = rmax*sinStart; }
|
|
if (etheta < CLHEP::halfpi) { rhomax = rmax*sinEnd; }
|
|
|
|
G4TwoVector xymin,xymax;
|
|
DiskExtent(rhomin,rhomax,
|
|
std::sin(startPhi),std::cos(startPhi),
|
|
std::sin(startPhi+delPhi),std::cos(startPhi+delPhi),
|
|
xymin,xymax);
|
|
|
|
G4double zmin = std::min(rmin*cosEnd,rmax*cosEnd);
|
|
G4double zmax = std::max(rmin*cosStart,rmax*cosStart);
|
|
pmin.set(xymin.x(),xymin.y(),zmin);
|
|
pmax.set(xymax.x(),xymax.y(),zmax);
|
|
return true;
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate hyperbolic surface stereo
|
|
|
|
G4double
|
|
G4GeomTools::HypeStereo(G4double r0, G4double r, G4double h)
|
|
{
|
|
static const G4double kCarTolerance =
|
|
G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
|
|
if (std::abs(r - r0) < kCarTolerance) { return 0.; }
|
|
return std::atan(std::sqrt((r - r0)*(r + r0))/std::abs(h));
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Find XY-coordinates of the corners of the bounding generic trap
|
|
// for the specified twisted tube
|
|
|
|
void
|
|
G4GeomTools::TwistedTubeBoundingTrap(G4double twistAng,
|
|
G4double endInnerRad,
|
|
G4double endOuterRad,
|
|
G4double dPhi,
|
|
G4TwoVectorList& vertices)
|
|
{
|
|
vertices.resize(8);
|
|
G4double rmin = std::abs(endInnerRad);
|
|
G4double rmax = std::abs(endOuterRad);
|
|
|
|
// Set untwisted vertices
|
|
G4double phi = dPhi/2.;
|
|
G4double sinphi = std::sin(phi);
|
|
G4double cosphi = std::cos(phi);
|
|
G4double tanphi = std::tan(phi);
|
|
vertices[0].set(rmin*cosphi, rmin*sinphi);
|
|
vertices[1].set(rmax, rmax*tanphi);
|
|
vertices[2].set(rmax,-rmax*tanphi);
|
|
vertices[3].set(rmin*cosphi,-rmin*sinphi);
|
|
vertices[4] = vertices[0];
|
|
vertices[5] = vertices[1];
|
|
vertices[6] = vertices[2];
|
|
vertices[7] = vertices[3];
|
|
|
|
// Twist vertices
|
|
G4double ang = twistAng/2.;
|
|
for(auto i = 0; i < 4; ++i)
|
|
{
|
|
vertices[i].rotate(-ang); // vertices at -halfz
|
|
vertices[i + 4].rotate(ang); // vertices at +halfz
|
|
}
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Calculate surface area of hyperboloid between zmin and zmax
|
|
|
|
G4double
|
|
G4GeomTools::HyperboloidSurfaceArea(G4double dphi, G4double r0, G4double tanstereo,
|
|
G4double zmin, G4double zmax)
|
|
{
|
|
static const G4double kCarTolerance =
|
|
G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
|
|
|
|
G4double a = std::abs(r0); // radius at z = 0
|
|
G4double t = std::abs(tanstereo); // tan(stereo)
|
|
G4double phi = std::abs(dphi); // delta phi
|
|
|
|
// Check spesial cases: cylindrical and conical surfaces
|
|
if (t < kCarTolerance) { return a*std::abs(zmax - zmin)*phi; } // cylinder
|
|
G4double rmin = std::hypot(t*zmin, a); // radius at zmin
|
|
G4double rmax = std::hypot(t*zmax, a); // radius at zmax
|
|
if (a < kCarTolerance) // cone
|
|
{
|
|
G4double smin = rmin*std::hypot(rmin, zmin);
|
|
G4double smax = rmax*std::hypot(rmax, zmax);
|
|
return (zmin*zmax < 0.) ? (smin + smax)*phi/2. : std::abs(smax - smin)*phi/2.;
|
|
}
|
|
// Find surface area
|
|
G4double tt = t*t;
|
|
G4double aa = a*a;
|
|
G4double cc = aa/tt;
|
|
G4double k = std::sqrt(aa + cc)/cc;
|
|
|
|
G4double hmin = std::abs(zmin);
|
|
G4double smin = a*(hmin*std::hypot(1., k*hmin) + std::asinh(k*hmin)/k);
|
|
if (zmax == -zmin) { return smin*phi; }
|
|
|
|
G4double hmax = std::abs(zmax);
|
|
G4double smax = a*(hmax*std::hypot(1., k*hmax) + std::asinh(k*hmax)/k);
|
|
return (zmin*zmax < 0.) ? (smin + smax)*phi/2. :std::abs(smax - smin)*phi/2.;
|
|
}
|