source: trunk/source/graphics_reps/src/HepPolyhedron.cc@ 1206

Last change on this file since 1206 was 1140, checked in by garnier, 16 years ago

update to CVS
File size: 84.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: HepPolyhedron.cc,v 1.34 2009/10/28 13:36:32 allison Exp $
28// GEANT4 tag $Name: $
29//
30//
31//
32// G4 Polyhedron library
33//
34// History:
35// 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
36//
37// 30.09.96 E.Chernyaev
38// - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
39// - added GetNextUnitNormal, GetNextEdgeIndeces, GetNextEdge
40//
41// 15.12.96 E.Chernyaev
42// - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
43// - rewritten G4PolyhedronCons;
44// - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
45//
46// 01.06.97 E.Chernyaev
47// - modified RotateAroundZ, added SetSideFacets
48//
49// 19.03.00 E.Chernyaev
50// - implemented boolean operations (add, subtract, intersect) on polyhedra;
51//
52// 25.05.01 E.Chernyaev
53// - added GetSurfaceArea() and GetVolume();
54//
55// 05.11.02 E.Chernyaev
56// - added createTwistedTrap() and createPolyhedron();
57//
58// 20.06.05 G.Cosmo
59// - added HepPolyhedronEllipsoid;
60//
61// 18.07.07 T.Nikitin
62// - added HepParaboloid;
63
64#include "HepPolyhedron.h"
65#include <CLHEP/Units/PhysicalConstants.h>
66#include <CLHEP/Geometry/Vector3D.h>
67
68#include <cstdlib> // Required on some compilers for std::abs(int) ...
69#include <cmath>
70
71using namespace HepGeom;
72using CLHEP::perMillion;
73using CLHEP::deg;
74using CLHEP::pi;
75using CLHEP::twopi;
76
77/***********************************************************************
78 * *
79 * Name: HepPolyhedron operator << Date: 09.05.96 *
80 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
81 * *
82 * Function: Print contents of G4 polyhedron *
83 * *
84 ***********************************************************************/
85std::ostream & operator<<(std::ostream & ostr, const G4Facet & facet) {
86 for (int k=0; k<4; k++) {
87 ostr << " " << facet.edge[k].v << "/" << facet.edge[k].f;
88 }
89 return ostr;
90}
91
92std::ostream & operator<<(std::ostream & ostr, const HepPolyhedron & ph) {
93 ostr << std::endl;
94 ostr << "Nverteces=" << ph.nvert << ", Nfacets=" << ph.nface << std::endl;
95 int i;
96 for (i=1; i<=ph.nvert; i++) {
97 ostr << "xyz(" << i << ")="
98 << ph.pV[i].x() << ' ' << ph.pV[i].y() << ' ' << ph.pV[i].z()
99 << std::endl;
100 }
101 for (i=1; i<=ph.nface; i++) {
102 ostr << "face(" << i << ")=" << ph.pF[i] << std::endl;
103 }
104 return ostr;
105}
106
107HepPolyhedron::HepPolyhedron(const HepPolyhedron &from)
108/***********************************************************************
109 * *
110 * Name: HepPolyhedron copy constructor Date: 23.07.96 *
111 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
112 * *
113 ***********************************************************************/
114: nvert(0), nface(0), pV(0), pF(0)
115{
116 AllocateMemory(from.nvert, from.nface);
117 for (int i=1; i<=nvert; i++) pV[i] = from.pV[i];
118 for (int k=1; k<=nface; k++) pF[k] = from.pF[k];
119}
120
121HepPolyhedron & HepPolyhedron::operator=(const HepPolyhedron &from)
122/***********************************************************************
123 * *
124 * Name: HepPolyhedron operator = Date: 23.07.96 *
125 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
126 * *
127 * Function: Copy contents of one polyhedron to another *
128 * *
129 ***********************************************************************/
130{
131 if (this != &from) {
132 AllocateMemory(from.nvert, from.nface);
133 for (int i=1; i<=nvert; i++) pV[i] = from.pV[i];
134 for (int k=1; k<=nface; k++) pF[k] = from.pF[k];
135 }
136 return *this;
137}
138
139int
140HepPolyhedron::FindNeighbour(int iFace, int iNode, int iOrder) const
141/***********************************************************************
142 * *
143 * Name: HepPolyhedron::FindNeighbour Date: 22.11.99 *
144 * Author: E.Chernyaev Revised: *
145 * *
146 * Function: Find neighbouring face *
147 * *
148 ***********************************************************************/
149{
150 int i;
151 for (i=0; i<4; i++) {
152 if (iNode == std::abs(pF[iFace].edge[i].v)) break;
153 }
154 if (i == 4) {
155 std::cerr
156 << "HepPolyhedron::FindNeighbour: face " << iFace
157 << " has no node " << iNode
158 << std::endl;
159 return 0;
160 }
161 if (iOrder < 0) {
162 if ( --i < 0) i = 3;
163 if (pF[iFace].edge[i].v == 0) i = 2;
164 }
165 return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
166}
167
168Normal3D<double> HepPolyhedron::FindNodeNormal(int iFace, int iNode) const
169/***********************************************************************
170 * *
171 * Name: HepPolyhedron::FindNodeNormal Date: 22.11.99 *
172 * Author: E.Chernyaev Revised: *
173 * *
174 * Function: Find normal at given node *
175 * *
176 ***********************************************************************/
177{
178 Normal3D<double> normal = GetUnitNormal(iFace);
179 int k = iFace, iOrder = 1, n = 1;
180
181 for(;;) {
182 k = FindNeighbour(k, iNode, iOrder);
183 if (k == iFace) break;
184 if (k > 0) {
185 n++;
186 normal += GetUnitNormal(k);
187 }else{
188 if (iOrder < 0) break;
189 k = iFace;
190 iOrder = -iOrder;
191 }
192 }
193 return normal.unit();
194}
195
196int HepPolyhedron::GetNumberOfRotationSteps()
197/***********************************************************************
198 * *
199 * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
200 * Author: J.Allison (Manchester University) Revised: *
201 * *
202 * Function: Get number of steps for whole circle *
203 * *
204 ***********************************************************************/
205{
206 return fNumberOfRotationSteps;
207}
208
209void HepPolyhedron::SetNumberOfRotationSteps(int n)
210/***********************************************************************
211 * *
212 * Name: HepPolyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
213 * Author: J.Allison (Manchester University) Revised: *
214 * *
215 * Function: Set number of steps for whole circle *
216 * *
217 ***********************************************************************/
218{
219 const int nMin = 3;
220 if (n < nMin) {
221 std::cerr
222 << "HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
223 << "number of steps per circle < " << nMin << "; forced to " << nMin
224 << std::endl;
225 fNumberOfRotationSteps = nMin;
226 }else{
227 fNumberOfRotationSteps = n;
228 }
229}
230
231void HepPolyhedron::ResetNumberOfRotationSteps()
232/***********************************************************************
233 * *
234 * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
235 * Author: J.Allison (Manchester University) Revised: *
236 * *
237 * Function: Reset number of steps for whole circle to default value *
238 * *
239 ***********************************************************************/
240{
241 fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
242}
243
244void HepPolyhedron::AllocateMemory(int Nvert, int Nface)
245/***********************************************************************
246 * *
247 * Name: HepPolyhedron::AllocateMemory Date: 19.06.96 *
248 * Author: E.Chernyaev (IHEP/Protvino) Revised: 05.11.02 *
249 * *
250 * Function: Allocate memory for GEANT4 polyhedron *
251 * *
252 * Input: Nvert - number of nodes *
253 * Nface - number of faces *
254 * *
255 ***********************************************************************/
256{
257 if (nvert == Nvert && nface == Nface) return;
258 if (pV != 0) delete [] pV;
259 if (pF != 0) delete [] pF;
260 if (Nvert > 0 && Nface > 0) {
261 nvert = Nvert;
262 nface = Nface;
263 pV = new Point3D<double>[nvert+1];
264 pF = new G4Facet[nface+1];
265 }else{
266 nvert = 0; nface = 0; pV = 0; pF = 0;
267 }
268}
269
270void HepPolyhedron::CreatePrism()
271/***********************************************************************
272 * *
273 * Name: HepPolyhedron::CreatePrism Date: 15.07.96 *
274 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
275 * *
276 * Function: Set facets for a prism *
277 * *
278 ***********************************************************************/
279{
280 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
281
282 pF[1] = G4Facet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
283 pF[2] = G4Facet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
284 pF[3] = G4Facet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
285 pF[4] = G4Facet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
286 pF[5] = G4Facet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
287 pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
288}
289
290void HepPolyhedron::RotateEdge(int k1, int k2, double r1, double r2,
291 int v1, int v2, int vEdge,
292 bool ifWholeCircle, int ns, int &kface)
293/***********************************************************************
294 * *
295 * Name: HepPolyhedron::RotateEdge Date: 05.12.96 *
296 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
297 * *
298 * Function: Create set of facets by rotation of an edge around Z-axis *
299 * *
300 * Input: k1, k2 - end vertices of the edge *
301 * r1, r2 - radiuses of the end vertices *
302 * v1, v2 - visibility of edges produced by rotation of the end *
303 * vertices *
304 * vEdge - visibility of the edge *
305 * ifWholeCircle - is true in case of whole circle rotation *
306 * ns - number of discrete steps *
307 * r[] - r-coordinates *
308 * kface - current free cell in the pF array *
309 * *
310 ***********************************************************************/
311{
312 if (r1 == 0. && r2 == 0) return;
313
314 int i;
315 int i1 = k1;
316 int i2 = k2;
317 int ii1 = ifWholeCircle ? i1 : i1+ns;
318 int ii2 = ifWholeCircle ? i2 : i2+ns;
319 int vv = ifWholeCircle ? vEdge : 1;
320
321 if (ns == 1) {
322 if (r1 == 0.) {
323 pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0);
324 }else if (r2 == 0.) {
325 pF[kface++] = G4Facet(i1,0, i2,0, v1*(i1+1),0);
326 }else{
327 pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
328 }
329 }else{
330 if (r1 == 0.) {
331 pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
332 for (i2++,i=1; i<ns-1; i2++,i++) {
333 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
334 }
335 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
336 }else if (r2 == 0.) {
337 pF[kface++] = G4Facet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
338 for (i1++,i=1; i<ns-1; i1++,i++) {
339 pF[kface++] = G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
340 }
341 pF[kface++] = G4Facet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
342 }else{
343 pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
344 for (i1++,i2++,i=1; i<ns-1; i1++,i2++,i++) {
345 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
346 }
347 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
348 }
349 }
350}
351
352void HepPolyhedron::SetSideFacets(int ii[4], int vv[4],
353 int *kk, double *r,
354 double dphi, int ns, int &kface)
355/***********************************************************************
356 * *
357 * Name: HepPolyhedron::SetSideFacets Date: 20.05.97 *
358 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
359 * *
360 * Function: Set side facets for the case of incomplete rotation *
361 * *
362 * Input: ii[4] - indeces of original verteces *
363 * vv[4] - visibility of edges *
364 * kk[] - indeces of nodes *
365 * r[] - radiuses *
366 * dphi - delta phi *
367 * ns - number of discrete steps *
368 * kface - current free cell in the pF array *
369 * *
370 ***********************************************************************/
371{
372 int k1, k2, k3, k4;
373
374 if (std::abs((double)(dphi-pi)) < perMillion) { // half a circle
375 for (int i=0; i<4; i++) {
376 k1 = ii[i];
377 k2 = (i == 3) ? ii[0] : ii[i+1];
378 if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
379 }
380 }
381
382 if (ii[1] == ii[2]) {
383 k1 = kk[ii[0]];
384 k2 = kk[ii[2]];
385 k3 = kk[ii[3]];
386 pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
387 if (r[ii[0]] != 0.) k1 += ns;
388 if (r[ii[2]] != 0.) k2 += ns;
389 if (r[ii[3]] != 0.) k3 += ns;
390 pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
391 }else if (kk[ii[0]] == kk[ii[1]]) {
392 k1 = kk[ii[0]];
393 k2 = kk[ii[2]];
394 k3 = kk[ii[3]];
395 pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
396 if (r[ii[0]] != 0.) k1 += ns;
397 if (r[ii[2]] != 0.) k2 += ns;
398 if (r[ii[3]] != 0.) k3 += ns;
399 pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
400 }else if (kk[ii[2]] == kk[ii[3]]) {
401 k1 = kk[ii[0]];
402 k2 = kk[ii[1]];
403 k3 = kk[ii[2]];
404 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
405 if (r[ii[0]] != 0.) k1 += ns;
406 if (r[ii[1]] != 0.) k2 += ns;
407 if (r[ii[2]] != 0.) k3 += ns;
408 pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
409 }else{
410 k1 = kk[ii[0]];
411 k2 = kk[ii[1]];
412 k3 = kk[ii[2]];
413 k4 = kk[ii[3]];
414 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
415 if (r[ii[0]] != 0.) k1 += ns;
416 if (r[ii[1]] != 0.) k2 += ns;
417 if (r[ii[2]] != 0.) k3 += ns;
418 if (r[ii[3]] != 0.) k4 += ns;
419 pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
420 }
421}
422
423void HepPolyhedron::RotateAroundZ(int nstep, double phi, double dphi,
424 int np1, int np2,
425 const double *z, double *r,
426 int nodeVis, int edgeVis)
427/***********************************************************************
428 * *
429 * Name: HepPolyhedron::RotateAroundZ Date: 27.11.96 *
430 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
431 * *
432 * Function: Create HepPolyhedron for a solid produced by rotation of *
433 * two polylines around Z-axis *
434 * *
435 * Input: nstep - number of discrete steps, if 0 then default *
436 * phi - starting phi angle *
437 * dphi - delta phi *
438 * np1 - number of points in external polyline *
439 * (must be negative in case of closed polyline) *
440 * np2 - number of points in internal polyline (may be 1) *
441 * z[] - z-coordinates (+z >>> -z for both polylines) *
442 * r[] - r-coordinates *
443 * nodeVis - how to Draw edges joing consecutive positions of *
444 * node during rotation *
445 * edgeVis - how to Draw edges *
446 * *
447 ***********************************************************************/
448{
449 static double wholeCircle = twopi;
450
451 // S E T R O T A T I O N P A R A M E T E R S
452
453 bool ifWholeCircle = (std::abs(dphi-wholeCircle) < perMillion) ? true : false;
454 double delPhi = ifWholeCircle ? wholeCircle : dphi;
455 int nSphi = (nstep > 0) ?
456 nstep : int(delPhi*GetNumberOfRotationSteps()/wholeCircle+.5);
457 if (nSphi == 0) nSphi = 1;
458 int nVphi = ifWholeCircle ? nSphi : nSphi+1;
459 bool ifClosed = np1 > 0 ? false : true;
460
461 // C O U N T V E R T E C E S
462
463 int absNp1 = std::abs(np1);
464 int absNp2 = std::abs(np2);
465 int i1beg = 0;
466 int i1end = absNp1-1;
467 int i2beg = absNp1;
468 int i2end = absNp1+absNp2-1;
469 int i, j, k;
470
471 for(i=i1beg; i<=i2end; i++) {
472 if (std::abs(r[i]) < perMillion) r[i] = 0.;
473 }
474
475 j = 0; // external nodes
476 for (i=i1beg; i<=i1end; i++) {
477 j += (r[i] == 0.) ? 1 : nVphi;
478 }
479
480 bool ifSide1 = false; // internal nodes
481 bool ifSide2 = false;
482
483 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
484 j += (r[i2beg] == 0.) ? 1 : nVphi;
485 ifSide1 = true;
486 }
487
488 for(i=i2beg+1; i<i2end; i++) {
489 j += (r[i] == 0.) ? 1 : nVphi;
490 }
491
492 if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
493 if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
494 ifSide2 = true;
495 }
496
497 // C O U N T F A C E S
498
499 k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi; // external faces
500
501 if (absNp2 > 1) { // internal faces
502 for(i=i2beg; i<i2end; i++) {
503 if (r[i] > 0. || r[i+1] > 0.) k += nSphi;
504 }
505
506 if (ifClosed) {
507 if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
508 }
509 }
510
511 if (!ifClosed) { // side faces
512 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
513 if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
514 }
515
516 if (!ifWholeCircle) { // phi_side faces
517 k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
518 }
519
520 // A L L O C A T E M E M O R Y
521
522 AllocateMemory(j, k);
523
524 // G E N E R A T E V E R T E C E S
525
526 int *kk;
527 kk = new int[absNp1+absNp2];
528
529 k = 1;
530 for(i=i1beg; i<=i1end; i++) {
531 kk[i] = k;
532 if (r[i] == 0.)
533 { pV[k++] = Point3D<double>(0, 0, z[i]); } else { k += nVphi; }
534 }
535
536 i = i2beg;
537 if (ifSide1) {
538 kk[i] = k;
539 if (r[i] == 0.)
540 { pV[k++] = Point3D<double>(0, 0, z[i]); } else { k += nVphi; }
541 }else{
542 kk[i] = kk[i1beg];
543 }
544
545 for(i=i2beg+1; i<i2end; i++) {
546 kk[i] = k;
547 if (r[i] == 0.)
548 { pV[k++] = Point3D<double>(0, 0, z[i]); } else { k += nVphi; }
549 }
550
551 if (absNp2 > 1) {
552 i = i2end;
553 if (ifSide2) {
554 kk[i] = k;
555 if (r[i] == 0.) pV[k] = Point3D<double>(0, 0, z[i]);
556 }else{
557 kk[i] = kk[i1end];
558 }
559 }
560
561 double cosPhi, sinPhi;
562
563 for(j=0; j<nVphi; j++) {
564 cosPhi = std::cos(phi+j*delPhi/nSphi);
565 sinPhi = std::sin(phi+j*delPhi/nSphi);
566 for(i=i1beg; i<=i2end; i++) {
567 if (r[i] != 0.)
568 pV[kk[i]+j] = Point3D<double>(r[i]*cosPhi,r[i]*sinPhi,z[i]);
569 }
570 }
571
572 // G E N E R A T E E X T E R N A L F A C E S
573
574 int v1,v2;
575
576 k = 1;
577 v2 = ifClosed ? nodeVis : 1;
578 for(i=i1beg; i<i1end; i++) {
579 v1 = v2;
580 if (!ifClosed && i == i1end-1) {
581 v2 = 1;
582 }else{
583 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
584 }
585 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
586 edgeVis, ifWholeCircle, nSphi, k);
587 }
588 if (ifClosed) {
589 RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
590 edgeVis, ifWholeCircle, nSphi, k);
591 }
592
593 // G E N E R A T E I N T E R N A L F A C E S
594
595 if (absNp2 > 1) {
596 v2 = ifClosed ? nodeVis : 1;
597 for(i=i2beg; i<i2end; i++) {
598 v1 = v2;
599 if (!ifClosed && i==i2end-1) {
600 v2 = 1;
601 }else{
602 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
603 }
604 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
605 edgeVis, ifWholeCircle, nSphi, k);
606 }
607 if (ifClosed) {
608 RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
609 edgeVis, ifWholeCircle, nSphi, k);
610 }
611 }
612
613 // G E N E R A T E S I D E F A C E S
614
615 if (!ifClosed) {
616 if (ifSide1) {
617 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
618 -1, ifWholeCircle, nSphi, k);
619 }
620 if (ifSide2) {
621 RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
622 -1, ifWholeCircle, nSphi, k);
623 }
624 }
625
626 // G E N E R A T E S I D E F A C E S for the case of incomplete circle
627
628 if (!ifWholeCircle) {
629
630 int ii[4], vv[4];
631
632 if (ifClosed) {
633 for (i=i1beg; i<=i1end; i++) {
634 ii[0] = i;
635 ii[3] = (i == i1end) ? i1beg : i+1;
636 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
637 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
638 vv[0] = -1;
639 vv[1] = 1;
640 vv[2] = -1;
641 vv[3] = 1;
642 SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
643 }
644 }else{
645 for (i=i1beg; i<i1end; i++) {
646 ii[0] = i;
647 ii[3] = i+1;
648 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
649 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
650 vv[0] = (i == i1beg) ? 1 : -1;
651 vv[1] = 1;
652 vv[2] = (i == i1end-1) ? 1 : -1;
653 vv[3] = 1;
654 SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
655 }
656 }
657 }
658
659 delete [] kk;
660
661 if (k-1 != nface) {
662 std::cerr
663 << "Polyhedron::RotateAroundZ: number of generated faces ("
664 << k-1 << ") is not equal to the number of allocated faces ("
665 << nface << ")"
666 << std::endl;
667 }
668}
669
670void HepPolyhedron::SetReferences()
671/***********************************************************************
672 * *
673 * Name: HepPolyhedron::SetReferences Date: 04.12.96 *
674 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
675 * *
676 * Function: For each edge set reference to neighbouring facet *
677 * *
678 ***********************************************************************/
679{
680 if (nface <= 0) return;
681
682 struct edgeListMember {
683 edgeListMember *next;
684 int v2;
685 int iface;
686 int iedge;
687 } *edgeList, *freeList, **headList;
688
689
690 // A L L O C A T E A N D I N I T I A T E L I S T S
691
692 edgeList = new edgeListMember[2*nface];
693 headList = new edgeListMember*[nvert];
694
695 int i;
696 for (i=0; i<nvert; i++) {
697 headList[i] = 0;
698 }
699 freeList = edgeList;
700 for (i=0; i<2*nface-1; i++) {
701 edgeList[i].next = &edgeList[i+1];
702 }
703 edgeList[2*nface-1].next = 0;
704
705 // L O O P A L O N G E D G E S
706
707 int iface, iedge, nedge, i1, i2, k1, k2;
708 edgeListMember *prev, *cur;
709
710 for(iface=1; iface<=nface; iface++) {
711 nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
712 for (iedge=0; iedge<nedge; iedge++) {
713 i1 = iedge;
714 i2 = (iedge < nedge-1) ? iedge+1 : 0;
715 i1 = std::abs(pF[iface].edge[i1].v);
716 i2 = std::abs(pF[iface].edge[i2].v);
717 k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
718 k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
719
720 // check head of the List corresponding to k1
721 cur = headList[k1];
722 if (cur == 0) {
723 headList[k1] = freeList;
724 freeList = freeList->next;
725 cur = headList[k1];
726 cur->next = 0;
727 cur->v2 = k2;
728 cur->iface = iface;
729 cur->iedge = iedge;
730 continue;
731 }
732
733 if (cur->v2 == k2) {
734 headList[k1] = cur->next;
735 cur->next = freeList;
736 freeList = cur;
737 pF[iface].edge[iedge].f = cur->iface;
738 pF[cur->iface].edge[cur->iedge].f = iface;
739 i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
740 i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
741 if (i1 != i2) {
742 std::cerr
743 << "Polyhedron::SetReferences: different edge visibility "
744 << iface << "/" << iedge << "/"
745 << pF[iface].edge[iedge].v << " and "
746 << cur->iface << "/" << cur->iedge << "/"
747 << pF[cur->iface].edge[cur->iedge].v
748 << std::endl;
749 }
750 continue;
751 }
752
753 // check List itself
754 for (;;) {
755 prev = cur;
756 cur = prev->next;
757 if (cur == 0) {
758 prev->next = freeList;
759 freeList = freeList->next;
760 cur = prev->next;
761 cur->next = 0;
762 cur->v2 = k2;
763 cur->iface = iface;
764 cur->iedge = iedge;
765 break;
766 }
767
768 if (cur->v2 == k2) {
769 prev->next = cur->next;
770 cur->next = freeList;
771 freeList = cur;
772 pF[iface].edge[iedge].f = cur->iface;
773 pF[cur->iface].edge[cur->iedge].f = iface;
774 i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
775 i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
776 if (i1 != i2) {
777 std::cerr
778 << "Polyhedron::SetReferences: different edge visibility "
779 << iface << "/" << iedge << "/"
780 << pF[iface].edge[iedge].v << " and "
781 << cur->iface << "/" << cur->iedge << "/"
782 << pF[cur->iface].edge[cur->iedge].v
783 << std::endl;
784 }
785 break;
786 }
787 }
788 }
789 }
790
791 // C H E C K T H A T A L L L I S T S A R E E M P T Y
792
793 for (i=0; i<nvert; i++) {
794 if (headList[i] != 0) {
795 std::cerr
796 << "Polyhedron::SetReferences: List " << i << " is not empty"
797 << std::endl;
798 }
799 }
800
801 // F R E E M E M O R Y
802
803 delete [] edgeList;
804 delete [] headList;
805}
806
807void HepPolyhedron::InvertFacets()
808/***********************************************************************
809 * *
810 * Name: HepPolyhedron::InvertFacets Date: 01.12.99 *
811 * Author: E.Chernyaev Revised: *
812 * *
813 * Function: Invert the order of the nodes in the facets *
814 * *
815 ***********************************************************************/
816{
817 if (nface <= 0) return;
818 int i, k, nnode, v[4],f[4];
819 for (i=1; i<=nface; i++) {
820 nnode = (pF[i].edge[3].v == 0) ? 3 : 4;
821 for (k=0; k<nnode; k++) {
822 v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
823 if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
824 f[k] = pF[i].edge[k].f;
825 }
826 for (k=0; k<nnode; k++) {
827 pF[i].edge[nnode-1-k].v = v[k];
828 pF[i].edge[nnode-1-k].f = f[k];
829 }
830 }
831}
832
833HepPolyhedron & HepPolyhedron::Transform(const Transform3D &t)
834/***********************************************************************
835 * *
836 * Name: HepPolyhedron::Transform Date: 01.12.99 *
837 * Author: E.Chernyaev Revised: *
838 * *
839 * Function: Make transformation of the polyhedron *
840 * *
841 ***********************************************************************/
842{
843 if (nvert > 0) {
844 for (int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
845
846 // C H E C K D E T E R M I N A N T A N D
847 // I N V E R T F A C E T S I F I T I S N E G A T I V E
848
849 Vector3D<double> d = t * Vector3D<double>(0,0,0);
850 Vector3D<double> x = t * Vector3D<double>(1,0,0) - d;
851 Vector3D<double> y = t * Vector3D<double>(0,1,0) - d;
852 Vector3D<double> z = t * Vector3D<double>(0,0,1) - d;
853 if ((x.cross(y))*z < 0) InvertFacets();
854 }
855 return *this;
856}
857
858bool HepPolyhedron::GetNextVertexIndex(int &index, int &edgeFlag) const
859/***********************************************************************
860 * *
861 * Name: HepPolyhedron::GetNextVertexIndex Date: 03.09.96 *
862 * Author: Yasuhide Sawada Revised: *
863 * *
864 * Function: *
865 * *
866 ***********************************************************************/
867{
868 static int iFace = 1;
869 static int iQVertex = 0;
870 int vIndex = pF[iFace].edge[iQVertex].v;
871
872 edgeFlag = (vIndex > 0) ? 1 : 0;
873 index = std::abs(vIndex);
874
875 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
876 iQVertex = 0;
877 if (++iFace > nface) iFace = 1;
878 return false; // Last Edge
879 }else{
880 ++iQVertex;
881 return true; // not Last Edge
882 }
883}
884
885Point3D<double> HepPolyhedron::GetVertex(int index) const
886/***********************************************************************
887 * *
888 * Name: HepPolyhedron::GetVertex Date: 03.09.96 *
889 * Author: Yasuhide Sawada Revised: 17.11.99 *
890 * *
891 * Function: Get vertex of the index. *
892 * *
893 ***********************************************************************/
894{
895 if (index <= 0 || index > nvert) {
896 std::cerr
897 << "HepPolyhedron::GetVertex: irrelevant index " << index
898 << std::endl;
899 return Point3D<double>();
900 }
901 return pV[index];
902}
903
904bool
905HepPolyhedron::GetNextVertex(Point3D<double> &vertex, int &edgeFlag) const
906/***********************************************************************
907 * *
908 * Name: HepPolyhedron::GetNextVertex Date: 22.07.96 *
909 * Author: John Allison Revised: *
910 * *
911 * Function: Get vertices of the quadrilaterals in order for each *
912 * face in face order. Returns false when finished each *
913 * face. *
914 * *
915 ***********************************************************************/
916{
917 int index;
918 bool rep = GetNextVertexIndex(index, edgeFlag);
919 vertex = pV[index];
920 return rep;
921}
922
923bool HepPolyhedron::GetNextVertex(Point3D<double> &vertex, int &edgeFlag,
924 Normal3D<double> &normal) const
925/***********************************************************************
926 * *
927 * Name: HepPolyhedron::GetNextVertex Date: 26.11.99 *
928 * Author: E.Chernyaev Revised: *
929 * *
930 * Function: Get vertices with normals of the quadrilaterals in order *
931 * for each face in face order. *
932 * Returns false when finished each face. *
933 * *
934 ***********************************************************************/
935{
936 static int iFace = 1;
937 static int iNode = 0;
938
939 if (nface == 0) return false; // empty polyhedron
940
941 int k = pF[iFace].edge[iNode].v;
942 if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
943 vertex = pV[k];
944 normal = FindNodeNormal(iFace,k);
945 if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
946 iNode = 0;
947 if (++iFace > nface) iFace = 1;
948 return false; // last node
949 }else{
950 ++iNode;
951 return true; // not last node
952 }
953}
954
955bool HepPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag,
956 int &iface1, int &iface2) const
957/***********************************************************************
958 * *
959 * Name: HepPolyhedron::GetNextEdgeIndeces Date: 30.09.96 *
960 * Author: E.Chernyaev Revised: 17.11.99 *
961 * *
962 * Function: Get indeces of the next edge together with indeces of *
963 * of the faces which share the edge. *
964 * Returns false when the last edge. *
965 * *
966 ***********************************************************************/
967{
968 static int iFace = 1;
969 static int iQVertex = 0;
970 static int iOrder = 1;
971 int k1, k2, kflag, kface1, kface2;
972
973 if (iFace == 1 && iQVertex == 0) {
974 k2 = pF[nface].edge[0].v;
975 k1 = pF[nface].edge[3].v;
976 if (k1 == 0) k1 = pF[nface].edge[2].v;
977 if (std::abs(k1) > std::abs(k2)) iOrder = -1;
978 }
979
980 do {
981 k1 = pF[iFace].edge[iQVertex].v;
982 kflag = k1;
983 k1 = std::abs(k1);
984 kface1 = iFace;
985 kface2 = pF[iFace].edge[iQVertex].f;
986 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
987 iQVertex = 0;
988 k2 = std::abs(pF[iFace].edge[iQVertex].v);
989 iFace++;
990 }else{
991 iQVertex++;
992 k2 = std::abs(pF[iFace].edge[iQVertex].v);
993 }
994 } while (iOrder*k1 > iOrder*k2);
995
996 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
997 iface1 = kface1; iface2 = kface2;
998
999 if (iFace > nface) {
1000 iFace = 1; iOrder = 1;
1001 return false;
1002 }else{
1003 return true;
1004 }
1005}
1006
1007bool
1008HepPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag) const
1009/***********************************************************************
1010 * *
1011 * Name: HepPolyhedron::GetNextEdgeIndeces Date: 17.11.99 *
1012 * Author: E.Chernyaev Revised: *
1013 * *
1014 * Function: Get indeces of the next edge. *
1015 * Returns false when the last edge. *
1016 * *
1017 ***********************************************************************/
1018{
1019 int kface1, kface2;
1020 return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1021}
1022
1023bool
1024HepPolyhedron::GetNextEdge(Point3D<double> &p1,
1025 Point3D<double> &p2,
1026 int &edgeFlag) const
1027/***********************************************************************
1028 * *
1029 * Name: HepPolyhedron::GetNextEdge Date: 30.09.96 *
1030 * Author: E.Chernyaev Revised: *
1031 * *
1032 * Function: Get next edge. *
1033 * Returns false when the last edge. *
1034 * *
1035 ***********************************************************************/
1036{
1037 int i1,i2;
1038 bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1039 p1 = pV[i1];
1040 p2 = pV[i2];
1041 return rep;
1042}
1043
1044bool
1045HepPolyhedron::GetNextEdge(Point3D<double> &p1, Point3D<double> &p2,
1046 int &edgeFlag, int &iface1, int &iface2) const
1047/***********************************************************************
1048 * *
1049 * Name: HepPolyhedron::GetNextEdge Date: 17.11.99 *
1050 * Author: E.Chernyaev Revised: *
1051 * *
1052 * Function: Get next edge with indeces of the faces which share *
1053 * the edge. *
1054 * Returns false when the last edge. *
1055 * *
1056 ***********************************************************************/
1057{
1058 int i1,i2;
1059 bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1060 p1 = pV[i1];
1061 p2 = pV[i2];
1062 return rep;
1063}
1064
1065void HepPolyhedron::GetFacet(int iFace, int &n, int *iNodes,
1066 int *edgeFlags, int *iFaces) const
1067/***********************************************************************
1068 * *
1069 * Name: HepPolyhedron::GetFacet Date: 15.12.99 *
1070 * Author: E.Chernyaev Revised: *
1071 * *
1072 * Function: Get face by index *
1073 * *
1074 ***********************************************************************/
1075{
1076 if (iFace < 1 || iFace > nface) {
1077 std::cerr
1078 << "HepPolyhedron::GetFacet: irrelevant index " << iFace
1079 << std::endl;
1080 n = 0;
1081 }else{
1082 int i, k;
1083 for (i=0; i<4; i++) {
1084 k = pF[iFace].edge[i].v;
1085 if (k == 0) break;
1086 if (iFaces != 0) iFaces[i] = pF[iFace].edge[i].f;
1087 if (k > 0) {
1088 iNodes[i] = k;
1089 if (edgeFlags != 0) edgeFlags[i] = 1;
1090 }else{
1091 iNodes[i] = -k;
1092 if (edgeFlags != 0) edgeFlags[i] = -1;
1093 }
1094 }
1095 n = i;
1096 }
1097}
1098
1099void HepPolyhedron::GetFacet(int index, int &n, Point3D<double> *nodes,
1100 int *edgeFlags, Normal3D<double> *normals) const
1101/***********************************************************************
1102 * *
1103 * Name: HepPolyhedron::GetFacet Date: 17.11.99 *
1104 * Author: E.Chernyaev Revised: *
1105 * *
1106 * Function: Get face by index *
1107 * *
1108 ***********************************************************************/
1109{
1110 int iNodes[4];
1111 GetFacet(index, n, iNodes, edgeFlags);
1112 if (n != 0) {
1113 for (int i=0; i<n; i++) {
1114 nodes[i] = pV[iNodes[i]];
1115 if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1116 }
1117 }
1118}
1119
1120bool
1121HepPolyhedron::GetNextFacet(int &n, Point3D<double> *nodes,
1122 int *edgeFlags, Normal3D<double> *normals) const
1123/***********************************************************************
1124 * *
1125 * Name: HepPolyhedron::GetNextFacet Date: 19.11.99 *
1126 * Author: E.Chernyaev Revised: *
1127 * *
1128 * Function: Get next face with normals of unit length at the nodes. *
1129 * Returns false when finished all faces. *
1130 * *
1131 ***********************************************************************/
1132{
1133 static int iFace = 1;
1134
1135 if (edgeFlags == 0) {
1136 GetFacet(iFace, n, nodes);
1137 }else if (normals == 0) {
1138 GetFacet(iFace, n, nodes, edgeFlags);
1139 }else{
1140 GetFacet(iFace, n, nodes, edgeFlags, normals);
1141 }
1142
1143 if (++iFace > nface) {
1144 iFace = 1;
1145 return false;
1146 }else{
1147 return true;
1148 }
1149}
1150
1151Normal3D<double> HepPolyhedron::GetNormal(int iFace) const
1152/***********************************************************************
1153 * *
1154 * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1155 * Author: E.Chernyaev Revised: *
1156 * *
1157 * Function: Get normal of the face given by index *
1158 * *
1159 ***********************************************************************/
1160{
1161 if (iFace < 1 || iFace > nface) {
1162 std::cerr
1163 << "HepPolyhedron::GetNormal: irrelevant index " << iFace
1164 << std::endl;
1165 return Normal3D<double>();
1166 }
1167
1168 int i0 = std::abs(pF[iFace].edge[0].v);
1169 int i1 = std::abs(pF[iFace].edge[1].v);
1170 int i2 = std::abs(pF[iFace].edge[2].v);
1171 int i3 = std::abs(pF[iFace].edge[3].v);
1172 if (i3 == 0) i3 = i0;
1173 return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1174}
1175
1176Normal3D<double> HepPolyhedron::GetUnitNormal(int iFace) const
1177/***********************************************************************
1178 * *
1179 * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1180 * Author: E.Chernyaev Revised: *
1181 * *
1182 * Function: Get unit normal of the face given by index *
1183 * *
1184 ***********************************************************************/
1185{
1186 if (iFace < 1 || iFace > nface) {
1187 std::cerr
1188 << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1189 << std::endl;
1190 return Normal3D<double>();
1191 }
1192
1193 int i0 = std::abs(pF[iFace].edge[0].v);
1194 int i1 = std::abs(pF[iFace].edge[1].v);
1195 int i2 = std::abs(pF[iFace].edge[2].v);
1196 int i3 = std::abs(pF[iFace].edge[3].v);
1197 if (i3 == 0) i3 = i0;
1198 return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1199}
1200
1201bool HepPolyhedron::GetNextNormal(Normal3D<double> &normal) const
1202/***********************************************************************
1203 * *
1204 * Name: HepPolyhedron::GetNextNormal Date: 22.07.96 *
1205 * Author: John Allison Revised: 19.11.99 *
1206 * *
1207 * Function: Get normals of each face in face order. Returns false *
1208 * when finished all faces. *
1209 * *
1210 ***********************************************************************/
1211{
1212 static int iFace = 1;
1213 normal = GetNormal(iFace);
1214 if (++iFace > nface) {
1215 iFace = 1;
1216 return false;
1217 }else{
1218 return true;
1219 }
1220}
1221
1222bool HepPolyhedron::GetNextUnitNormal(Normal3D<double> &normal) const
1223/***********************************************************************
1224 * *
1225 * Name: HepPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1226 * Author: E.Chernyaev Revised: *
1227 * *
1228 * Function: Get normals of unit length of each face in face order. *
1229 * Returns false when finished all faces. *
1230 * *
1231 ***********************************************************************/
1232{
1233 bool rep = GetNextNormal(normal);
1234 normal = normal.unit();
1235 return rep;
1236}
1237
1238double HepPolyhedron::GetSurfaceArea() const
1239/***********************************************************************
1240 * *
1241 * Name: HepPolyhedron::GetSurfaceArea Date: 25.05.01 *
1242 * Author: E.Chernyaev Revised: *
1243 * *
1244 * Function: Returns area of the surface of the polyhedron. *
1245 * *
1246 ***********************************************************************/
1247{
1248 double s = 0.;
1249 for (int iFace=1; iFace<=nface; iFace++) {
1250 int i0 = std::abs(pF[iFace].edge[0].v);
1251 int i1 = std::abs(pF[iFace].edge[1].v);
1252 int i2 = std::abs(pF[iFace].edge[2].v);
1253 int i3 = std::abs(pF[iFace].edge[3].v);
1254 if (i3 == 0) i3 = i0;
1255 s += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).mag();
1256 }
1257 return s/2.;
1258}
1259
1260double HepPolyhedron::GetVolume() const
1261/***********************************************************************
1262 * *
1263 * Name: HepPolyhedron::GetVolume Date: 25.05.01 *
1264 * Author: E.Chernyaev Revised: *
1265 * *
1266 * Function: Returns volume of the polyhedron. *
1267 * *
1268 ***********************************************************************/
1269{
1270 double v = 0.;
1271 for (int iFace=1; iFace<=nface; iFace++) {
1272 int i0 = std::abs(pF[iFace].edge[0].v);
1273 int i1 = std::abs(pF[iFace].edge[1].v);
1274 int i2 = std::abs(pF[iFace].edge[2].v);
1275 int i3 = std::abs(pF[iFace].edge[3].v);
1276 Point3D<double> g;
1277 if (i3 == 0) {
1278 i3 = i0;
1279 g = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
1280 }else{
1281 g = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
1282 }
1283 v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).dot(g);
1284 }
1285 return v/6.;
1286}
1287
1288int
1289HepPolyhedron::createTwistedTrap(double Dz,
1290 const double xy1[][2],
1291 const double xy2[][2])
1292/***********************************************************************
1293 * *
1294 * Name: createTwistedTrap Date: 05.11.02 *
1295 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1296 * *
1297 * Function: Creates polyhedron for twisted trapezoid *
1298 * *
1299 * Input: Dz - half-length along Z 8----7 *
1300 * xy1[2,4] - quadrilateral at Z=-Dz 5----6 ! *
1301 * xy2[2,4] - quadrilateral at Z=+Dz ! 4-!--3 *
1302 * 1----2 *
1303 * *
1304 ***********************************************************************/
1305{
1306 AllocateMemory(12,18);
1307
1308 pV[ 1] = Point3D<double>(xy1[0][0],xy1[0][1],-Dz);
1309 pV[ 2] = Point3D<double>(xy1[1][0],xy1[1][1],-Dz);
1310 pV[ 3] = Point3D<double>(xy1[2][0],xy1[2][1],-Dz);
1311 pV[ 4] = Point3D<double>(xy1[3][0],xy1[3][1],-Dz);
1312
1313 pV[ 5] = Point3D<double>(xy2[0][0],xy2[0][1], Dz);
1314 pV[ 6] = Point3D<double>(xy2[1][0],xy2[1][1], Dz);
1315 pV[ 7] = Point3D<double>(xy2[2][0],xy2[2][1], Dz);
1316 pV[ 8] = Point3D<double>(xy2[3][0],xy2[3][1], Dz);
1317
1318 pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
1319 pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
1320 pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
1321 pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
1322
1323 enum {DUMMY, BOTTOM,
1324 LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP, LEFT_BACK,
1325 BACK_BOTTOM, BACK_LEFT, BACK_TOP, BACK_RIGHT,
1326 RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP, RIGHT_FRONT,
1327 FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP, FRONT_LEFT,
1328 TOP};
1329
1330 pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
1331
1332 pF[ 2]=G4Facet(4,BOTTOM, -1,LEFT_FRONT, -12,LEFT_BACK, 0,0);
1333 pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP, -12,LEFT_BOTTOM, 0,0);
1334 pF[ 4]=G4Facet(5,TOP, -8,LEFT_BACK, -12,LEFT_FRONT, 0,0);
1335 pF[ 5]=G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM, -12,LEFT_TOP, 0,0);
1336
1337 pF[ 6]=G4Facet(3,BOTTOM, -4,BACK_LEFT, -11,BACK_RIGHT, 0,0);
1338 pF[ 7]=G4Facet(4,LEFT_BACK, -8,BACK_TOP, -11,BACK_BOTTOM, 0,0);
1339 pF[ 8]=G4Facet(8,TOP, -7,BACK_RIGHT, -11,BACK_LEFT, 0,0);
1340 pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP, 0,0);
1341
1342 pF[10]=G4Facet(2,BOTTOM, -3,RIGHT_BACK, -10,RIGHT_FRONT, 0,0);
1343 pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP, -10,RIGHT_BOTTOM, 0,0);
1344 pF[12]=G4Facet(7,TOP, -6,RIGHT_FRONT, -10,RIGHT_BACK, 0,0);
1345 pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP, 0,0);
1346
1347 pF[14]=G4Facet(1,BOTTOM, -2,FRONT_RIGHT, -9,FRONT_LEFT, 0,0);
1348 pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP, -9,FRONT_BOTTOM, 0,0);
1349 pF[16]=G4Facet(6,TOP, -5,FRONT_LEFT, -9,FRONT_RIGHT, 0,0);
1350 pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP, 0,0);
1351
1352 pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1353
1354 return 0;
1355}
1356
1357int
1358HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
1359 const double xyz[][3],
1360 const int faces[][4])
1361/***********************************************************************
1362 * *
1363 * Name: createPolyhedron Date: 05.11.02 *
1364 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1365 * *
1366 * Function: Creates user defined polyhedron *
1367 * *
1368 * Input: Nnodes - number of nodes *
1369 * Nfaces - number of faces *
1370 * nodes[][3] - node coordinates *
1371 * faces[][4] - faces *
1372 * *
1373 ***********************************************************************/
1374{
1375 AllocateMemory(Nnodes, Nfaces);
1376 if (nvert == 0) return 1;
1377
1378 for (int i=0; i<Nnodes; i++) {
1379 pV[i+1] = Point3D<double>(xyz[i][0], xyz[i][1], xyz[i][2]);
1380 }
1381 for (int k=0; k<Nfaces; k++) {
1382 pF[k+1] = G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
1383 }
1384 SetReferences();
1385 return 0;
1386}
1387
1388HepPolyhedronTrd2::HepPolyhedronTrd2(double Dx1, double Dx2,
1389 double Dy1, double Dy2,
1390 double Dz)
1391/***********************************************************************
1392 * *
1393 * Name: HepPolyhedronTrd2 Date: 22.07.96 *
1394 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1395 * *
1396 * Function: Create GEANT4 TRD2-trapezoid *
1397 * *
1398 * Input: Dx1 - half-length along X at -Dz 8----7 *
1399 * Dx2 - half-length along X ay +Dz 5----6 ! *
1400 * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1401 * Dy2 - half-length along Y ay +Dz 1----2 *
1402 * Dz - half-length along Z *
1403 * *
1404 ***********************************************************************/
1405{
1406 AllocateMemory(8,6);
1407
1408 pV[1] = Point3D<double>(-Dx1,-Dy1,-Dz);
1409 pV[2] = Point3D<double>( Dx1,-Dy1,-Dz);
1410 pV[3] = Point3D<double>( Dx1, Dy1,-Dz);
1411 pV[4] = Point3D<double>(-Dx1, Dy1,-Dz);
1412 pV[5] = Point3D<double>(-Dx2,-Dy2, Dz);
1413 pV[6] = Point3D<double>( Dx2,-Dy2, Dz);
1414 pV[7] = Point3D<double>( Dx2, Dy2, Dz);
1415 pV[8] = Point3D<double>(-Dx2, Dy2, Dz);
1416
1417 CreatePrism();
1418}
1419
1420HepPolyhedronTrd2::~HepPolyhedronTrd2() {}
1421
1422HepPolyhedronTrd1::HepPolyhedronTrd1(double Dx1, double Dx2,
1423 double Dy, double Dz)
1424 : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1425
1426HepPolyhedronTrd1::~HepPolyhedronTrd1() {}
1427
1428HepPolyhedronBox::HepPolyhedronBox(double Dx, double Dy, double Dz)
1429 : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1430
1431HepPolyhedronBox::~HepPolyhedronBox() {}
1432
1433HepPolyhedronTrap::HepPolyhedronTrap(double Dz,
1434 double Theta,
1435 double Phi,
1436 double Dy1,
1437 double Dx1,
1438 double Dx2,
1439 double Alp1,
1440 double Dy2,
1441 double Dx3,
1442 double Dx4,
1443 double Alp2)
1444/***********************************************************************
1445 * *
1446 * Name: HepPolyhedronTrap Date: 20.11.96 *
1447 * Author: E.Chernyaev Revised: *
1448 * *
1449 * Function: Create GEANT4 TRAP-trapezoid *
1450 * *
1451 * Input: DZ - half-length in Z *
1452 * Theta,Phi - polar angles of the line joining centres of the *
1453 * faces at Z=-Dz and Z=+Dz *
1454 * Dy1 - half-length in Y of the face at Z=-Dz *
1455 * Dx1 - half-length in X of low edge of the face at Z=-Dz *
1456 * Dx2 - half-length in X of top edge of the face at Z=-Dz *
1457 * Alp1 - angle between Y-axis and the median joining top and *
1458 * low edges of the face at Z=-Dz *
1459 * Dy2 - half-length in Y of the face at Z=+Dz *
1460 * Dx3 - half-length in X of low edge of the face at Z=+Dz *
1461 * Dx4 - half-length in X of top edge of the face at Z=+Dz *
1462 * Alp2 - angle between Y-axis and the median joining top and *
1463 * low edges of the face at Z=+Dz *
1464 * *
1465 ***********************************************************************/
1466{
1467 double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
1468 double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
1469 double Dy1Talp1 = Dy1*std::tan(Alp1);
1470 double Dy2Talp2 = Dy2*std::tan(Alp2);
1471
1472 AllocateMemory(8,6);
1473
1474 pV[1] = Point3D<double>(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
1475 pV[2] = Point3D<double>(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
1476 pV[3] = Point3D<double>(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
1477 pV[4] = Point3D<double>(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
1478 pV[5] = Point3D<double>( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
1479 pV[6] = Point3D<double>( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
1480 pV[7] = Point3D<double>( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
1481 pV[8] = Point3D<double>( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
1482
1483 CreatePrism();
1484}
1485
1486HepPolyhedronTrap::~HepPolyhedronTrap() {}
1487
1488HepPolyhedronPara::HepPolyhedronPara(double Dx, double Dy, double Dz,
1489 double Alpha, double Theta,
1490 double Phi)
1491 : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1492
1493HepPolyhedronPara::~HepPolyhedronPara() {}
1494
1495HepPolyhedronParaboloid::HepPolyhedronParaboloid(double r1,
1496 double r2,
1497 double dz,
1498 double sPhi,
1499 double dPhi)
1500/***********************************************************************
1501 * *
1502 * Name: HepPolyhedronParaboloid Date: 28.06.07 *
1503 * Author: L.Lindroos, T.Nikitina (CERN), July 2007 Revised: 28.06.07 *
1504 * *
1505 * Function: Constructor for paraboloid *
1506 * *
1507 * Input: r1 - inside and outside radiuses at -Dz *
1508 * r2 - inside and outside radiuses at +Dz *
1509 * dz - half length in Z *
1510 * sPhi - starting angle of the segment *
1511 * dPhi - segment range *
1512 * *
1513 ***********************************************************************/
1514{
1515 static double wholeCircle=twopi;
1516
1517 // C H E C K I N P U T P A R A M E T E R S
1518
1519 int k = 0;
1520 if (r1 < 0. || r2 <= 0.) k = 1;
1521
1522 if (dz <= 0.) k += 2;
1523
1524 double phi1, phi2, dphi;
1525
1526 if(dPhi < 0.)
1527 {
1528 phi2 = sPhi; phi1 = phi2 + dPhi;
1529 }
1530 else if(dPhi == 0.)
1531 {
1532 phi1 = sPhi; phi2 = phi1 + wholeCircle;
1533 }
1534 else
1535 {
1536 phi1 = sPhi; phi2 = phi1 + dPhi;
1537 }
1538 dphi = phi2 - phi1;
1539
1540 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1541 if (dphi > wholeCircle) k += 4;
1542
1543 if (k != 0) {
1544 std::cerr << "HepPolyhedronParaboloid: error in input parameters";
1545 if ((k & 1) != 0) std::cerr << " (radiuses)";
1546 if ((k & 2) != 0) std::cerr << " (half-length)";
1547 if ((k & 4) != 0) std::cerr << " (angles)";
1548 std::cerr << std::endl;
1549 std::cerr << " r1=" << r1;
1550 std::cerr << " r2=" << r2;
1551 std::cerr << " dz=" << dz << " sPhi=" << sPhi << " dPhi=" << dPhi
1552 << std::endl;
1553 return;
1554 }
1555
1556 // P R E P A R E T W O P O L Y L I N E S
1557
1558 int n = GetNumberOfRotationSteps();
1559 double dl = (r2 - r1) / n;
1560 double k1 = (r2*r2 - r1*r1) / 2 / dz;
1561 double k2 = (r2*r2 + r1*r1) / 2;
1562
1563 double *zz = new double[n + 2], *rr = new double[n + 2];
1564
1565 zz[0] = dz;
1566 rr[0] = r2;
1567
1568 for(int i = 1; i < n - 1; i++)
1569 {
1570 rr[i] = rr[i-1] - dl;
1571 zz[i] = (rr[i]*rr[i] - k2) / k1;
1572 if(rr[i] < 0)
1573 {
1574 rr[i] = 0;
1575 zz[i] = 0;
1576 }
1577 }
1578
1579 zz[n-1] = -dz;
1580 rr[n-1] = r1;
1581
1582 zz[n] = dz;
1583 rr[n] = 0;
1584
1585 zz[n+1] = -dz;
1586 rr[n+1] = 0;
1587
1588 // R O T A T E P O L Y L I N E S
1589
1590 RotateAroundZ(0, phi1, dphi, n, 2, zz, rr, -1, -1);
1591 SetReferences();
1592
1593 delete zz;
1594 delete rr;
1595}
1596
1597HepPolyhedronParaboloid::~HepPolyhedronParaboloid() {}
1598
1599HepPolyhedronHype::HepPolyhedronHype(double r1,
1600 double r2,
1601 double sqrtan1,
1602 double sqrtan2,
1603 double halfZ)
1604/***********************************************************************
1605 * *
1606 * Name: HepPolyhedronHype Date: 14.04.08 *
1607 * Author: Tatiana Nikitina (CERN) Revised: 14.04.08 *
1608 * *
1609 * Function: Constructor for Hype *
1610 * *
1611 * Input: r1 - inside radius at z=0 *
1612 * r2 - outside radiuses at z=0 *
1613 * sqrtan1 - sqr of tan of Inner Stereo Angle *
1614 * sqrtan2 - sqr of tan of Outer Stereo Angle *
1615 * halfZ - half length in Z *
1616 * *
1617 ***********************************************************************/
1618{
1619 static double wholeCircle=twopi;
1620
1621 // C H E C K I N P U T P A R A M E T E R S
1622
1623 int k = 0;
1624 if (r2 < 0. || r1 < 0. ) k = 1;
1625 if (r1 > r2 ) k = 1;
1626 if (r1 == r2) k = 1;
1627
1628 if (halfZ <= 0.) k += 2;
1629
1630 if (sqrtan1<0.||sqrtan2<0.) k += 4;
1631
1632 if (k != 0)
1633 {
1634 std::cerr << "HepPolyhedronHype: error in input parameters";
1635 if ((k & 1) != 0) std::cerr << " (radiuses)";
1636 if ((k & 2) != 0) std::cerr << " (half-length)";
1637 if ((k & 4) != 0) std::cerr << " (angles)";
1638 std::cerr << std::endl;
1639 std::cerr << " r1=" << r1 << " r2=" << r2;
1640 std::cerr << " halfZ=" << halfZ << " sqrTan1=" << sqrtan1
1641 << " sqrTan2=" << sqrtan2
1642 << std::endl;
1643 return;
1644 }
1645
1646 // P R E P A R E T W O P O L Y L I N E S
1647
1648 int n = GetNumberOfRotationSteps();
1649 double dz = 2.*halfZ / n;
1650 double k1 = r1*r1;
1651 double k2 = r2*r2;
1652
1653 double *zz = new double[n+n+1], *rr = new double[n+n+1];
1654
1655 zz[0] = halfZ;
1656 rr[0] = std::sqrt(sqrtan2*halfZ*halfZ+k2);
1657
1658 for(int i = 1; i < n-1; i++)
1659 {
1660 zz[i] = zz[i-1] - dz;
1661 rr[i] =std::sqrt(sqrtan2*zz[i]*zz[i]+k2);
1662 }
1663
1664 zz[n-1] = -halfZ;
1665 rr[n-1] = rr[0];
1666
1667 zz[n] = halfZ;
1668 rr[n] = std::sqrt(sqrtan1*halfZ*halfZ+k1);
1669
1670 for(int i = n+1; i < n+n; i++)
1671 {
1672 zz[i] = zz[i-1] - dz;
1673 rr[i] =std::sqrt(sqrtan1*zz[i]*zz[i]+k1);
1674 }
1675 zz[n+n] = -halfZ;
1676 rr[n+n] = rr[n];
1677
1678 // R O T A T E P O L Y L I N E S
1679
1680 RotateAroundZ(0, 0., wholeCircle, n, n, zz, rr, -1, -1);
1681 SetReferences();
1682}
1683
1684HepPolyhedronHype::~HepPolyhedronHype() {}
1685
1686HepPolyhedronCons::HepPolyhedronCons(double Rmn1,
1687 double Rmx1,
1688 double Rmn2,
1689 double Rmx2,
1690 double Dz,
1691 double Phi1,
1692 double Dphi)
1693/***********************************************************************
1694 * *
1695 * Name: HepPolyhedronCons::HepPolyhedronCons Date: 15.12.96 *
1696 * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
1697 * *
1698 * Function: Constructor for CONS, TUBS, CONE, TUBE *
1699 * *
1700 * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
1701 * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
1702 * Dz - half length in Z *
1703 * Phi1 - starting angle of the segment *
1704 * Dphi - segment range *
1705 * *
1706 ***********************************************************************/
1707{
1708 static double wholeCircle=twopi;
1709
1710 // C H E C K I N P U T P A R A M E T E R S
1711
1712 int k = 0;
1713 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
1714 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
1715 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
1716
1717 if (Dz <= 0.) k += 2;
1718
1719 double phi1, phi2, dphi;
1720 if (Dphi < 0.) {
1721 phi2 = Phi1; phi1 = phi2 - Dphi;
1722 }else if (Dphi == 0.) {
1723 phi1 = Phi1; phi2 = phi1 + wholeCircle;
1724 }else{
1725 phi1 = Phi1; phi2 = phi1 + Dphi;
1726 }
1727 dphi = phi2 - phi1;
1728 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1729 if (dphi > wholeCircle) k += 4;
1730
1731 if (k != 0) {
1732 std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
1733 if ((k & 1) != 0) std::cerr << " (radiuses)";
1734 if ((k & 2) != 0) std::cerr << " (half-length)";
1735 if ((k & 4) != 0) std::cerr << " (angles)";
1736 std::cerr << std::endl;
1737 std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
1738 std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
1739 std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
1740 << std::endl;
1741 return;
1742 }
1743
1744 // P R E P A R E T W O P O L Y L I N E S
1745
1746 double zz[4], rr[4];
1747 zz[0] = Dz;
1748 zz[1] = -Dz;
1749 zz[2] = Dz;
1750 zz[3] = -Dz;
1751 rr[0] = Rmx2;
1752 rr[1] = Rmx1;
1753 rr[2] = Rmn2;
1754 rr[3] = Rmn1;
1755
1756 // R O T A T E P O L Y L I N E S
1757
1758 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
1759 SetReferences();
1760}
1761
1762HepPolyhedronCons::~HepPolyhedronCons() {}
1763
1764HepPolyhedronCone::HepPolyhedronCone(double Rmn1, double Rmx1,
1765 double Rmn2, double Rmx2,
1766 double Dz) :
1767 HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
1768
1769HepPolyhedronCone::~HepPolyhedronCone() {}
1770
1771HepPolyhedronTubs::HepPolyhedronTubs(double Rmin, double Rmax,
1772 double Dz,
1773 double Phi1, double Dphi)
1774 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1775
1776HepPolyhedronTubs::~HepPolyhedronTubs() {}
1777
1778HepPolyhedronTube::HepPolyhedronTube (double Rmin, double Rmax,
1779 double Dz)
1780 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
1781
1782HepPolyhedronTube::~HepPolyhedronTube () {}
1783
1784HepPolyhedronPgon::HepPolyhedronPgon(double phi,
1785 double dphi,
1786 int npdv,
1787 int nz,
1788 const double *z,
1789 const double *rmin,
1790 const double *rmax)
1791/***********************************************************************
1792 * *
1793 * Name: HepPolyhedronPgon Date: 09.12.96 *
1794 * Author: E.Chernyaev Revised: *
1795 * *
1796 * Function: Constructor of polyhedron for PGON, PCON *
1797 * *
1798 * Input: phi - initial phi *
1799 * dphi - delta phi *
1800 * npdv - number of steps along phi *
1801 * nz - number of z-planes (at least two) *
1802 * z[] - z coordinates of the slices *
1803 * rmin[] - smaller r at the slices *
1804 * rmax[] - bigger r at the slices *
1805 * *
1806 ***********************************************************************/
1807{
1808 // C H E C K I N P U T P A R A M E T E R S
1809
1810 if (dphi <= 0. || dphi > twopi) {
1811 std::cerr
1812 << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1813 << std::endl;
1814 return;
1815 }
1816
1817 if (nz < 2) {
1818 std::cerr
1819 << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1820 << std::endl;
1821 return;
1822 }
1823
1824 if (npdv < 0) {
1825 std::cerr
1826 << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1827 << std::endl;
1828 return;
1829 }
1830
1831 int i;
1832 for (i=0; i<nz; i++) {
1833 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1834 std::cerr
1835 << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
1836 << rmin[i] << " rmax[" << i << "]=" << rmax[i]
1837 << std::endl;
1838 return;
1839 }
1840 }
1841
1842 // P R E P A R E T W O P O L Y L I N E S
1843
1844 double *zz, *rr;
1845 zz = new double[2*nz];
1846 rr = new double[2*nz];
1847
1848 if (z[0] > z[nz-1]) {
1849 for (i=0; i<nz; i++) {
1850 zz[i] = z[i];
1851 rr[i] = rmax[i];
1852 zz[i+nz] = z[i];
1853 rr[i+nz] = rmin[i];
1854 }
1855 }else{
1856 for (i=0; i<nz; i++) {
1857 zz[i] = z[nz-i-1];
1858 rr[i] = rmax[nz-i-1];
1859 zz[i+nz] = z[nz-i-1];
1860 rr[i+nz] = rmin[nz-i-1];
1861 }
1862 }
1863
1864 // R O T A T E P O L Y L I N E S
1865
1866 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1867 SetReferences();
1868
1869 delete [] zz;
1870 delete [] rr;
1871}
1872
1873HepPolyhedronPgon::~HepPolyhedronPgon() {}
1874
1875HepPolyhedronPcon::HepPolyhedronPcon(double phi, double dphi, int nz,
1876 const double *z,
1877 const double *rmin,
1878 const double *rmax)
1879 : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
1880
1881HepPolyhedronPcon::~HepPolyhedronPcon() {}
1882
1883HepPolyhedronSphere::HepPolyhedronSphere(double rmin, double rmax,
1884 double phi, double dphi,
1885 double the, double dthe)
1886/***********************************************************************
1887 * *
1888 * Name: HepPolyhedronSphere Date: 11.12.96 *
1889 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1890 * *
1891 * Function: Constructor of polyhedron for SPHERE *
1892 * *
1893 * Input: rmin - internal radius *
1894 * rmax - external radius *
1895 * phi - initial phi *
1896 * dphi - delta phi *
1897 * the - initial theta *
1898 * dthe - delta theta *
1899 * *
1900 ***********************************************************************/
1901{
1902 // C H E C K I N P U T P A R A M E T E R S
1903
1904 if (dphi <= 0. || dphi > twopi) {
1905 std::cerr
1906 << "HepPolyhedronSphere: wrong delta phi = " << dphi
1907 << std::endl;
1908 return;
1909 }
1910
1911 if (the < 0. || the > pi) {
1912 std::cerr
1913 << "HepPolyhedronSphere: wrong theta = " << the
1914 << std::endl;
1915 return;
1916 }
1917
1918 if (dthe <= 0. || dthe > pi) {
1919 std::cerr
1920 << "HepPolyhedronSphere: wrong delta theta = " << dthe
1921 << std::endl;
1922 return;
1923 }
1924
1925 if (the+dthe > pi) {
1926 std::cerr
1927 << "HepPolyhedronSphere: wrong theta + delta theta = "
1928 << the << " " << dthe
1929 << std::endl;
1930 return;
1931 }
1932
1933 if (rmin < 0. || rmin >= rmax) {
1934 std::cerr
1935 << "HepPolyhedronSphere: error in radiuses"
1936 << " rmin=" << rmin << " rmax=" << rmax
1937 << std::endl;
1938 return;
1939 }
1940
1941 // P R E P A R E T W O P O L Y L I N E S
1942
1943 int ns = (GetNumberOfRotationSteps() + 1) / 2;
1944 int np1 = int(dthe*ns/pi+.5) + 1;
1945 if (np1 <= 1) np1 = 2;
1946 int np2 = rmin < perMillion ? 1 : np1;
1947
1948 double *zz, *rr;
1949 zz = new double[np1+np2];
1950 rr = new double[np1+np2];
1951
1952 double a = dthe/(np1-1);
1953 double cosa, sina;
1954 for (int i=0; i<np1; i++) {
1955 cosa = std::cos(the+i*a);
1956 sina = std::sin(the+i*a);
1957 zz[i] = rmax*cosa;
1958 rr[i] = rmax*sina;
1959 if (np2 > 1) {
1960 zz[i+np1] = rmin*cosa;
1961 rr[i+np1] = rmin*sina;
1962 }
1963 }
1964 if (np2 == 1) {
1965 zz[np1] = 0.;
1966 rr[np1] = 0.;
1967 }
1968
1969 // R O T A T E P O L Y L I N E S
1970
1971 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
1972 SetReferences();
1973
1974 delete [] zz;
1975 delete [] rr;
1976}
1977
1978HepPolyhedronSphere::~HepPolyhedronSphere() {}
1979
1980HepPolyhedronTorus::HepPolyhedronTorus(double rmin,
1981 double rmax,
1982 double rtor,
1983 double phi,
1984 double dphi)
1985/***********************************************************************
1986 * *
1987 * Name: HepPolyhedronTorus Date: 11.12.96 *
1988 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1989 * *
1990 * Function: Constructor of polyhedron for TORUS *
1991 * *
1992 * Input: rmin - internal radius *
1993 * rmax - external radius *
1994 * rtor - radius of torus *
1995 * phi - initial phi *
1996 * dphi - delta phi *
1997 * *
1998 ***********************************************************************/
1999{
2000 // C H E C K I N P U T P A R A M E T E R S
2001
2002 if (dphi <= 0. || dphi > twopi) {
2003 std::cerr
2004 << "HepPolyhedronTorus: wrong delta phi = " << dphi
2005 << std::endl;
2006 return;
2007 }
2008
2009 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2010 std::cerr
2011 << "HepPolyhedronTorus: error in radiuses"
2012 << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2013 << std::endl;
2014 return;
2015 }
2016
2017 // P R E P A R E T W O P O L Y L I N E S
2018
2019 int np1 = GetNumberOfRotationSteps();
2020 int np2 = rmin < perMillion ? 1 : np1;
2021
2022 double *zz, *rr;
2023 zz = new double[np1+np2];
2024 rr = new double[np1+np2];
2025
2026 double a = twopi/np1;
2027 double cosa, sina;
2028 for (int i=0; i<np1; i++) {
2029 cosa = std::cos(i*a);
2030 sina = std::sin(i*a);
2031 zz[i] = rmax*cosa;
2032 rr[i] = rtor+rmax*sina;
2033 if (np2 > 1) {
2034 zz[i+np1] = rmin*cosa;
2035 rr[i+np1] = rtor+rmin*sina;
2036 }
2037 }
2038 if (np2 == 1) {
2039 zz[np1] = 0.;
2040 rr[np1] = rtor;
2041 np2 = -1;
2042 }
2043
2044 // R O T A T E P O L Y L I N E S
2045
2046 RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
2047 SetReferences();
2048
2049 delete [] zz;
2050 delete [] rr;
2051}
2052
2053HepPolyhedronTorus::~HepPolyhedronTorus() {}
2054
2055HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(double ax, double by,
2056 double cz, double zCut1,
2057 double zCut2)
2058/***********************************************************************
2059 * *
2060 * Name: HepPolyhedronEllipsoid Date: 25.02.05 *
2061 * Author: G.Guerrieri Revised: *
2062 * *
2063 * Function: Constructor of polyhedron for ELLIPSOID *
2064 * *
2065 * Input: ax - semiaxis x *
2066 * by - semiaxis y *
2067 * cz - semiaxis z *
2068 * zCut1 - lower cut plane level (solid lies above this plane) *
2069 * zCut2 - upper cut plane level (solid lies below this plane) *
2070 * *
2071 ***********************************************************************/
2072{
2073 // C H E C K I N P U T P A R A M E T E R S
2074
2075 if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2076 std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2077 << " zCut2 = " << zCut2
2078 << " for given cz = " << cz << std::endl;
2079 return;
2080 }
2081 if (cz <= 0.0) {
2082 std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2083 << std::endl;
2084 return;
2085 }
2086
2087 double dthe;
2088 double sthe;
2089 int cutflag;
2090 cutflag= 0;
2091 if (zCut2 >= cz)
2092 {
2093 sthe= 0.0;
2094 }
2095 else
2096 {
2097 sthe= std::acos(zCut2/cz);
2098 cutflag++;
2099 }
2100 if (zCut1 <= -cz)
2101 {
2102 dthe= pi - sthe;
2103 }
2104 else
2105 {
2106 dthe= std::acos(zCut1/cz)-sthe;
2107 cutflag++;
2108 }
2109
2110 // P R E P A R E T W O P O L Y L I N E S
2111 // generate sphere of radius cz first, then rescale x and y later
2112
2113 int ns = (GetNumberOfRotationSteps() + 1) / 2;
2114 int np1 = int(dthe*ns/pi) + 2 + cutflag;
2115
2116 double *zz, *rr;
2117 zz = new double[np1+1];
2118 rr = new double[np1+1];
2119 if (!zz || !rr)
2120 {
2121 std::cerr << "Out of memory in HepPolyhedronEllipsoid!" << std::endl;
2122 //Exception("Out of memory in HepPolyhedronEllipsoid!");
2123 }
2124
2125 double a = dthe/(np1-cutflag-1);
2126 double cosa, sina;
2127 int j=0;
2128 if (sthe > 0.0)
2129 {
2130 zz[j]= zCut2;
2131 rr[j]= 0.;
2132 j++;
2133 }
2134 for (int i=0; i<np1-cutflag; i++) {
2135 cosa = std::cos(sthe+i*a);
2136 sina = std::sin(sthe+i*a);
2137 zz[j] = cz*cosa;
2138 rr[j] = cz*sina;
2139 j++;
2140 }
2141 if (j < np1)
2142 {
2143 zz[j]= zCut1;
2144 rr[j]= 0.;
2145 j++;
2146 }
2147 if (j > np1)
2148 {
2149 std::cerr << "Logic error in HepPolyhedronEllipsoid, memory corrupted!"
2150 << std::endl;
2151 }
2152 if (j < np1)
2153 {
2154 std::cerr << "Warning: logic error in HepPolyhedronEllipsoid."
2155 << std::endl;
2156 np1= j;
2157 }
2158 zz[j] = 0.;
2159 rr[j] = 0.;
2160
2161
2162 // R O T A T E P O L Y L I N E S
2163
2164 RotateAroundZ(0, 0.0, twopi, np1, 1, zz, rr, -1, 1);
2165 SetReferences();
2166
2167 delete [] zz;
2168 delete [] rr;
2169
2170 // rescale x and y vertex coordinates
2171 {
2172 Point3D<double> * p= pV;
2173 for (int i=0; i<nvert; i++, p++) {
2174 p->setX( p->x() * ax/cz );
2175 p->setY( p->y() * by/cz );
2176 }
2177 }
2178}
2179
2180HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
2181
2182HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(double ax,
2183 double ay,
2184 double h,
2185 double zTopCut)
2186/***********************************************************************
2187 * *
2188 * Name: HepPolyhedronEllipticalCone Date: 8.9.2005 *
2189 * Author: D.Anninos Revised: 9.9.2005 *
2190 * *
2191 * Function: Constructor for EllipticalCone *
2192 * *
2193 * Input: ax, ay - X & Y semi axes at z = 0 *
2194 * h - height of full cone *
2195 * zTopCut - Top Cut in Z Axis *
2196 * *
2197 ***********************************************************************/
2198{
2199 // C H E C K I N P U T P A R A M E T E R S
2200
2201 int k = 0;
2202 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2203
2204 if (k != 0) {
2205 std::cerr << "HepPolyhedronCone: error in input parameters";
2206 std::cerr << std::endl;
2207 return;
2208 }
2209
2210 // P R E P A R E T W O P O L Y L I N E S
2211
2212 zTopCut = (h >= zTopCut ? zTopCut : h);
2213
2214 double *zz, *rr;
2215 zz = new double[4];
2216 rr = new double[4];
2217 zz[0] = zTopCut;
2218 zz[1] = -zTopCut;
2219 zz[2] = zTopCut;
2220 zz[3] = -zTopCut;
2221 rr[0] = (h-zTopCut);
2222 rr[1] = (h+zTopCut);
2223 rr[2] = 0.;
2224 rr[3] = 0.;
2225
2226 // R O T A T E P O L Y L I N E S
2227
2228 RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, -1);
2229 SetReferences();
2230
2231 delete [] zz;
2232 delete [] rr;
2233
2234 // rescale x and y vertex coordinates
2235 {
2236 Point3D<double> * p= pV;
2237 for (int i=0; i<nvert; i++, p++) {
2238 p->setX( p->x() * ax );
2239 p->setY( p->y() * ay );
2240 }
2241 }
2242}
2243
2244HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
2245
2246int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
2247/***********************************************************************
2248 * *
2249 * Name: HepPolyhedron::fNumberOfRotationSteps Date: 24.06.97 *
2250 * Author: J.Allison (Manchester University) Revised: *
2251 * *
2252 * Function: Number of steps for whole circle *
2253 * *
2254 ***********************************************************************/
2255
2256#include "BooleanProcessor.src"
2257
2258HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const
2259/***********************************************************************
2260 * *
2261 * Name: HepPolyhedron::add Date: 19.03.00 *
2262 * Author: E.Chernyaev Revised: *
2263 * *
2264 * Function: Boolean "union" of two polyhedra *
2265 * *
2266 ***********************************************************************/
2267{
2268 int ierr;
2269 BooleanProcessor processor;
2270 return processor.execute(OP_UNION, *this, p,ierr);
2271}
2272
2273HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const
2274/***********************************************************************
2275 * *
2276 * Name: HepPolyhedron::intersect Date: 19.03.00 *
2277 * Author: E.Chernyaev Revised: *
2278 * *
2279 * Function: Boolean "intersection" of two polyhedra *
2280 * *
2281 ***********************************************************************/
2282{
2283 int ierr;
2284 BooleanProcessor processor;
2285 return processor.execute(OP_INTERSECTION, *this, p,ierr);
2286}
2287
2288HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const
2289/***********************************************************************
2290 * *
2291 * Name: HepPolyhedron::add Date: 19.03.00 *
2292 * Author: E.Chernyaev Revised: *
2293 * *
2294 * Function: Boolean "subtraction" of "p" from "this" *
2295 * *
2296 ***********************************************************************/
2297{
2298 int ierr;
2299 BooleanProcessor processor;
2300 return processor.execute(OP_SUBTRACTION, *this, p,ierr);
2301}
2302
2303//NOTE : include the code of HepPolyhedronProcessor here
2304// since there is no BooleanProcessor.h
2305
2306#undef INTERSECTION
2307
2308#include "HepPolyhedronProcessor.src"
2309
Note: See TracBrowser for help on using the repository browser.