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