Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NURBStubesector.cc
Go to the documentation of this file.
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$
28 //
29 //
30 // Olivier Crumeyrolle 12 September 1996
31 
32 // Tubesector builder implementation
33 // OC 290896
34 
35 #include <sstream>
36 
37 #include "G4NURBStubesector.hh"
38 #include "G4PhysicalConstants.hh"
39 
41  G4double DZ, G4double PHI1, G4double PHI2)
42  : G4NURBS( 2, 3, // linear along U, quadratic along V
43  5, DecideNbrCtrlPts(PHI1, PHI2),
44  // rectangle along U, required stuff along V
45  // we must use a static function which
46  // take the two angles because the
47  // mother constructor is initialised
48  // before everything
49  Regular, // the knot vector along U will be generated
50  RegularRep ) // circular like knot vector also
51 {
52  // check angles
53  G4double deltaPHI = PHI2-PHI1;
54  while (deltaPHI <= 0) { PHI2 += twopi; deltaPHI += twopi; };
55 
56  G4int f = (int)floor(deltaPHI / (halfpi)); //number of pi/2 arcs
57 
58  const G4double mr = (r+R)/2;
59 
60  const G4double cp1 = std::cos(PHI1);
61  const G4double sp1 = std::sin(PHI1);
62  const G4double cp2 = std::cos(PHI2);
63  const G4double sp2 = std::sin(PHI2);
64 
65 
66  // define control points
67  CP(mpCtrlPts[ 0] , cp1*mr, sp1*mr, 0, 1 );
68  CP(mpCtrlPts[ 1] , cp1*mr, sp1*mr, 0, 1 );
69  CP(mpCtrlPts[ 2] , cp1*mr, sp1*mr, 0, 1 );
70  CP(mpCtrlPts[ 3] , cp1*mr, sp1*mr, 0, 1 );
71  CP(mpCtrlPts[ 4] , cp1*mr, sp1*mr, 0, 1 );
72 
73  CP(mpCtrlPts[ 5] , cp1*mr, sp1*mr, 0, 1 );
74  CP(mpCtrlPts[ 6] , cp1*mr, sp1*mr, 0, 1 );
75  CP(mpCtrlPts[ 7] , cp1*mr, sp1*mr, 0, 1 );
76  CP(mpCtrlPts[ 8] , cp1*mr, sp1*mr, 0, 1 );
77  CP(mpCtrlPts[ 9] , cp1*mr, sp1*mr, 0, 1 );
78 
79  CP(mpCtrlPts[10] , cp1*r, sp1*r, DZ, 1 );
80  CP(mpCtrlPts[11] , cp1*R, sp1*R, DZ, 1 );
81  CP(mpCtrlPts[12] , cp1*R, sp1*R, -DZ, 1 );
82  CP(mpCtrlPts[13] , cp1*r, sp1*r, -DZ, 1 );
83  CP(mpCtrlPts[14] , cp1*r, sp1*r, DZ, 1 );
84 
85  t_indCtrlPt i = 15;
86  G4double srcAngle = PHI1;
87  G4double deltaAngleo2;
88 
89  G4double destAngle = halfpi + PHI1;
90 
91  for(; f > 0; f--)
92  {
93  // the first arc CP is already Done
94 
95  deltaAngleo2 = (destAngle - srcAngle) / 2;
96  const G4double csa = std::cos(srcAngle);
97  const G4double ssa = std::sin(srcAngle);
98  const G4double tdao2 = std::tan(deltaAngleo2);
99 
100  // to calculate the intermediate CP :
101  // rotate by srcAngle the (1, tdao2) point
102  const t_Coord x = csa - ssa*tdao2;
103  const t_Coord y = ssa + csa*tdao2;
104 
105  // weight of the CP
106  const G4Float weight = (std::cos(deltaAngleo2));
107 
108  // initialization. postfix ++ because i initialized to 15
109  CP(mpCtrlPts[i++], x*r, y*r, DZ, 1, weight);
110  CP(mpCtrlPts[i++], x*R, y*R, DZ, 1, weight);
111  CP(mpCtrlPts[i++], x*R, y*R, -DZ, 1, weight);
112  CP(mpCtrlPts[i++], x*r, y*r, -DZ, 1, weight);
113  CP(mpCtrlPts[i++], x*r, y*r, DZ, 1, weight);
114 
115  // end CP (which is the first CP of the next arc)
116  const G4double cda = std::cos(destAngle);
117  const G4double sda = std::sin(destAngle);
118  CP(mpCtrlPts[i++], cda*r, sda*r, DZ, 1);
119  CP(mpCtrlPts[i++], cda*R, sda*R, DZ, 1);
120  CP(mpCtrlPts[i++], cda*R, sda*R, -DZ, 1);
121  CP(mpCtrlPts[i++], cda*r, sda*r, -DZ, 1);
122  CP(mpCtrlPts[i++], cda*r, sda*r, DZ, 1);
123 
124  // prepare next arc
125  srcAngle = destAngle;
126  destAngle += halfpi;
127  }
128 
129  // f == 0, final Arc
130  // could be handled in the loops
131 
132  destAngle = PHI2;
133  deltaAngleo2 = (destAngle - srcAngle) / 2;
134  const G4double csa = std::cos(srcAngle);
135  const G4double ssa = std::sin(srcAngle);
136  const G4double tdao2 = std::tan(deltaAngleo2);
137 
138  // to calculate the intermediate CP :
139  // rotate by srcAngle the (1, tdao2) point
140  const t_Coord x = csa - ssa*tdao2;
141  const t_Coord y = ssa + csa*tdao2;
142 
143  // weight of the CP
144  const G4Float weight = (std::cos(deltaAngleo2));
145 
146  // initialization.
147  CP(mpCtrlPts[i++], x*r, y*r, DZ, 1, weight);
148  CP(mpCtrlPts[i++], x*R, y*R, DZ, 1, weight);
149  CP(mpCtrlPts[i++], x*R, y*R, -DZ, 1, weight);
150  CP(mpCtrlPts[i++], x*r, y*r, -DZ, 1, weight);
151  CP(mpCtrlPts[i++], x*r, y*r, DZ, 1, weight);
152 
153  // end CP
154  const G4double cda = std::cos(destAngle);
155  const G4double sda = std::sin(destAngle);
156  CP(mpCtrlPts[i++], cda*r, sda*r, DZ, 1);
157  CP(mpCtrlPts[i++], cda*R, sda*R, DZ, 1);
158  CP(mpCtrlPts[i++], cda*R, sda*R, -DZ, 1);
159  CP(mpCtrlPts[i++], cda*r, sda*r, -DZ, 1);
160  CP(mpCtrlPts[i++], cda*r, sda*r, DZ, 1);
161 
162  if (i != (mtotnbrCtrlPts - 10) )
163  {
164  G4cerr << "\nERROR: G4NURBStubesector::G4NURBStubesector: wrong index,"
165  << i << " instead of " << (mtotnbrCtrlPts - 10)
166  << "\n\tThe tubesector won't be correct."
167  << G4endl;
168  }
169 
170  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
171  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
172  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
173  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
174  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
175 
176  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
177  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
178  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
179  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
180  CP(mpCtrlPts[i++] , cp2*mr, sp2*mr, 0, 1);
181 
182  // possible to put a DZ DZ -DZ -DZ DZ column to scratch
183  // to a line instead of a point
184 
185  // creating the nurbs identity
186  std::ostringstream tmpstr;
187  tmpstr << "Tubs" << " \tPHI1=" << PHI1 << " ; PHI2=" << PHI2;
188  mpwhoami = new char [tmpstr.str().length() + 1];
189  mpwhoami = std::strcpy(mpwhoami, tmpstr.str().c_str());
190 }
191 
192 const char* G4NURBStubesector::Whoami() const
193 {
194  return mpwhoami;
195 }
196 
198 {
199  if (mpwhoami) { delete [] mpwhoami; mpwhoami = 0; }
200 }
201 
203 G4NURBStubesector::DecideNbrCtrlPts(G4double PHI1, G4double PHI2)
204 {
205  // check angles
206  G4double deltaPHI = PHI2-PHI1;
207  while (deltaPHI <= 0) { PHI2 += twopi; deltaPHI += twopi; }
208  G4double k = deltaPHI / (halfpi);
209 
210  // G4cerr << " k " << k << G4endl;
211  // G4cerr << " fk " << std::floor(k) << G4endl;
212  // G4cerr << " ifk " << ((int)(std::floor(k))) << G4endl;
213  // G4cerr << " n " << (2*((int)(std::floor(k))) + 7) << G4endl;
214 
215  return ( 2*((int)(std::floor(k))) + 7 );
216 }