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

Last change on this file since 872 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 84.8 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.31 2008/04/28 16:06:06 allison Exp $
28// GEANT4 tag $Name: HEAD $
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 std::cerr << "HepPolyhedronHype: error in input parameters";
1634 if ((k & 1) != 0) std::cerr << " (radiuses)";
1635 if ((k & 2) != 0) std::cerr << " (half-length)";
1636 if ((k & 4) != 0) std::cerr << " (angles)";
1637 std::cerr << std::endl;
1638 std::cerr << " r1=" << r1 << " r2=" << r2;
1639 std::cerr << " halfZ=" << halfZ << " sqrTan1=" << sqrtan1
1640 << " sqrTan2=" << sqrtan2
1641 << std::endl;
1642 return;
1643 }
1644
1645 // P R E P A R E T W O P O L Y L I N E S
1646
1647 int n = GetNumberOfRotationSteps();
1648 double dz = 2.*halfZ / n;
1649 double k1 = r1*r1;
1650 double k2 = r2*r2;
1651
1652 double *zz = new double[n + n], *rr = new double[n + n];
1653
1654 zz[0] = halfZ;
1655 rr[0] = std::sqrt(sqrtan2*halfZ*halfZ+k2);
1656
1657 for(int i = 1; i < n - 1; i++)
1658 {
1659 zz[i] = zz[i-1] - dz;
1660 rr[i] =std::sqrt(sqrtan2*zz[i]*zz[i]+k2) ;
1661
1662 }
1663
1664 zz[n-1] = -halfZ;
1665 rr[n-1] = rr[0];
1666
1667
1668 zz[n] = halfZ;
1669 rr[n] = std::sqrt(sqrtan1*halfZ*halfZ+k1);
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 }
1676 zz[n+n] = -halfZ;
1677 rr[n+n] = rr[n];
1678
1679 // R O T A T E P O L Y L I N E S
1680
1681 RotateAroundZ(0, 0., wholeCircle, n, n, zz, rr, -1, -1);
1682 SetReferences();
1683}
1684
1685HepPolyhedronHype::~HepPolyhedronHype() {}
1686
1687HepPolyhedronCons::HepPolyhedronCons(double Rmn1,
1688 double Rmx1,
1689 double Rmn2,
1690 double Rmx2,
1691 double Dz,
1692 double Phi1,
1693 double Dphi)
1694/***********************************************************************
1695 * *
1696 * Name: HepPolyhedronCons::HepPolyhedronCons Date: 15.12.96 *
1697 * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
1698 * *
1699 * Function: Constructor for CONS, TUBS, CONE, TUBE *
1700 * *
1701 * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
1702 * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
1703 * Dz - half length in Z *
1704 * Phi1 - starting angle of the segment *
1705 * Dphi - segment range *
1706 * *
1707 ***********************************************************************/
1708{
1709 static double wholeCircle=twopi;
1710
1711 // C H E C K I N P U T P A R A M E T E R S
1712
1713 int k = 0;
1714 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
1715 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
1716 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
1717
1718 if (Dz <= 0.) k += 2;
1719
1720 double phi1, phi2, dphi;
1721 if (Dphi < 0.) {
1722 phi2 = Phi1; phi1 = phi2 - Dphi;
1723 }else if (Dphi == 0.) {
1724 phi1 = Phi1; phi2 = phi1 + wholeCircle;
1725 }else{
1726 phi1 = Phi1; phi2 = phi1 + Dphi;
1727 }
1728 dphi = phi2 - phi1;
1729 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1730 if (dphi > wholeCircle) k += 4;
1731
1732 if (k != 0) {
1733 std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
1734 if ((k & 1) != 0) std::cerr << " (radiuses)";
1735 if ((k & 2) != 0) std::cerr << " (half-length)";
1736 if ((k & 4) != 0) std::cerr << " (angles)";
1737 std::cerr << std::endl;
1738 std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
1739 std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
1740 std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
1741 << std::endl;
1742 return;
1743 }
1744
1745 // P R E P A R E T W O P O L Y L I N E S
1746
1747 double zz[4], rr[4];
1748 zz[0] = Dz;
1749 zz[1] = -Dz;
1750 zz[2] = Dz;
1751 zz[3] = -Dz;
1752 rr[0] = Rmx2;
1753 rr[1] = Rmx1;
1754 rr[2] = Rmn2;
1755 rr[3] = Rmn1;
1756
1757 // R O T A T E P O L Y L I N E S
1758
1759 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
1760 SetReferences();
1761}
1762
1763HepPolyhedronCons::~HepPolyhedronCons() {}
1764
1765HepPolyhedronCone::HepPolyhedronCone(double Rmn1, double Rmx1,
1766 double Rmn2, double Rmx2,
1767 double Dz) :
1768 HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
1769
1770HepPolyhedronCone::~HepPolyhedronCone() {}
1771
1772HepPolyhedronTubs::HepPolyhedronTubs(double Rmin, double Rmax,
1773 double Dz,
1774 double Phi1, double Dphi)
1775 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1776
1777HepPolyhedronTubs::~HepPolyhedronTubs() {}
1778
1779HepPolyhedronTube::HepPolyhedronTube (double Rmin, double Rmax,
1780 double Dz)
1781 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
1782
1783HepPolyhedronTube::~HepPolyhedronTube () {}
1784
1785HepPolyhedronPgon::HepPolyhedronPgon(double phi,
1786 double dphi,
1787 int npdv,
1788 int nz,
1789 const double *z,
1790 const double *rmin,
1791 const double *rmax)
1792/***********************************************************************
1793 * *
1794 * Name: HepPolyhedronPgon Date: 09.12.96 *
1795 * Author: E.Chernyaev Revised: *
1796 * *
1797 * Function: Constructor of polyhedron for PGON, PCON *
1798 * *
1799 * Input: phi - initial phi *
1800 * dphi - delta phi *
1801 * npdv - number of steps along phi *
1802 * nz - number of z-planes (at least two) *
1803 * z[] - z coordinates of the slices *
1804 * rmin[] - smaller r at the slices *
1805 * rmax[] - bigger r at the slices *
1806 * *
1807 ***********************************************************************/
1808{
1809 // C H E C K I N P U T P A R A M E T E R S
1810
1811 if (dphi <= 0. || dphi > twopi) {
1812 std::cerr
1813 << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1814 << std::endl;
1815 return;
1816 }
1817
1818 if (nz < 2) {
1819 std::cerr
1820 << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1821 << std::endl;
1822 return;
1823 }
1824
1825 if (npdv < 0) {
1826 std::cerr
1827 << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1828 << std::endl;
1829 return;
1830 }
1831
1832 int i;
1833 for (i=0; i<nz; i++) {
1834 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1835 std::cerr
1836 << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
1837 << rmin[i] << " rmax[" << i << "]=" << rmax[i]
1838 << std::endl;
1839 return;
1840 }
1841 }
1842
1843 // P R E P A R E T W O P O L Y L I N E S
1844
1845 double *zz, *rr;
1846 zz = new double[2*nz];
1847 rr = new double[2*nz];
1848
1849 if (z[0] > z[nz-1]) {
1850 for (i=0; i<nz; i++) {
1851 zz[i] = z[i];
1852 rr[i] = rmax[i];
1853 zz[i+nz] = z[i];
1854 rr[i+nz] = rmin[i];
1855 }
1856 }else{
1857 for (i=0; i<nz; i++) {
1858 zz[i] = z[nz-i-1];
1859 rr[i] = rmax[nz-i-1];
1860 zz[i+nz] = z[nz-i-1];
1861 rr[i+nz] = rmin[nz-i-1];
1862 }
1863 }
1864
1865 // R O T A T E P O L Y L I N E S
1866
1867 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1868 SetReferences();
1869
1870 delete [] zz;
1871 delete [] rr;
1872}
1873
1874HepPolyhedronPgon::~HepPolyhedronPgon() {}
1875
1876HepPolyhedronPcon::HepPolyhedronPcon(double phi, double dphi, int nz,
1877 const double *z,
1878 const double *rmin,
1879 const double *rmax)
1880 : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
1881
1882HepPolyhedronPcon::~HepPolyhedronPcon() {}
1883
1884HepPolyhedronSphere::HepPolyhedronSphere(double rmin, double rmax,
1885 double phi, double dphi,
1886 double the, double dthe)
1887/***********************************************************************
1888 * *
1889 * Name: HepPolyhedronSphere Date: 11.12.96 *
1890 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1891 * *
1892 * Function: Constructor of polyhedron for SPHERE *
1893 * *
1894 * Input: rmin - internal radius *
1895 * rmax - external radius *
1896 * phi - initial phi *
1897 * dphi - delta phi *
1898 * the - initial theta *
1899 * dthe - delta theta *
1900 * *
1901 ***********************************************************************/
1902{
1903 // C H E C K I N P U T P A R A M E T E R S
1904
1905 if (dphi <= 0. || dphi > twopi) {
1906 std::cerr
1907 << "HepPolyhedronSphere: wrong delta phi = " << dphi
1908 << std::endl;
1909 return;
1910 }
1911
1912 if (the < 0. || the > pi) {
1913 std::cerr
1914 << "HepPolyhedronSphere: wrong theta = " << the
1915 << std::endl;
1916 return;
1917 }
1918
1919 if (dthe <= 0. || dthe > pi) {
1920 std::cerr
1921 << "HepPolyhedronSphere: wrong delta theta = " << dthe
1922 << std::endl;
1923 return;
1924 }
1925
1926 if (the+dthe > pi) {
1927 std::cerr
1928 << "HepPolyhedronSphere: wrong theta + delta theta = "
1929 << the << " " << dthe
1930 << std::endl;
1931 return;
1932 }
1933
1934 if (rmin < 0. || rmin >= rmax) {
1935 std::cerr
1936 << "HepPolyhedronSphere: error in radiuses"
1937 << " rmin=" << rmin << " rmax=" << rmax
1938 << std::endl;
1939 return;
1940 }
1941
1942 // P R E P A R E T W O P O L Y L I N E S
1943
1944 int ns = (GetNumberOfRotationSteps() + 1) / 2;
1945 int np1 = int(dthe*ns/pi+.5) + 1;
1946 if (np1 <= 1) np1 = 2;
1947 int np2 = rmin < perMillion ? 1 : np1;
1948
1949 double *zz, *rr;
1950 zz = new double[np1+np2];
1951 rr = new double[np1+np2];
1952
1953 double a = dthe/(np1-1);
1954 double cosa, sina;
1955 for (int i=0; i<np1; i++) {
1956 cosa = std::cos(the+i*a);
1957 sina = std::sin(the+i*a);
1958 zz[i] = rmax*cosa;
1959 rr[i] = rmax*sina;
1960 if (np2 > 1) {
1961 zz[i+np1] = rmin*cosa;
1962 rr[i+np1] = rmin*sina;
1963 }
1964 }
1965 if (np2 == 1) {
1966 zz[np1] = 0.;
1967 rr[np1] = 0.;
1968 }
1969
1970 // R O T A T E P O L Y L I N E S
1971
1972 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
1973 SetReferences();
1974
1975 delete [] zz;
1976 delete [] rr;
1977}
1978
1979HepPolyhedronSphere::~HepPolyhedronSphere() {}
1980
1981HepPolyhedronTorus::HepPolyhedronTorus(double rmin,
1982 double rmax,
1983 double rtor,
1984 double phi,
1985 double dphi)
1986/***********************************************************************
1987 * *
1988 * Name: HepPolyhedronTorus Date: 11.12.96 *
1989 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1990 * *
1991 * Function: Constructor of polyhedron for TORUS *
1992 * *
1993 * Input: rmin - internal radius *
1994 * rmax - external radius *
1995 * rtor - radius of torus *
1996 * phi - initial phi *
1997 * dphi - delta phi *
1998 * *
1999 ***********************************************************************/
2000{
2001 // C H E C K I N P U T P A R A M E T E R S
2002
2003 if (dphi <= 0. || dphi > twopi) {
2004 std::cerr
2005 << "HepPolyhedronTorus: wrong delta phi = " << dphi
2006 << std::endl;
2007 return;
2008 }
2009
2010 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2011 std::cerr
2012 << "HepPolyhedronTorus: error in radiuses"
2013 << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2014 << std::endl;
2015 return;
2016 }
2017
2018 // P R E P A R E T W O P O L Y L I N E S
2019
2020 int np1 = GetNumberOfRotationSteps();
2021 int np2 = rmin < perMillion ? 1 : np1;
2022
2023 double *zz, *rr;
2024 zz = new double[np1+np2];
2025 rr = new double[np1+np2];
2026
2027 double a = twopi/np1;
2028 double cosa, sina;
2029 for (int i=0; i<np1; i++) {
2030 cosa = std::cos(i*a);
2031 sina = std::sin(i*a);
2032 zz[i] = rmax*cosa;
2033 rr[i] = rtor+rmax*sina;
2034 if (np2 > 1) {
2035 zz[i+np1] = rmin*cosa;
2036 rr[i+np1] = rtor+rmin*sina;
2037 }
2038 }
2039 if (np2 == 1) {
2040 zz[np1] = 0.;
2041 rr[np1] = rtor;
2042 np2 = -1;
2043 }
2044
2045 // R O T A T E P O L Y L I N E S
2046
2047 RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
2048 SetReferences();
2049
2050 delete [] zz;
2051 delete [] rr;
2052}
2053
2054HepPolyhedronTorus::~HepPolyhedronTorus() {}
2055
2056HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(double ax, double by,
2057 double cz, double zCut1,
2058 double zCut2)
2059/***********************************************************************
2060 * *
2061 * Name: HepPolyhedronEllipsoid Date: 25.02.05 *
2062 * Author: G.Guerrieri Revised: *
2063 * *
2064 * Function: Constructor of polyhedron for ELLIPSOID *
2065 * *
2066 * Input: ax - semiaxis x *
2067 * by - semiaxis y *
2068 * cz - semiaxis z *
2069 * zCut1 - lower cut plane level (solid lies above this plane) *
2070 * zCut2 - upper cut plane level (solid lies below this plane) *
2071 * *
2072 ***********************************************************************/
2073{
2074 // C H E C K I N P U T P A R A M E T E R S
2075
2076 if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2077 std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2078 << " zCut2 = " << zCut2
2079 << " for given cz = " << cz << std::endl;
2080 return;
2081 }
2082 if (cz <= 0.0) {
2083 std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2084 << std::endl;
2085 return;
2086 }
2087
2088 double dthe;
2089 double sthe;
2090 int cutflag;
2091 cutflag= 0;
2092 if (zCut2 >= cz)
2093 {
2094 sthe= 0.0;
2095 }
2096 else
2097 {
2098 sthe= std::acos(zCut2/cz);
2099 cutflag++;
2100 }
2101 if (zCut1 <= -cz)
2102 {
2103 dthe= pi - sthe;
2104 }
2105 else
2106 {
2107 dthe= std::acos(zCut1/cz)-sthe;
2108 cutflag++;
2109 }
2110
2111 // P R E P A R E T W O P O L Y L I N E S
2112 // generate sphere of radius cz first, then rescale x and y later
2113
2114 int ns = (GetNumberOfRotationSteps() + 1) / 2;
2115 int np1 = int(dthe*ns/pi) + 2 + cutflag;
2116
2117 double *zz, *rr;
2118 zz = new double[np1+1];
2119 rr = new double[np1+1];
2120 if (!zz || !rr)
2121 {
2122 std::cerr << "Out of memory in HepPolyhedronEllipsoid!" << std::endl;
2123 //Exception("Out of memory in HepPolyhedronEllipsoid!");
2124 }
2125
2126 double a = dthe/(np1-cutflag-1);
2127 double cosa, sina;
2128 int j=0;
2129 if (sthe > 0.0)
2130 {
2131 zz[j]= zCut2;
2132 rr[j]= 0.;
2133 j++;
2134 }
2135 for (int i=0; i<np1-cutflag; i++) {
2136 cosa = std::cos(sthe+i*a);
2137 sina = std::sin(sthe+i*a);
2138 zz[j] = cz*cosa;
2139 rr[j] = cz*sina;
2140 j++;
2141 }
2142 if (j < np1)
2143 {
2144 zz[j]= zCut1;
2145 rr[j]= 0.;
2146 j++;
2147 }
2148 if (j > np1)
2149 {
2150 std::cerr << "Logic error in HepPolyhedronEllipsoid, memory corrupted!"
2151 << std::endl;
2152 }
2153 if (j < np1)
2154 {
2155 std::cerr << "Warning: logic error in HepPolyhedronEllipsoid."
2156 << std::endl;
2157 np1= j;
2158 }
2159 zz[j] = 0.;
2160 rr[j] = 0.;
2161
2162
2163 // R O T A T E P O L Y L I N E S
2164
2165 RotateAroundZ(0, 0.0, twopi, np1, 1, zz, rr, -1, 1);
2166 SetReferences();
2167
2168 delete [] zz;
2169 delete [] rr;
2170
2171 // rescale x and y vertex coordinates
2172 {
2173 Point3D<double> * p= pV;
2174 for (int i=0; i<nvert; i++, p++) {
2175 p->setX( p->x() * ax/cz );
2176 p->setY( p->y() * by/cz );
2177 }
2178 }
2179}
2180
2181HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
2182
2183HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(double ax,
2184 double ay,
2185 double h,
2186 double zTopCut)
2187/***********************************************************************
2188 * *
2189 * Name: HepPolyhedronEllipticalCone Date: 8.9.2005 *
2190 * Author: D.Anninos Revised: 9.9.2005 *
2191 * *
2192 * Function: Constructor for EllipticalCone *
2193 * *
2194 * Input: ax, ay - X & Y semi axes at z = 0 *
2195 * h - height of full cone *
2196 * zTopCut - Top Cut in Z Axis *
2197 * *
2198 ***********************************************************************/
2199{
2200 // C H E C K I N P U T P A R A M E T E R S
2201
2202 int k = 0;
2203 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2204
2205 if (k != 0) {
2206 std::cerr << "HepPolyhedronCone: error in input parameters";
2207 std::cerr << std::endl;
2208 return;
2209 }
2210
2211 // P R E P A R E T W O P O L Y L I N E S
2212
2213 zTopCut = (h >= zTopCut ? zTopCut : h);
2214
2215 double *zz, *rr;
2216 zz = new double[4];
2217 rr = new double[4];
2218 zz[0] = zTopCut;
2219 zz[1] = -zTopCut;
2220 zz[2] = zTopCut;
2221 zz[3] = -zTopCut;
2222 rr[0] = (h-zTopCut);
2223 rr[1] = (h+zTopCut);
2224 rr[2] = 0.;
2225 rr[3] = 0.;
2226
2227 // R O T A T E P O L Y L I N E S
2228
2229 RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, -1);
2230 SetReferences();
2231
2232 delete [] zz;
2233 delete [] rr;
2234
2235 // rescale x and y vertex coordinates
2236 {
2237 Point3D<double> * p= pV;
2238 for (int i=0; i<nvert; i++, p++) {
2239 p->setX( p->x() * ax );
2240 p->setY( p->y() * ay );
2241 }
2242 }
2243}
2244
2245HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
2246
2247int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
2248/***********************************************************************
2249 * *
2250 * Name: HepPolyhedron::fNumberOfRotationSteps Date: 24.06.97 *
2251 * Author: J.Allison (Manchester University) Revised: *
2252 * *
2253 * Function: Number of steps for whole circle *
2254 * *
2255 ***********************************************************************/
2256
2257#include "BooleanProcessor.src"
2258static BooleanProcessor processor;
2259
2260HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const
2261/***********************************************************************
2262 * *
2263 * Name: HepPolyhedron::add Date: 19.03.00 *
2264 * Author: E.Chernyaev Revised: *
2265 * *
2266 * Function: Boolean "union" of two polyhedra *
2267 * *
2268 ***********************************************************************/
2269{
2270 return processor.execute(OP_UNION, *this, p);
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 return processor.execute(OP_INTERSECTION, *this, p);
2284}
2285
2286HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const
2287/***********************************************************************
2288 * *
2289 * Name: HepPolyhedron::add Date: 19.03.00 *
2290 * Author: E.Chernyaev Revised: *
2291 * *
2292 * Function: Boolean "subtraction" of "p" from "this" *
2293 * *
2294 ***********************************************************************/
2295{
2296 return processor.execute(OP_SUBTRACTION, *this, p);
2297}
2298
2299bool HepPolyhedron::IsErrorBooleanProcess() const {
2300 return processor.get_processor_error();
2301}
Note: See TracBrowser for help on using the repository browser.