Context Navigation

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

Last change on this file since 1058 was 1058, checked in by garnier, 15 years ago
file release 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.32 2008/11/13 09:05:27 gcosmo Exp $
28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
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
71	using namespace HepGeom;
72	using CLHEP::perMillion;
73	using CLHEP::deg;
74	using CLHEP::pi;
75	using 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	***********************************************************************/
85	std::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
92	std::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
107	HepPolyhedron::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
121	HepPolyhedron & 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
139	int
140	HepPolyhedron::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
168	Normal3D<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
196	int 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
209	void 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
231	void 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
244	void 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
270	void 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
290	void 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, v2i2,0, (i2+1),0, v1(i1+1),0);
328	}
329	}else{
330	if (r1 == 0.) {
331	pF[kface++] = G4Facet(vvi1,0, v2i2,0, vEdge*(i2+1),0);
332	for (i2++,i=1; i<ns-1; i2++,i++) {
333	pF[kface++] = G4Facet(vEdgei1,0, v2i2,0, vEdge*(i2+1),0);
334	}
335	pF[kface++] = G4Facet(vEdgei1,0, v2i2,0, vv*ii2,0);
336	}else if (r2 == 0.) {
337	pF[kface++] = G4Facet(vvi1,0, vEdgei2,0, v1*(i1+1),0);
338	for (i1++,i=1; i<ns-1; i1++,i++) {
339	pF[kface++] = G4Facet(vEdgei1,0, vEdgei2,0, v1*(i1+1),0);
340	}
341	pF[kface++] = G4Facet(vEdgei1,0, vvi2,0, v1*ii1,0);
342	}else{
343	pF[kface++] = G4Facet(vvi1,0, v2i2,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(vEdgei1,0, v2i2,0, vEdge(i2+1),0,v1(i1+1),0);
346	}
347	pF[kface++] = G4Facet(vEdgei1,0, v2i2,0, vvii2,0, v1ii1,0);
348	}
349	}
350	}
351
352	void 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
423	void 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 ? absNp1nSphi : (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 ? 2absNp1 : 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
670	void 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
807	void 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
833	HepPolyhedron & 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
858	bool 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
885	Point3D<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
904	bool
905	HepPolyhedron::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
923	bool 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
955	bool 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 (iOrderk1 > iOrderk2);
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
1007	bool
1008	HepPolyhedron::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
1023	bool
1024	HepPolyhedron::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
1044	bool
1045	HepPolyhedron::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
1065	void 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
1099	void 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
1120	bool
1121	HepPolyhedron::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
1151	Normal3D<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
1176	Normal3D<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
1201	bool 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
1222	bool 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
1238	double 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
1260	double 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
1288	int
1289	HepPolyhedron::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
1357	int
1358	HepPolyhedron::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
1388	HepPolyhedronTrd2::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
1420	HepPolyhedronTrd2::~HepPolyhedronTrd2() {}
1421
1422	HepPolyhedronTrd1::HepPolyhedronTrd1(double Dx1, double Dx2,
1423	double Dy, double Dz)
1424	: HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1425
1426	HepPolyhedronTrd1::~HepPolyhedronTrd1() {}
1427
1428	HepPolyhedronBox::HepPolyhedronBox(double Dx, double Dy, double Dz)
1429	: HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1430
1431	HepPolyhedronBox::~HepPolyhedronBox() {}
1432
1433	HepPolyhedronTrap::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 = Dzstd::tan(Theta)std::cos(Phi);
1468	double DzTthetaSphi = Dzstd::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
1486	HepPolyhedronTrap::~HepPolyhedronTrap() {}
1487
1488	HepPolyhedronPara::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
1493	HepPolyhedronPara::~HepPolyhedronPara() {}
1494
1495	HepPolyhedronParaboloid::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 = (r2r2 - r1r1) / 2 / dz;
1561	double k2 = (r2r2 + r1r1) / 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
1597	HepPolyhedronParaboloid::~HepPolyhedronParaboloid() {}
1598
1599	HepPolyhedronHype::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(sqrtan2halfZhalfZ+k2);
1657
1658	for(int i = 1; i < n-1; i++)
1659	{
1660	zz[i] = zz[i-1] - dz;
1661	rr[i] =std::sqrt(sqrtan2zz[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(sqrtan1halfZhalfZ+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(sqrtan1zz[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
1684	HepPolyhedronHype::~HepPolyhedronHype() {}
1685
1686	HepPolyhedronCons::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
1762	HepPolyhedronCons::~HepPolyhedronCons() {}
1763
1764	HepPolyhedronCone::HepPolyhedronCone(double Rmn1, double Rmx1,
1765	double Rmn2, double Rmx2,
1766	double Dz) :
1767	HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0deg, 360deg) {}
1768
1769	HepPolyhedronCone::~HepPolyhedronCone() {}
1770
1771	HepPolyhedronTubs::HepPolyhedronTubs(double Rmin, double Rmax,
1772	double Dz,
1773	double Phi1, double Dphi)
1774	: HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1775
1776	HepPolyhedronTubs::~HepPolyhedronTubs() {}
1777
1778	HepPolyhedronTube::HepPolyhedronTube (double Rmin, double Rmax,
1779	double Dz)
1780	: HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0deg, 360deg) {}
1781
1782	HepPolyhedronTube::~HepPolyhedronTube () {}
1783
1784	HepPolyhedronPgon::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
1873	HepPolyhedronPgon::~HepPolyhedronPgon() {}
1874
1875	HepPolyhedronPcon::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
1881	HepPolyhedronPcon::~HepPolyhedronPcon() {}
1882
1883	HepPolyhedronSphere::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
1978	HepPolyhedronSphere::~HepPolyhedronSphere() {}
1979
1980	HepPolyhedronTorus::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
2053	HepPolyhedronTorus::~HepPolyhedronTorus() {}
2054
2055	HepPolyhedronEllipsoid::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
2180	HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
2181
2182	HepPolyhedronEllipticalCone::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
2244	HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
2245
2246	int 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	static BooleanProcessor processor;
2258
2259	HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const
2260	/***********************************************************************
2261	* *
2262	* Name: HepPolyhedron::add Date: 19.03.00 *
2263	* Author: E.Chernyaev Revised: *
2264	* *
2265	* Function: Boolean "union" of two polyhedra *
2266	* *
2267	***********************************************************************/
2268	{
2269	return processor.execute(OP_UNION, *this, p);
2270	}
2271
2272	HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const
2273	/***********************************************************************
2274	* *
2275	* Name: HepPolyhedron::intersect Date: 19.03.00 *
2276	* Author: E.Chernyaev Revised: *
2277	* *
2278	* Function: Boolean "intersection" of two polyhedra *
2279	* *
2280	***********************************************************************/
2281	{
2282	return processor.execute(OP_INTERSECTION, *this, p);
2283	}
2284
2285	HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const
2286	/***********************************************************************
2287	* *
2288	* Name: HepPolyhedron::add Date: 19.03.00 *
2289	* Author: E.Chernyaev Revised: *
2290	* *
2291	* Function: Boolean "subtraction" of "p" from "this" *
2292	* *
2293	***********************************************************************/
2294	{
2295	return processor.execute(OP_SUBTRACTION, *this, p);
2296	}
2297
2298	bool HepPolyhedron::IsErrorBooleanProcess() const {
2299	return processor.get_processor_error();
2300	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: