Context Navigation

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

Last change on this file since 850 was 850, checked in by garnier, 16 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
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	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(sqrtan2halfZhalfZ+k2);
1656
1657	for(int i = 1; i < n - 1; i++)
1658	{
1659	zz[i] = zz[i-1] - dz;
1660	rr[i] =std::sqrt(sqrtan2zz[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(sqrtan1halfZhalfZ+k1);
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	}
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
1685	HepPolyhedronHype::~HepPolyhedronHype() {}
1686
1687	HepPolyhedronCons::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
1763	HepPolyhedronCons::~HepPolyhedronCons() {}
1764
1765	HepPolyhedronCone::HepPolyhedronCone(double Rmn1, double Rmx1,
1766	double Rmn2, double Rmx2,
1767	double Dz) :
1768	HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0deg, 360deg) {}
1769
1770	HepPolyhedronCone::~HepPolyhedronCone() {}
1771
1772	HepPolyhedronTubs::HepPolyhedronTubs(double Rmin, double Rmax,
1773	double Dz,
1774	double Phi1, double Dphi)
1775	: HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1776
1777	HepPolyhedronTubs::~HepPolyhedronTubs() {}
1778
1779	HepPolyhedronTube::HepPolyhedronTube (double Rmin, double Rmax,
1780	double Dz)
1781	: HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0deg, 360deg) {}
1782
1783	HepPolyhedronTube::~HepPolyhedronTube () {}
1784
1785	HepPolyhedronPgon::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
1874	HepPolyhedronPgon::~HepPolyhedronPgon() {}
1875
1876	HepPolyhedronPcon::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
1882	HepPolyhedronPcon::~HepPolyhedronPcon() {}
1883
1884	HepPolyhedronSphere::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
1979	HepPolyhedronSphere::~HepPolyhedronSphere() {}
1980
1981	HepPolyhedronTorus::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
2054	HepPolyhedronTorus::~HepPolyhedronTorus() {}
2055
2056	HepPolyhedronEllipsoid::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
2181	HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() {}
2182
2183	HepPolyhedronEllipticalCone::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
2245	HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() {}
2246
2247	int 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"
2258	static BooleanProcessor processor;
2259
2260	HepPolyhedron 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
2273	HepPolyhedron 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
2286	HepPolyhedron 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
2299	bool HepPolyhedron::IsErrorBooleanProcess() const {
2300	return processor.get_processor_error();
2301	}

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

Download in other formats: