| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4TwistBoxSide implementation 26 // G4TwistBoxSide implementation
27 // 27 //
28 // Author: 18/03/2005 - O.Link (Oliver.Link@ce 28 // Author: 18/03/2005 - O.Link (Oliver.Link@cern.ch)
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include <cmath> 31 #include <cmath>
32 32
33 #include "G4TwistBoxSide.hh" 33 #include "G4TwistBoxSide.hh"
34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh"
35 #include "G4JTPolynomialSolver.hh" 35 #include "G4JTPolynomialSolver.hh"
36 36
37 //============================================ 37 //=====================================================================
38 //* constructors ----------------------------- 38 //* constructors ------------------------------------------------------
39 39
40 G4TwistBoxSide::G4TwistBoxSide(const G4String& 40 G4TwistBoxSide::G4TwistBoxSide(const G4String& name,
41 G4double PhiTw 41 G4double PhiTwist, // twist angle
42 G4double pDz, 42 G4double pDz, // half z lenght
43 G4double pThet 43 G4double pTheta, // direction between end planes
44 G4double pPhi, 44 G4double pPhi, // defined by polar and azimutal angles.
45 G4double pDy1, 45 G4double pDy1, // half y length at -pDz
46 G4double pDx1, 46 G4double pDx1, // half x length at -pDz,-pDy
47 G4double pDx2, 47 G4double pDx2, // half x length at -pDz,+pDy
48 G4double pDy2, 48 G4double pDy2, // half y length at +pDz
49 G4double pDx3, 49 G4double pDx3, // half x length at +pDz,-pDy
50 G4double pDx4, 50 G4double pDx4, // half x length at +pDz,+pDy
51 G4double pAlph 51 G4double pAlph, // tilt angle at +pDz
52 G4double Angle 52 G4double AngleSide // parity
53 53 ) : G4VTwistSurface(name)
54 { 54 {
55 55
56 56
57 fAxis[0] = kYAxis; // in local coordinate 57 fAxis[0] = kYAxis; // in local coordinate system
58 fAxis[1] = kZAxis; 58 fAxis[1] = kZAxis;
59 fAxisMin[0] = -kInfinity ; // Y Axis bounda 59 fAxisMin[0] = -kInfinity ; // Y Axis boundary
60 fAxisMax[0] = kInfinity ; // depends on 60 fAxisMax[0] = kInfinity ; // depends on z !!
61 fAxisMin[1] = -pDz ; // Z Axis boundary 61 fAxisMin[1] = -pDz ; // Z Axis boundary
62 fAxisMax[1] = pDz ; 62 fAxisMax[1] = pDz ;
63 63
64 fDx1 = pDx1 ; 64 fDx1 = pDx1 ;
65 fDx2 = pDx2 ; // box 65 fDx2 = pDx2 ; // box
66 fDx3 = pDx3 ; 66 fDx3 = pDx3 ;
67 fDx4 = pDx4 ; // box 67 fDx4 = pDx4 ; // box
68 68
69 // this is an overhead. But the parameter na 69 // this is an overhead. But the parameter naming scheme fits to the other surfaces.
70 70
71 if ( fDx1 != fDx2 || fDx3 != fDx4 ) << 71 if ( ! (fDx1 == fDx2 && fDx3 == fDx4 ) )
72 { 72 {
73 std::ostringstream message; 73 std::ostringstream message;
74 message << "TwistedTrapBoxSide is not used 74 message << "TwistedTrapBoxSide is not used as a the side of a box: "
75 << GetName() << G4endl 75 << GetName() << G4endl
76 << " Not a box !"; 76 << " Not a box !";
77 G4Exception("G4TwistBoxSide::G4TwistBoxSid 77 G4Exception("G4TwistBoxSide::G4TwistBoxSide()", "GeomSolids0002",
78 FatalException, message); 78 FatalException, message);
79 } 79 }
80 80
81 fDy1 = pDy1 ; 81 fDy1 = pDy1 ;
82 fDy2 = pDy2 ; 82 fDy2 = pDy2 ;
83 83
84 fDz = pDz ; 84 fDz = pDz ;
85 85
86 fAlph = pAlph ; 86 fAlph = pAlph ;
87 fTAlph = std::tan(fAlph) ; 87 fTAlph = std::tan(fAlph) ;
88 88
89 fTheta = pTheta ; 89 fTheta = pTheta ;
90 fPhi = pPhi ; 90 fPhi = pPhi ;
91 91
92 // precalculate frequently used parameters 92 // precalculate frequently used parameters
93 93
94 fDx4plus2 = fDx4 + fDx2 ; 94 fDx4plus2 = fDx4 + fDx2 ;
95 fDx4minus2 = fDx4 - fDx2 ; 95 fDx4minus2 = fDx4 - fDx2 ;
96 fDx3plus1 = fDx3 + fDx1 ; 96 fDx3plus1 = fDx3 + fDx1 ;
97 fDx3minus1 = fDx3 - fDx1 ; 97 fDx3minus1 = fDx3 - fDx1 ;
98 fDy2plus1 = fDy2 + fDy1 ; 98 fDy2plus1 = fDy2 + fDy1 ;
99 fDy2minus1 = fDy2 - fDy1 ; 99 fDy2minus1 = fDy2 - fDy1 ;
100 100
101 fa1md1 = 2*fDx2 - 2*fDx1 ; 101 fa1md1 = 2*fDx2 - 2*fDx1 ;
102 fa2md2 = 2*fDx4 - 2*fDx3 ; 102 fa2md2 = 2*fDx4 - 2*fDx3 ;
103 103
104 104
105 fPhiTwist = PhiTwist ; // dphi 105 fPhiTwist = PhiTwist ; // dphi
106 fAngleSide = AngleSide ; // 0,90,180,270 de 106 fAngleSide = AngleSide ; // 0,90,180,270 deg
107 107
108 fdeltaX = 2*fDz * std::tan(fTheta) * std::co 108 fdeltaX = 2*fDz * std::tan(fTheta) * std::cos(fPhi); // dx in surface equation
109 fdeltaY = 2*fDz * std::tan(fTheta) * std::si 109 fdeltaY = 2*fDz * std::tan(fTheta) * std::sin(fPhi); // dy in surface equation
110 110
111 fRot.rotateZ( AngleSide ) ; 111 fRot.rotateZ( AngleSide ) ;
112 112
113 fTrans.set(0, 0, 0); // No Translation 113 fTrans.set(0, 0, 0); // No Translation
114 fIsValidNorm = false; 114 fIsValidNorm = false;
115 115
116 SetCorners(); 116 SetCorners();
117 SetBoundaries(); 117 SetBoundaries();
118 } 118 }
119 119
120 120
121 //============================================ 121 //=====================================================================
122 //* Fake default constructor ----------------- 122 //* Fake default constructor ------------------------------------------
123 123
124 G4TwistBoxSide::G4TwistBoxSide( __void__& a ) 124 G4TwistBoxSide::G4TwistBoxSide( __void__& a )
125 : G4VTwistSurface(a) << 125 : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
>> 126 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
>> 127 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.),
>> 128 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.),
>> 129 fa2md2(0.)
126 { 130 {
127 } 131 }
128 132
129 133
130 //============================================ 134 //=====================================================================
131 //* destructor ------------------------------- 135 //* destructor --------------------------------------------------------
132 136
133 G4TwistBoxSide::~G4TwistBoxSide() = default; << 137 G4TwistBoxSide::~G4TwistBoxSide()
>> 138 {
>> 139 }
134 140
135 //============================================ 141 //=====================================================================
136 //* GetNormal -------------------------------- 142 //* GetNormal ---------------------------------------------------------
137 143
138 G4ThreeVector G4TwistBoxSide::GetNormal(const 144 G4ThreeVector G4TwistBoxSide::GetNormal(const G4ThreeVector& tmpxx,
139 145 G4bool isGlobal)
140 { 146 {
141 // GetNormal returns a normal vector at a s 147 // GetNormal returns a normal vector at a surface (or very close
142 // to surface) point at tmpxx. 148 // to surface) point at tmpxx.
143 // If isGlobal=true, it returns the normal 149 // If isGlobal=true, it returns the normal in global coordinate.
144 // 150 //
145 151
146 G4ThreeVector xx; 152 G4ThreeVector xx;
147 if (isGlobal) 153 if (isGlobal)
148 { 154 {
149 xx = ComputeLocalPoint(tmpxx); 155 xx = ComputeLocalPoint(tmpxx);
150 if ((xx - fCurrentNormal.p).mag() < 0.5 156 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
151 { 157 {
152 return ComputeGlobalDirection(fCurren 158 return ComputeGlobalDirection(fCurrentNormal.normal);
153 } 159 }
154 } 160 }
155 else 161 else
156 { 162 {
157 xx = tmpxx; 163 xx = tmpxx;
158 if (xx == fCurrentNormal.p) 164 if (xx == fCurrentNormal.p)
159 { 165 {
160 return fCurrentNormal.normal; 166 return fCurrentNormal.normal;
161 } 167 }
162 } 168 }
163 169
164 G4double phi ; 170 G4double phi ;
165 G4double u ; 171 G4double u ;
166 172
167 GetPhiUAtX(xx,phi,u) ; // phi,u for point 173 GetPhiUAtX(xx,phi,u) ; // phi,u for point xx close to surface
168 174
169 G4ThreeVector normal = NormAng(phi,u) ; // 175 G4ThreeVector normal = NormAng(phi,u) ; // the normal vector at phi,u
170 176
171 #ifdef G4TWISTDEBUG 177 #ifdef G4TWISTDEBUG
172 G4cout << "normal vector = " << normal << 178 G4cout << "normal vector = " << normal << G4endl ;
173 G4cout << "phi = " << phi << " , u = " << u 179 G4cout << "phi = " << phi << " , u = " << u << G4endl ;
174 #endif 180 #endif
175 181
176 // normal = normal/normal.mag() ; 182 // normal = normal/normal.mag() ;
177 183
178 if (isGlobal) 184 if (isGlobal)
179 { 185 {
180 fCurrentNormal.normal = ComputeGlobalDir 186 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
181 } 187 }
182 else 188 else
183 { 189 {
184 fCurrentNormal.normal = normal.unit(); 190 fCurrentNormal.normal = normal.unit();
185 } 191 }
186 return fCurrentNormal.normal; 192 return fCurrentNormal.normal;
187 } 193 }
188 194
189 //============================================ 195 //=====================================================================
190 //* DistanceToSurface ------------------------ 196 //* DistanceToSurface -------------------------------------------------
191 197
192 G4int G4TwistBoxSide::DistanceToSurface(const 198 G4int G4TwistBoxSide::DistanceToSurface(const G4ThreeVector& gp,
193 const 199 const G4ThreeVector& gv,
194 200 G4ThreeVector gxx[],
195 201 G4double distance[],
196 202 G4int areacode[],
197 203 G4bool isvalid[],
198 204 EValidate validate)
199 { 205 {
200 206
201 static const G4double pihalf = pi/2 ; 207 static const G4double pihalf = pi/2 ;
202 const G4double ctol = 0.5 * kCarTolerance; 208 const G4double ctol = 0.5 * kCarTolerance;
203 209
204 G4bool IsParallel = false ; 210 G4bool IsParallel = false ;
205 G4bool IsConverged = false ; 211 G4bool IsConverged = false ;
206 212
207 G4int nxx = 0 ; // number of physical solut 213 G4int nxx = 0 ; // number of physical solutions
208 214
209 fCurStatWithV.ResetfDone(validate, &gp, &gv) 215 fCurStatWithV.ResetfDone(validate, &gp, &gv);
210 216
211 if (fCurStatWithV.IsDone()) 217 if (fCurStatWithV.IsDone())
212 { 218 {
213 for (G4int i=0; i<fCurStatWithV.GetNXX(); 219 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
214 { 220 {
215 gxx[i] = fCurStatWithV.GetXX(i); 221 gxx[i] = fCurStatWithV.GetXX(i);
216 distance[i] = fCurStatWithV.GetDistance( 222 distance[i] = fCurStatWithV.GetDistance(i);
217 areacode[i] = fCurStatWithV.GetAreacode( 223 areacode[i] = fCurStatWithV.GetAreacode(i);
218 isvalid[i] = fCurStatWithV.IsValid(i); 224 isvalid[i] = fCurStatWithV.IsValid(i);
219 } 225 }
220 return fCurStatWithV.GetNXX(); 226 return fCurStatWithV.GetNXX();
221 } 227 }
222 else // initialize 228 else // initialize
223 { 229 {
224 for (G4int i=0; i<G4VSURFACENXX ; ++i) 230 for (G4int i=0; i<G4VSURFACENXX ; ++i)
225 { 231 {
226 distance[i] = kInfinity; 232 distance[i] = kInfinity;
227 areacode[i] = sOutside; 233 areacode[i] = sOutside;
228 isvalid[i] = false; 234 isvalid[i] = false;
229 gxx[i].set(kInfinity, kInfinity, kInfini 235 gxx[i].set(kInfinity, kInfinity, kInfinity);
230 } 236 }
231 } 237 }
232 238
233 G4ThreeVector p = ComputeLocalPoint(gp); 239 G4ThreeVector p = ComputeLocalPoint(gp);
234 G4ThreeVector v = ComputeLocalDirection(gv); 240 G4ThreeVector v = ComputeLocalDirection(gv);
235 241
236 #ifdef G4TWISTDEBUG 242 #ifdef G4TWISTDEBUG
237 G4cout << "Local point p = " << p << G4endl 243 G4cout << "Local point p = " << p << G4endl ;
238 G4cout << "Local direction v = " << v << G4e 244 G4cout << "Local direction v = " << v << G4endl ;
239 #endif 245 #endif
240 246
241 G4double phi=0.,u=0.,q=0.; // parameters 247 G4double phi=0.,u=0.,q=0.; // parameters
242 248
243 // temporary variables 249 // temporary variables
244 250
245 G4double tmpdist = kInfinity ; 251 G4double tmpdist = kInfinity ;
246 G4ThreeVector tmpxx; 252 G4ThreeVector tmpxx;
247 G4int tmpareacode = sOutside ; 253 G4int tmpareacode = sOutside ;
248 G4bool tmpisvalid = false ; 254 G4bool tmpisvalid = false ;
249 255
250 std::vector<Intersection> xbuf ; 256 std::vector<Intersection> xbuf ;
251 Intersection xbuftmp ; 257 Intersection xbuftmp ;
252 258
253 // prepare some variables for the intersecti 259 // prepare some variables for the intersection finder
254 260
255 G4double L = 2*fDz ; 261 G4double L = 2*fDz ;
256 262
257 G4double phixz = fPhiTwist * ( p.x() * v.z() 263 G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ;
258 G4double phiyz = fPhiTwist * ( p.y() * v.z() 264 G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ;
259 265
260 // special case vz = 0 266 // special case vz = 0
261 267
262 if ( v.z() == 0. ) 268 if ( v.z() == 0. )
263 { 269 {
264 if ( (std::fabs(p.z()) <= L) && (fPhiTwist << 270 if ( (std::fabs(p.z()) <= L) && fPhiTwist ) // intersection possible in z
265 { 271 {
266 phi = p.z() * fPhiTwist / L ; // phi is 272 phi = p.z() * fPhiTwist / L ; // phi is determined by the z-position
267 273
268 q = (2.* fPhiTwist*((v.x() - fTAlph*v.y( 274 q = (2.* fPhiTwist*((v.x() - fTAlph*v.y())*std::cos(phi)
269 + (fTAlph*v.x() + 275 + (fTAlph*v.x() + v.y())*std::sin(phi)));
270 276
271 if (q != 0.0) << 277 if (q)
272 { 278 {
273 u = (2*(-(fdeltaY*phi*v.x()) + fPhiTwi 279 u = (2*(-(fdeltaY*phi*v.x()) + fPhiTwist*p.y()*v.x()
274 + fdeltaX*phi*v.y() - fPhiTwis 280 + fdeltaX*phi*v.y() - fPhiTwist*p.x()*v.y())
275 + (fDx4plus2*fPhiTwist + 2*fDx4mi 281 + (fDx4plus2*fPhiTwist + 2*fDx4minus2*phi)
276 * (v.y()*std::cos(phi) - v.x()*st 282 * (v.y()*std::cos(phi) - v.x()*std::sin(phi))) / q;
277 } 283 }
278 xbuftmp.phi = phi ; 284 xbuftmp.phi = phi ;
279 xbuftmp.u = u ; 285 xbuftmp.u = u ;
280 xbuftmp.areacode = sOutside ; 286 xbuftmp.areacode = sOutside ;
281 xbuftmp.distance = kInfinity ; 287 xbuftmp.distance = kInfinity ;
282 xbuftmp.isvalid = false ; 288 xbuftmp.isvalid = false ;
283 289
284 xbuf.push_back(xbuftmp) ; // store it t 290 xbuf.push_back(xbuftmp) ; // store it to xbuf
285 } 291 }
286 else // no interse 292 else // no intersection possible
287 { 293 {
288 distance[0] = kInfinity; 294 distance[0] = kInfinity;
289 gxx[0].set(kInfinity,kInfinity,kInfinity 295 gxx[0].set(kInfinity,kInfinity,kInfinity);
290 isvalid[0] = false ; 296 isvalid[0] = false ;
291 areacode[0] = sOutside ; 297 areacode[0] = sOutside ;
292 fCurStatWithV.SetCurrentStatus(0, gxx[0] 298 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
293 areacode[ 299 areacode[0], isvalid[0],
294 0, valida 300 0, validate, &gp, &gv);
295 301
296 return 0; 302 return 0;
297 } // end std::fabs(p.z() <= L 303 } // end std::fabs(p.z() <= L
298 } // end v.z() == 0 304 } // end v.z() == 0
299 else // general solution for non-zero vz 305 else // general solution for non-zero vz
300 { 306 {
301 G4double c[8],srd[7],si[7] ; 307 G4double c[8],srd[7],si[7] ;
302 308
303 c[7] = -14400*(-2*phixz + 2*fTAlph*phiyz + 309 c[7] = -14400*(-2*phixz + 2*fTAlph*phiyz + fDx4plus2*fPhiTwist*v.z()) ;
304 c[6] = 28800*(phiyz + 2*fDz*v.x() - (fdelt 310 c[6] = 28800*(phiyz + 2*fDz*v.x() - (fdeltaX + fDx4minus2)*v.z() + fTAlph*(phixz - 2*fDz*v.y() + fdeltaY*v.z())) ;
305 c[5] = -1200*(10*phixz - 48*fDz*v.y() + 24 311 c[5] = -1200*(10*phixz - 48*fDz*v.y() + 24*fdeltaY*v.z() + fDx4plus2*fPhiTwist*v.z() - 2*fTAlph*(5*phiyz + 24*fDz*v.x() - 12*fdeltaX*v.z())) ;
306 c[4] = -2400*(phiyz + 10*fDz*v.x() - 5*fde 312 c[4] = -2400*(phiyz + 10*fDz*v.x() - 5*fdeltaX*v.z() + fDx4minus2*v.z() + fTAlph*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())) ;
307 c[3] = 24*(2*phixz - 200*fDz*v.y() + 100*f 313 c[3] = 24*(2*phixz - 200*fDz*v.y() + 100*fdeltaY*v.z() - fDx4plus2*fPhiTwist*v.z() - 2*fTAlph*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z())) ;
308 c[2] = -16*(7*fTAlph* phixz + 7*phiyz - 6* 314 c[2] = -16*(7*fTAlph* phixz + 7*phiyz - 6*fDz*v.x() + 6*fDz*fTAlph*v.y() + 3*(fdeltaX + fDx4minus2 - fdeltaY*fTAlph)*v.z()) ;
309 c[1] = 4*(9*phixz - 9*fTAlph*phiyz - 56*fD 315 c[1] = 4*(9*phixz - 9*fTAlph*phiyz - 56*fDz*fTAlph*v.x() - 56*fDz*v.y() + 28*(fdeltaY + fdeltaX*fTAlph)*v.z()) ;
310 c[0] = 36*(2* fDz*(v.x() - fTAlph*v.y()) - 316 c[0] = 36*(2* fDz*(v.x() - fTAlph*v.y()) - fdeltaX*v.z() + fdeltaY*fTAlph*v.z()) ;
311 317
312 #ifdef G4TWISTDEBUG 318 #ifdef G4TWISTDEBUG
313 G4cout << "coef = " << c[0] << " " 319 G4cout << "coef = " << c[0] << " "
314 << c[1] << " " 320 << c[1] << " "
315 << c[2] << " " 321 << c[2] << " "
316 << c[3] << " " 322 << c[3] << " "
317 << c[4] << " " 323 << c[4] << " "
318 << c[5] << " " 324 << c[5] << " "
319 << c[6] << " " 325 << c[6] << " "
320 << c[7] << G4endl ; 326 << c[7] << G4endl ;
321 #endif 327 #endif
322 328
323 G4JTPolynomialSolver trapEq ; 329 G4JTPolynomialSolver trapEq ;
324 G4int num = trapEq.FindRoots(c,7,srd,si); 330 G4int num = trapEq.FindRoots(c,7,srd,si);
325 331
326 for (G4int i = 0 ; i<num ; ++i ) // loo 332 for (G4int i = 0 ; i<num ; ++i ) // loop over all mathematical solutions
327 { 333 {
328 if ( (si[i]==0.0) && (fPhiTwist != 0.0) << 334 if ( (si[i]==0.0) && fPhiTwist ) // only real solutions
329 { 335 {
330 #ifdef G4TWISTDEBUG 336 #ifdef G4TWISTDEBUG
331 G4cout << "Solution " << i << " : " << 337 G4cout << "Solution " << i << " : " << srd[i] << G4endl ;
332 #endif 338 #endif
333 phi = std::fmod(srd[i], pihalf) ; 339 phi = std::fmod(srd[i], pihalf) ;
334 340
335 u = (2*phiyz + 4*fDz*phi*v.y() - 2*f 341 u = (2*phiyz + 4*fDz*phi*v.y() - 2*fdeltaY*phi*v.z()
336 - fDx4plus2*fPhiTwist*v.z()*std:: 342 - fDx4plus2*fPhiTwist*v.z()*std::sin(phi)
337 - 2*fDx4minus2*phi*v.z()*std::sin 343 - 2*fDx4minus2*phi*v.z()*std::sin(phi))
338 /(2*fPhiTwist*v.z()*std::cos(phi) 344 /(2*fPhiTwist*v.z()*std::cos(phi)
339 + 2*fPhiTwist*fTAlph*v.z()*std: 345 + 2*fPhiTwist*fTAlph*v.z()*std::sin(phi)) ;
340 346
341 xbuftmp.phi = phi ; 347 xbuftmp.phi = phi ;
342 xbuftmp.u = u ; 348 xbuftmp.u = u ;
343 xbuftmp.areacode = sOutside ; 349 xbuftmp.areacode = sOutside ;
344 xbuftmp.distance = kInfinity ; 350 xbuftmp.distance = kInfinity ;
345 xbuftmp.isvalid = false ; 351 xbuftmp.isvalid = false ;
346 352
347 xbuf.push_back(xbuftmp) ; // store it 353 xbuf.push_back(xbuftmp) ; // store it to xbuf
348 354
349 #ifdef G4TWISTDEBUG 355 #ifdef G4TWISTDEBUG
350 G4cout << "solution " << i << " = " << 356 G4cout << "solution " << i << " = " << phi << " , " << u << G4endl ;
351 #endif 357 #endif
352 } // end if real solution 358 } // end if real solution
353 } // end loop i 359 } // end loop i
354 } // end general case 360 } // end general case
355 361
356 362
357 nxx = (G4int)xbuf.size() ; // save the numb << 363 nxx = xbuf.size() ; // save the number of solutions
358 364
359 G4ThreeVector xxonsurface ; // point 365 G4ThreeVector xxonsurface ; // point on surface
360 G4ThreeVector surfacenormal ; // normal 366 G4ThreeVector surfacenormal ; // normal vector
361 G4double deltaX ; // distan 367 G4double deltaX ; // distance between intersection point and
362 // point 368 // point on surface
363 G4double theta ; // angle 369 G4double theta ; // angle between track and surfacenormal
364 G4double factor ; // a scal 370 G4double factor ; // a scaling factor
365 G4int maxint = 30 ; // number 371 G4int maxint = 30 ; // number of iterations
366 372
367 373
368 for (auto & k : xbuf) << 374 for ( size_t k = 0 ; k<xbuf.size() ; ++k )
369 { 375 {
370 #ifdef G4TWISTDEBUG 376 #ifdef G4TWISTDEBUG
371 G4cout << "Solution " << k << " : " 377 G4cout << "Solution " << k << " : "
372 << "reconstructed phiR = " << xbuf[ 378 << "reconstructed phiR = " << xbuf[k].phi
373 << ", uR = " << xbuf[k].u << G4endl 379 << ", uR = " << xbuf[k].u << G4endl ;
374 #endif 380 #endif
375 381
376 phi = k.phi ; // get the stored values fo << 382 phi = xbuf[k].phi ; // get the stored values for phi and u
377 u = k.u ; << 383 u = xbuf[k].u ;
378 384
379 IsConverged = false ; // no convergence 385 IsConverged = false ; // no convergence at the beginning
380 386
381 for ( G4int i = 1 ; i<maxint ; ++i ) 387 for ( G4int i = 1 ; i<maxint ; ++i )
382 { 388 {
383 xxonsurface = SurfacePoint(phi,u) ; 389 xxonsurface = SurfacePoint(phi,u) ;
384 surfacenormal = NormAng(phi,u) ; 390 surfacenormal = NormAng(phi,u) ;
385 tmpdist = DistanceToPlaneWithV(p, v, xxo 391 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx);
386 deltaX = ( tmpxx - xxonsurface ).mag() ; 392 deltaX = ( tmpxx - xxonsurface ).mag() ;
387 theta = std::fabs(std::acos(v*surfacenor 393 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
388 if ( theta < 0.001 ) 394 if ( theta < 0.001 )
389 { 395 {
390 factor = 50 ; 396 factor = 50 ;
391 IsParallel = true ; 397 IsParallel = true ;
392 } 398 }
393 else 399 else
394 { 400 {
395 factor = 1 ; 401 factor = 1 ;
396 } 402 }
397 403
398 #ifdef G4TWISTDEBUG 404 #ifdef G4TWISTDEBUG
399 G4cout << "Step i = " << i << ", distanc 405 G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ;
400 G4cout << "X = " << tmpxx << G4endl ; 406 G4cout << "X = " << tmpxx << G4endl ;
401 #endif 407 #endif
402 408
403 GetPhiUAtX(tmpxx, phi, u) ; // the new p 409 GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced
404 410
405 #ifdef G4TWISTDEBUG 411 #ifdef G4TWISTDEBUG
406 G4cout << "approximated phi = " << phi < 412 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ;
407 #endif 413 #endif
408 414
409 if ( deltaX <= factor*ctol ) { IsConverg 415 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
410 416
411 } // end iterative loop (i) 417 } // end iterative loop (i)
412 418
413 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 419 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0. ;
414 420
415 #ifdef G4TWISTDEBUG 421 #ifdef G4TWISTDEBUG
416 G4cout << "refined solution " << phi << " 422 G4cout << "refined solution " << phi << " , " << u << G4endl ;
417 G4cout << "distance = " << tmpdist << G4en 423 G4cout << "distance = " << tmpdist << G4endl ;
418 G4cout << "local X = " << tmpxx << G4endl 424 G4cout << "local X = " << tmpxx << G4endl ;
419 #endif 425 #endif
420 426
421 tmpisvalid = false ; // init 427 tmpisvalid = false ; // init
422 428
423 if ( IsConverged ) 429 if ( IsConverged )
424 { 430 {
425 if (validate == kValidateWithTol) 431 if (validate == kValidateWithTol)
426 { 432 {
427 tmpareacode = GetAreaCode(tmpxx); 433 tmpareacode = GetAreaCode(tmpxx);
428 if (!IsOutside(tmpareacode)) 434 if (!IsOutside(tmpareacode))
429 { 435 {
430 if (tmpdist >= 0) tmpisvalid = true; 436 if (tmpdist >= 0) tmpisvalid = true;
431 } 437 }
432 } 438 }
433 else if (validate == kValidateWithoutTol 439 else if (validate == kValidateWithoutTol)
434 { 440 {
435 tmpareacode = GetAreaCode(tmpxx, false 441 tmpareacode = GetAreaCode(tmpxx, false);
436 if (IsInside(tmpareacode)) 442 if (IsInside(tmpareacode))
437 { 443 {
438 if (tmpdist >= 0) tmpisvalid = true; 444 if (tmpdist >= 0) tmpisvalid = true;
439 } 445 }
440 } 446 }
441 else // kDontValidate 447 else // kDontValidate
442 { 448 {
443 G4Exception("G4TwistBoxSide::DistanceT 449 G4Exception("G4TwistBoxSide::DistanceToSurface()",
444 "GeomSolids0001", FatalExc 450 "GeomSolids0001", FatalException,
445 "Feature NOT implemented ! 451 "Feature NOT implemented !");
446 } 452 }
447 } 453 }
448 else 454 else
449 { 455 {
450 tmpdist = kInfinity; // no convergen 456 tmpdist = kInfinity; // no convergence after 10 steps
451 tmpisvalid = false ; // solution is 457 tmpisvalid = false ; // solution is not vaild
452 } 458 }
453 459
454 // store the found values 460 // store the found values
455 k.xx = tmpxx ; << 461 xbuf[k].xx = tmpxx ;
456 k.distance = tmpdist ; << 462 xbuf[k].distance = tmpdist ;
457 k.areacode = tmpareacode ; << 463 xbuf[k].areacode = tmpareacode ;
458 k.isvalid = tmpisvalid ; << 464 xbuf[k].isvalid = tmpisvalid ;
459 465
460 } // end loop over physical solutions (vari 466 } // end loop over physical solutions (variable k)
461 467
462 std::sort(xbuf.begin() , xbuf.end(), Distanc 468 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting
463 469
464 #ifdef G4TWISTDEBUG 470 #ifdef G4TWISTDEBUG
465 G4cout << G4endl << "list xbuf after sorting 471 G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
466 G4cout << G4endl << G4endl ; 472 G4cout << G4endl << G4endl ;
467 #endif 473 #endif
468 474
469 // erase identical intersection (within kCar 475 // erase identical intersection (within kCarTolerance)
470 xbuf.erase(std::unique(xbuf.begin(),xbuf.end 476 xbuf.erase(std::unique(xbuf.begin(),xbuf.end(),EqualIntersection),xbuf.end());
471 477
472 // add guesses 478 // add guesses
473 479
474 auto nxxtmp = (G4int)xbuf.size() ; << 480 G4int nxxtmp = xbuf.size() ;
475 481
476 if ( nxxtmp<2 || IsParallel ) 482 if ( nxxtmp<2 || IsParallel )
477 { 483 {
478 // positive end 484 // positive end
479 #ifdef G4TWISTDEBUG 485 #ifdef G4TWISTDEBUG
480 G4cout << "add guess at +z/2 .. " << G4end 486 G4cout << "add guess at +z/2 .. " << G4endl ;
481 #endif 487 #endif
482 488
483 phi = fPhiTwist/2 ; 489 phi = fPhiTwist/2 ;
484 u = 0. ; 490 u = 0. ;
485 xbuftmp.phi = phi ; 491 xbuftmp.phi = phi ;
486 xbuftmp.u = u ; 492 xbuftmp.u = u ;
487 xbuftmp.areacode = sOutside ; 493 xbuftmp.areacode = sOutside ;
488 xbuftmp.distance = kInfinity ; 494 xbuftmp.distance = kInfinity ;
489 xbuftmp.isvalid = false ; 495 xbuftmp.isvalid = false ;
490 496
491 xbuf.push_back(xbuftmp) ; // store it to 497 xbuf.push_back(xbuftmp) ; // store it to xbuf
492 498
493 #ifdef G4TWISTDEBUG 499 #ifdef G4TWISTDEBUG
494 G4cout << "add guess at -z/2 .. " << G4end 500 G4cout << "add guess at -z/2 .. " << G4endl ;
495 #endif 501 #endif
496 502
497 phi = -fPhiTwist/2 ; 503 phi = -fPhiTwist/2 ;
498 u = 0. ; 504 u = 0. ;
499 505
500 xbuftmp.phi = phi ; 506 xbuftmp.phi = phi ;
501 xbuftmp.u = u ; 507 xbuftmp.u = u ;
502 xbuftmp.areacode = sOutside ; 508 xbuftmp.areacode = sOutside ;
503 xbuftmp.distance = kInfinity ; 509 xbuftmp.distance = kInfinity ;
504 xbuftmp.isvalid = false ; 510 xbuftmp.isvalid = false ;
505 511
506 xbuf.push_back(xbuftmp) ; // store it to 512 xbuf.push_back(xbuftmp) ; // store it to xbuf
507 513
508 for ( std::size_t k = nxxtmp ; k<xbuf.size << 514 for ( size_t k = nxxtmp ; k<xbuf.size() ; ++k )
509 { 515 {
510 #ifdef G4TWISTDEBUG 516 #ifdef G4TWISTDEBUG
511 G4cout << "Solution " << k << " : " 517 G4cout << "Solution " << k << " : "
512 << "reconstructed phiR = " << xbu 518 << "reconstructed phiR = " << xbuf[k].phi
513 << ", uR = " << xbuf[k].u << G4en 519 << ", uR = " << xbuf[k].u << G4endl ;
514 #endif 520 #endif
515 521
516 phi = xbuf[k].phi ; // get the stored v 522 phi = xbuf[k].phi ; // get the stored values for phi and u
517 u = xbuf[k].u ; 523 u = xbuf[k].u ;
518 524
519 IsConverged = false ; // no convergenc 525 IsConverged = false ; // no convergence at the beginning
520 526
521 for ( G4int i = 1 ; i<maxint ; ++i ) 527 for ( G4int i = 1 ; i<maxint ; ++i )
522 { 528 {
523 xxonsurface = SurfacePoint(phi,u) ; 529 xxonsurface = SurfacePoint(phi,u) ;
524 surfacenormal = NormAng(phi,u) ; 530 surfacenormal = NormAng(phi,u) ;
525 tmpdist = DistanceToPlaneWithV(p, v, x 531 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx);
526 deltaX = ( tmpxx - xxonsurface ).mag() 532 deltaX = ( tmpxx - xxonsurface ).mag() ;
527 theta = std::fabs(std::acos(v*surfacen 533 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
528 if ( theta < 0.001 ) 534 if ( theta < 0.001 )
529 { 535 {
530 factor = 50 ; 536 factor = 50 ;
531 } 537 }
532 else 538 else
533 { 539 {
534 factor = 1 ; 540 factor = 1 ;
535 } 541 }
536 542
537 #ifdef G4TWISTDEBUG 543 #ifdef G4TWISTDEBUG
538 G4cout << "Step i = " << i << ", dista 544 G4cout << "Step i = " << i << ", distance = "
539 << tmpdist << ", " << deltaX << 545 << tmpdist << ", " << deltaX << G4endl ;
540 G4cout << "X = " << tmpxx << G4endl ; 546 G4cout << "X = " << tmpxx << G4endl ;
541 #endif 547 #endif
542 548
543 GetPhiUAtX(tmpxx, phi, u) ; 549 GetPhiUAtX(tmpxx, phi, u) ;
544 // the new point xx is accepted and 550 // the new point xx is accepted and phi/u replaced
545 551
546 #ifdef G4TWISTDEBUG 552 #ifdef G4TWISTDEBUG
547 G4cout << "approximated phi = " << phi 553 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ;
548 #endif 554 #endif
549 555
550 if ( deltaX <= factor*ctol ) { IsConve 556 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
551 557
552 } // end iterative loop (i) 558 } // end iterative loop (i)
553 559
554 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 560 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0. ;
555 561
556 #ifdef G4TWISTDEBUG 562 #ifdef G4TWISTDEBUG
557 G4cout << "refined solution " << phi << 563 G4cout << "refined solution " << phi << " , " << u << G4endl ;
558 G4cout << "distance = " << tmpdist << G4 564 G4cout << "distance = " << tmpdist << G4endl ;
559 G4cout << "local X = " << tmpxx << G4end 565 G4cout << "local X = " << tmpxx << G4endl ;
560 #endif 566 #endif
561 567
562 tmpisvalid = false ; // init 568 tmpisvalid = false ; // init
563 569
564 if ( IsConverged ) 570 if ( IsConverged )
565 { 571 {
566 if (validate == kValidateWithTol) 572 if (validate == kValidateWithTol)
567 { 573 {
568 tmpareacode = GetAreaCode(tmpxx); 574 tmpareacode = GetAreaCode(tmpxx);
569 if (!IsOutside(tmpareacode)) 575 if (!IsOutside(tmpareacode))
570 { 576 {
571 if (tmpdist >= 0) tmpisvalid = tru 577 if (tmpdist >= 0) tmpisvalid = true;
572 } 578 }
573 } 579 }
574 else if (validate == kValidateWithoutT 580 else if (validate == kValidateWithoutTol)
575 { 581 {
576 tmpareacode = GetAreaCode(tmpxx, fal 582 tmpareacode = GetAreaCode(tmpxx, false);
577 if (IsInside(tmpareacode)) 583 if (IsInside(tmpareacode))
578 { 584 {
579 if (tmpdist >= 0) tmpisvalid = tru 585 if (tmpdist >= 0) tmpisvalid = true;
580 } 586 }
581 } 587 }
582 else // kDontValidate 588 else // kDontValidate
583 { 589 {
584 G4Exception("G4TwistedBoxSide::Dista 590 G4Exception("G4TwistedBoxSide::DistanceToSurface()",
585 "GeomSolids0001", FatalE 591 "GeomSolids0001", FatalException,
586 "Feature NOT implemented 592 "Feature NOT implemented !");
587 } 593 }
588 } 594 }
589 else 595 else
590 { 596 {
591 tmpdist = kInfinity; // no converg 597 tmpdist = kInfinity; // no convergence after 10 steps
592 tmpisvalid = false ; // solution i 598 tmpisvalid = false ; // solution is not vaild
593 } 599 }
594 600
595 // store the found values 601 // store the found values
596 xbuf[k].xx = tmpxx ; 602 xbuf[k].xx = tmpxx ;
597 xbuf[k].distance = tmpdist ; 603 xbuf[k].distance = tmpdist ;
598 xbuf[k].areacode = tmpareacode ; 604 xbuf[k].areacode = tmpareacode ;
599 xbuf[k].isvalid = tmpisvalid ; 605 xbuf[k].isvalid = tmpisvalid ;
600 } // end loop over physical solutions 606 } // end loop over physical solutions
601 } // end less than 2 solutions 607 } // end less than 2 solutions
602 608
603 // sort again 609 // sort again
604 std::sort(xbuf.begin() , xbuf.end(), Distanc 610 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting
605 611
606 // erase identical intersection (within kCar 612 // erase identical intersection (within kCarTolerance)
607 xbuf.erase(std::unique(xbuf.begin(),xbuf.end 613 xbuf.erase(std::unique(xbuf.begin(),xbuf.end(),EqualIntersection),xbuf.end());
608 614
609 #ifdef G4TWISTDEBUG 615 #ifdef G4TWISTDEBUG
610 G4cout << G4endl << "list xbuf after sorting 616 G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
611 G4cout << G4endl << G4endl ; 617 G4cout << G4endl << G4endl ;
612 #endif 618 #endif
613 619
614 nxx = (G4int)xbuf.size() ; // determine nu << 620 nxx = xbuf.size() ; // determine number of solutions again.
615 621
616 for ( G4int i = 0 ; i<(G4int)xbuf.size() ; + << 622 for ( size_t i = 0 ; i<xbuf.size() ; ++i )
617 { 623 {
618 distance[i] = xbuf[i].distance; 624 distance[i] = xbuf[i].distance;
619 gxx[i] = ComputeGlobalPoint(xbuf[i].x 625 gxx[i] = ComputeGlobalPoint(xbuf[i].xx);
620 areacode[i] = xbuf[i].areacode ; 626 areacode[i] = xbuf[i].areacode ;
621 isvalid[i] = xbuf[i].isvalid ; 627 isvalid[i] = xbuf[i].isvalid ;
622 628
623 fCurStatWithV.SetCurrentStatus(i, gxx[i], 629 fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i],
624 isvalid[i 630 isvalid[i], nxx, validate, &gp, &gv);
625 631
626 #ifdef G4TWISTDEBUG 632 #ifdef G4TWISTDEBUG
627 G4cout << "element Nr. " << i 633 G4cout << "element Nr. " << i
628 << ", local Intersection = " << xbu 634 << ", local Intersection = " << xbuf[i].xx
629 << ", distance = " << xbuf[i].dista 635 << ", distance = " << xbuf[i].distance
630 << ", u = " << xbuf[i].u 636 << ", u = " << xbuf[i].u
631 << ", phi = " << xbuf[i].phi 637 << ", phi = " << xbuf[i].phi
632 << ", isvalid = " << xbuf[i].isvali 638 << ", isvalid = " << xbuf[i].isvalid
633 << G4endl ; 639 << G4endl ;
634 #endif 640 #endif
635 641
636 } // end for( i ) loop 642 } // end for( i ) loop
637 643
638 #ifdef G4TWISTDEBUG 644 #ifdef G4TWISTDEBUG
639 G4cout << "G4TwistBoxSide finished " << G4en 645 G4cout << "G4TwistBoxSide finished " << G4endl ;
640 G4cout << nxx << " possible physical solutio 646 G4cout << nxx << " possible physical solutions found" << G4endl ;
641 for ( G4int k= 0 ; k< nxx ; ++k ) 647 for ( G4int k= 0 ; k< nxx ; ++k )
642 { 648 {
643 G4cout << "global intersection Point found 649 G4cout << "global intersection Point found: " << gxx[k] << G4endl ;
644 G4cout << "distance = " << distance[k] << 650 G4cout << "distance = " << distance[k] << G4endl ;
645 G4cout << "isvalid = " << isvalid[k] << G4 651 G4cout << "isvalid = " << isvalid[k] << G4endl ;
646 } 652 }
647 #endif 653 #endif
648 654
649 return nxx ; 655 return nxx ;
650 } 656 }
651 657
652 //============================================ 658 //=====================================================================
653 //* DistanceToSurface ------------------------ 659 //* DistanceToSurface -------------------------------------------------
654 660
655 G4int G4TwistBoxSide::DistanceToSurface(const 661 G4int G4TwistBoxSide::DistanceToSurface(const G4ThreeVector& gp,
656 662 G4ThreeVector gxx[],
657 663 G4double distance[],
658 664 G4int areacode[])
659 { 665 {
660 const G4double ctol = 0.5 * kCarTolerance; 666 const G4double ctol = 0.5 * kCarTolerance;
661 667
662 fCurStat.ResetfDone(kDontValidate, &gp); 668 fCurStat.ResetfDone(kDontValidate, &gp);
663 669
664 if (fCurStat.IsDone()) 670 if (fCurStat.IsDone())
665 { 671 {
666 for (G4int i=0; i<fCurStat.GetNXX(); ++i 672 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
667 { 673 {
668 gxx[i] = fCurStat.GetXX(i); 674 gxx[i] = fCurStat.GetXX(i);
669 distance[i] = fCurStat.GetDistance(i) 675 distance[i] = fCurStat.GetDistance(i);
670 areacode[i] = fCurStat.GetAreacode(i) 676 areacode[i] = fCurStat.GetAreacode(i);
671 } 677 }
672 return fCurStat.GetNXX(); 678 return fCurStat.GetNXX();
673 } 679 }
674 else // initialize 680 else // initialize
675 { 681 {
676 for (G4int i=0; i<G4VSURFACENXX; ++i) 682 for (G4int i=0; i<G4VSURFACENXX; ++i)
677 { 683 {
678 distance[i] = kInfinity; 684 distance[i] = kInfinity;
679 areacode[i] = sOutside; 685 areacode[i] = sOutside;
680 gxx[i].set(kInfinity, kInfinity, kInf 686 gxx[i].set(kInfinity, kInfinity, kInfinity);
681 } 687 }
682 } 688 }
683 689
684 G4ThreeVector p = ComputeLocalPoint(gp); 690 G4ThreeVector p = ComputeLocalPoint(gp);
685 G4ThreeVector xx; // intersection point 691 G4ThreeVector xx; // intersection point
686 G4ThreeVector xxonsurface ; // interpolated 692 G4ThreeVector xxonsurface ; // interpolated intersection point
687 693
688 // the surfacenormal at that surface point 694 // the surfacenormal at that surface point
689 G4double phiR = 0. ; 695 G4double phiR = 0. ;
690 G4double uR = 0. ; 696 G4double uR = 0. ;
691 697
692 G4ThreeVector surfacenormal ; 698 G4ThreeVector surfacenormal ;
693 G4double deltaX ; 699 G4double deltaX ;
694 700
695 G4int maxint = 20 ; 701 G4int maxint = 20 ;
696 702
697 for ( G4int i = 1 ; i<maxint ; ++i ) 703 for ( G4int i = 1 ; i<maxint ; ++i )
698 { 704 {
699 xxonsurface = SurfacePoint(phiR,uR) ; 705 xxonsurface = SurfacePoint(phiR,uR) ;
700 surfacenormal = NormAng(phiR,uR) ; 706 surfacenormal = NormAng(phiR,uR) ;
701 distance[0] = DistanceToPlane(p, xxonsurf 707 distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX
702 deltaX = ( xx - xxonsurface ).mag() ; 708 deltaX = ( xx - xxonsurface ).mag() ;
703 709
704 #ifdef G4TWISTDEBUG 710 #ifdef G4TWISTDEBUG
705 G4cout << "i = " << i << ", distance = " 711 G4cout << "i = " << i << ", distance = " << distance[0]
706 << ", " << deltaX << G4endl ; 712 << ", " << deltaX << G4endl ;
707 G4cout << "X = " << xx << G4endl ; 713 G4cout << "X = " << xx << G4endl ;
708 #endif 714 #endif
709 715
710 // the new point xx is accepted and phi/p 716 // the new point xx is accepted and phi/psi replaced
711 GetPhiUAtX(xx, phiR, uR) ; 717 GetPhiUAtX(xx, phiR, uR) ;
712 718
713 if ( deltaX <= ctol ) { break ; } 719 if ( deltaX <= ctol ) { break ; }
714 } 720 }
715 721
716 // check validity of solution ( valid phi,p 722 // check validity of solution ( valid phi,psi )
717 723
718 G4double halfphi = 0.5*fPhiTwist ; 724 G4double halfphi = 0.5*fPhiTwist ;
719 G4double uMax = GetBoundaryMax(phiR) ; 725 G4double uMax = GetBoundaryMax(phiR) ;
720 726
721 if ( phiR > halfphi ) phiR = halfphi ; 727 if ( phiR > halfphi ) phiR = halfphi ;
722 if ( phiR < -halfphi ) phiR = -halfphi ; 728 if ( phiR < -halfphi ) phiR = -halfphi ;
723 if ( uR > uMax ) uR = uMax ; 729 if ( uR > uMax ) uR = uMax ;
724 if ( uR < -uMax ) uR = -uMax ; 730 if ( uR < -uMax ) uR = -uMax ;
725 731
726 xxonsurface = SurfacePoint(phiR,uR) ; 732 xxonsurface = SurfacePoint(phiR,uR) ;
727 distance[0] = ( p - xx ).mag() ; 733 distance[0] = ( p - xx ).mag() ;
728 if ( distance[0] <= ctol ) { distance[0] = 734 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
729 735
730 // end of validity 736 // end of validity
731 737
732 #ifdef G4TWISTDEBUG 738 #ifdef G4TWISTDEBUG
733 G4cout << "refined solution " << phiR << " 739 G4cout << "refined solution " << phiR << " , " << uR << " , " << G4endl ;
734 G4cout << "distance = " << distance[0] << G 740 G4cout << "distance = " << distance[0] << G4endl ;
735 G4cout << "X = " << xx << G4endl ; 741 G4cout << "X = " << xx << G4endl ;
736 #endif 742 #endif
737 743
738 G4bool isvalid = true; 744 G4bool isvalid = true;
739 gxx[0] = ComputeGlobalPoint(xx); 745 gxx[0] = ComputeGlobalPoint(xx);
740 746
741 #ifdef G4TWISTDEBUG 747 #ifdef G4TWISTDEBUG
742 G4cout << "intersection Point found: " << g 748 G4cout << "intersection Point found: " << gxx[0] << G4endl ;
743 G4cout << "distance = " << distance[0] << G 749 G4cout << "distance = " << distance[0] << G4endl ;
744 #endif 750 #endif
745 751
746 fCurStat.SetCurrentStatus(0, gxx[0], distan 752 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
747 isvalid, 1, kDontV 753 isvalid, 1, kDontValidate, &gp);
748 return 1; 754 return 1;
749 } 755 }
750 756
751 //============================================ 757 //=====================================================================
752 //* GetAreaCode ------------------------------ 758 //* GetAreaCode -------------------------------------------------------
753 759
754 G4int G4TwistBoxSide::GetAreaCode(const G4Thre 760 G4int G4TwistBoxSide::GetAreaCode(const G4ThreeVector& xx,
755 G4bool 761 G4bool withTol)
756 { 762 {
757 // We must use the function in local coordi 763 // We must use the function in local coordinate system.
758 // See the description of DistanceToSurface 764 // See the description of DistanceToSurface(p,v).
759 765
760 const G4double ctol = 0.5 * kCarTolerance; 766 const G4double ctol = 0.5 * kCarTolerance;
761 767
762 G4double phi ; 768 G4double phi ;
763 G4double yprime ; 769 G4double yprime ;
764 GetPhiUAtX(xx, phi,yprime ) ; 770 GetPhiUAtX(xx, phi,yprime ) ;
765 771
766 G4double fYAxisMax = GetBoundaryMax(phi) ; 772 G4double fYAxisMax = GetBoundaryMax(phi) ; // Boundaries are symmetric
767 G4double fYAxisMin = - fYAxisMax ; 773 G4double fYAxisMin = - fYAxisMax ;
768 774
769 #ifdef G4TWISTDEBUG 775 #ifdef G4TWISTDEBUG
770 G4cout << "GetAreaCode: phi = " << phi << G 776 G4cout << "GetAreaCode: phi = " << phi << G4endl ;
771 G4cout << "GetAreaCode: yprime = " << yprim 777 G4cout << "GetAreaCode: yprime = " << yprime << G4endl ;
772 G4cout << "Intervall is " << fYAxisMin << " 778 G4cout << "Intervall is " << fYAxisMin << " to " << fYAxisMax << G4endl ;
773 #endif 779 #endif
774 780
775 G4int areacode = sInside; 781 G4int areacode = sInside;
776 782
777 if (fAxis[0] == kYAxis && fAxis[1] == kZAxi 783 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
778 { 784 {
779 G4int zaxis = 1; 785 G4int zaxis = 1;
780 786
781 if (withTol) 787 if (withTol)
782 { 788 {
783 G4bool isoutside = false; 789 G4bool isoutside = false;
784 790
785 // test boundary of yaxis 791 // test boundary of yaxis
786 792
787 if (yprime < fYAxisMin + ctol) 793 if (yprime < fYAxisMin + ctol)
788 { 794 {
789 areacode |= (sAxis0 & (sAxisY | sA 795 areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
790 if (yprime <= fYAxisMin - ctol) is 796 if (yprime <= fYAxisMin - ctol) isoutside = true;
791 797
792 } 798 }
793 else if (yprime > fYAxisMax - ctol) 799 else if (yprime > fYAxisMax - ctol)
794 { 800 {
795 areacode |= (sAxis0 & (sAxisY | sA 801 areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
796 if (yprime >= fYAxisMax + ctol) i 802 if (yprime >= fYAxisMax + ctol) isoutside = true;
797 } 803 }
798 804
799 // test boundary of z-axis 805 // test boundary of z-axis
800 806
801 if (xx.z() < fAxisMin[zaxis] + ctol) 807 if (xx.z() < fAxisMin[zaxis] + ctol)
802 { 808 {
803 areacode |= (sAxis1 & (sAxisZ | sA 809 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
804 810
805 if ((areacode & sBoundary) != 0) << 811 if (areacode & sBoundary) areacode |= sCorner; // xx on corner.
806 else areaco 812 else areacode |= sBoundary;
807 if (xx.z() <= fAxisMin[zaxis] - ct 813 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
808 814
809 } 815 }
810 else if (xx.z() > fAxisMax[zaxis] - c 816 else if (xx.z() > fAxisMax[zaxis] - ctol)
811 { 817 {
812 areacode |= (sAxis1 & (sAxisZ | sA 818 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
813 819
814 if ((areacode & sBoundary) != 0) << 820 if (areacode & sBoundary) areacode |= sCorner; // xx on corner.
815 else areaco 821 else areacode |= sBoundary;
816 if (xx.z() >= fAxisMax[zaxis] + ct 822 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
817 } 823 }
818 824
819 // if isoutside = true, clear inside 825 // if isoutside = true, clear inside bit.
820 // if not on boundary, add axis infor 826 // if not on boundary, add axis information.
821 827
822 if (isoutside) 828 if (isoutside)
823 { 829 {
824 G4int tmpareacode = areacode & (~s 830 G4int tmpareacode = areacode & (~sInside);
825 areacode = tmpareacode; 831 areacode = tmpareacode;
826 } 832 }
827 else if ((areacode & sBoundary) != sB 833 else if ((areacode & sBoundary) != sBoundary)
828 { 834 {
829 areacode |= (sAxis0 & sAxisY) | (s 835 areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
830 } 836 }
831 837
832 } 838 }
833 else 839 else
834 { 840 {
835 // boundary of y-axis 841 // boundary of y-axis
836 842
837 if (yprime < fYAxisMin ) 843 if (yprime < fYAxisMin )
838 { 844 {
839 areacode |= (sAxis0 & (sAxisY | sA 845 areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
840 } 846 }
841 else if (yprime > fYAxisMax) 847 else if (yprime > fYAxisMax)
842 { 848 {
843 areacode |= (sAxis0 & (sAxisY | sA 849 areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
844 } 850 }
845 851
846 // boundary of z-axis 852 // boundary of z-axis
847 853
848 if (xx.z() < fAxisMin[zaxis]) 854 if (xx.z() < fAxisMin[zaxis])
849 { 855 {
850 areacode |= (sAxis1 & (sAxisZ | sA 856 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
851 if ((areacode & sBoundary) != 0) << 857 if (areacode & sBoundary) areacode |= sCorner; // xx on corner.
852 else areaco 858 else areacode |= sBoundary;
853 859
854 } 860 }
855 else if (xx.z() > fAxisMax[zaxis]) 861 else if (xx.z() > fAxisMax[zaxis])
856 { 862 {
857 areacode |= (sAxis1 & (sAxisZ | sA 863 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
858 if ((areacode & sBoundary) != 0) << 864 if (areacode & sBoundary) areacode |= sCorner; // xx on corner.
859 else areaco 865 else areacode |= sBoundary;
860 } 866 }
861 867
862 if ((areacode & sBoundary) != sBounda 868 if ((areacode & sBoundary) != sBoundary)
863 { 869 {
864 areacode |= (sAxis0 & sAxisY) | (s 870 areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
865 } 871 }
866 } 872 }
867 return areacode; 873 return areacode;
868 } 874 }
869 else 875 else
870 { 876 {
871 G4Exception("G4TwistBoxSide::GetAreaCode 877 G4Exception("G4TwistBoxSide::GetAreaCode()",
872 "GeomSolids0001", FatalExcep 878 "GeomSolids0001", FatalException,
873 "Feature NOT implemented !") 879 "Feature NOT implemented !");
874 } 880 }
875 return areacode; 881 return areacode;
876 } 882 }
877 883
878 //============================================ 884 //=====================================================================
879 //* SetCorners() ----------------------------- 885 //* SetCorners() ------------------------------------------------------
880 886
881 void G4TwistBoxSide::SetCorners() 887 void G4TwistBoxSide::SetCorners()
882 { 888 {
883 889
884 // Set Corner points in local coodinate. 890 // Set Corner points in local coodinate.
885 891
886 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis 892 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
887 { 893 {
888 G4double x, y, z; 894 G4double x, y, z;
889 895
890 // corner of Axis0min and Axis1min 896 // corner of Axis0min and Axis1min
891 897
892 x = -fdeltaX/2. + (fDx2 - fDy1*fTAlph)*std 898 x = -fdeltaX/2. + (fDx2 - fDy1*fTAlph)*std::cos(fPhiTwist/2.) - fDy1*std::sin(fPhiTwist/2.) ;
893 y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/ 899 y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/2.) + (-fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
894 z = -fDz ; 900 z = -fDz ;
895 901
896 SetCorner(sC0Min1Min, x, y, z); 902 SetCorner(sC0Min1Min, x, y, z);
897 903
898 // corner of Axis0max and Axis1min 904 // corner of Axis0max and Axis1min
899 905
900 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std 906 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
901 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/ 907 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
902 z = -fDz ; 908 z = -fDz ;
903 909
904 SetCorner(sC0Max1Min, x, y, z); 910 SetCorner(sC0Max1Min, x, y, z);
905 911
906 // corner of Axis0max and Axis1max 912 // corner of Axis0max and Axis1max
907 913
908 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std: 914 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
909 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2 915 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
910 z = fDz ; 916 z = fDz ;
911 917
912 SetCorner(sC0Max1Max, x, y, z); 918 SetCorner(sC0Max1Max, x, y, z);
913 919
914 // corner of Axis0min and Axis1max 920 // corner of Axis0min and Axis1max
915 921
916 x = fdeltaX/2. + (fDx4 - fDy2*fTAlph)*std: 922 x = fdeltaX/2. + (fDx4 - fDy2*fTAlph)*std::cos(fPhiTwist/2.) + fDy2*std::sin(fPhiTwist/2.) ;
917 y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2 923 y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2.) + (fDx4 - fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
918 z = fDz ; 924 z = fDz ;
919 925
920 SetCorner(sC0Min1Max, x, y, z); 926 SetCorner(sC0Min1Max, x, y, z);
921 927
922 } 928 }
923 else 929 else
924 { 930 {
925 G4Exception("G4TwistBoxSide::SetCorners()" 931 G4Exception("G4TwistBoxSide::SetCorners()",
926 "GeomSolids0001", FatalExcepti 932 "GeomSolids0001", FatalException,
927 "Method NOT implemented !"); 933 "Method NOT implemented !");
928 } 934 }
929 } 935 }
930 936
931 //============================================ 937 //=====================================================================
932 //* SetBoundaries() -------------------------- 938 //* SetBoundaries() ---------------------------------------------------
933 939
934 void G4TwistBoxSide::SetBoundaries() 940 void G4TwistBoxSide::SetBoundaries()
935 { 941 {
936 // Set direction-unit vector of boundary-li 942 // Set direction-unit vector of boundary-lines in local coodinates.
937 943
938 G4ThreeVector direction; 944 G4ThreeVector direction;
939 945
940 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis 946 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
941 { 947 {
942 // sAxis0 & sAxisMin 948 // sAxis0 & sAxisMin
943 direction = GetCorner(sC0Min1Max) - GetCor 949 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
944 direction = direction.unit(); 950 direction = direction.unit();
945 SetBoundary(sAxis0 & (sAxisY | sAxisMin), 951 SetBoundary(sAxis0 & (sAxisY | sAxisMin), direction,
946 GetCorner(sC0Min1Min), sAxisZ) 952 GetCorner(sC0Min1Min), sAxisZ) ;
947 953
948 // sAxis0 & sAxisMax 954 // sAxis0 & sAxisMax
949 direction = GetCorner(sC0Max1Max) - GetCor 955 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
950 direction = direction.unit(); 956 direction = direction.unit();
951 SetBoundary(sAxis0 & (sAxisY | sAxisMax), 957 SetBoundary(sAxis0 & (sAxisY | sAxisMax), direction,
952 GetCorner(sC0Max1Min), sAxisZ) 958 GetCorner(sC0Max1Min), sAxisZ);
953 959
954 // sAxis1 & sAxisMin 960 // sAxis1 & sAxisMin
955 direction = GetCorner(sC0Max1Min) - GetCor 961 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
956 direction = direction.unit(); 962 direction = direction.unit();
957 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), 963 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
958 GetCorner(sC0Min1Min), sAxisY) 964 GetCorner(sC0Min1Min), sAxisY);
959 965
960 // sAxis1 & sAxisMax 966 // sAxis1 & sAxisMax
961 direction = GetCorner(sC0Max1Max) - GetCor 967 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
962 direction = direction.unit(); 968 direction = direction.unit();
963 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), 969 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
964 GetCorner(sC0Min1Max), sAxisY) 970 GetCorner(sC0Min1Max), sAxisY);
965 971
966 } 972 }
967 else 973 else
968 { 974 {
969 G4Exception("G4TwistBoxSide::SetCorners()" 975 G4Exception("G4TwistBoxSide::SetCorners()",
970 "GeomSolids0001", FatalExcepti 976 "GeomSolids0001", FatalException,
971 "Feature NOT implemented !"); 977 "Feature NOT implemented !");
972 } 978 }
973 } 979 }
974 980
975 981
976 void G4TwistBoxSide::GetPhiUAtX( const G4Three << 982 void G4TwistBoxSide::GetPhiUAtX( G4ThreeVector p, G4double& phi, G4double& u)
977 { 983 {
978 // find closest point XX on surface for a gi 984 // find closest point XX on surface for a given point p
979 // X0 is a point on the surface, d is the d 985 // X0 is a point on the surface, d is the direction
980 // ( both for a fixed z = pz) 986 // ( both for a fixed z = pz)
981 987
982 // phi is given by the z coordinate of p 988 // phi is given by the z coordinate of p
983 989
984 phi = p.z()/(2*fDz)*fPhiTwist ; 990 phi = p.z()/(2*fDz)*fPhiTwist ;
985 991
986 u = -(fTAlph*(fDx4plus2*fPhiTwist + 2*fDx4m 992 u = -(fTAlph*(fDx4plus2*fPhiTwist + 2*fDx4minus2*phi)
987 + 2*(fdeltaY*phi + fdeltaX*fTAlph*phi 993 + 2*(fdeltaY*phi + fdeltaX*fTAlph*phi
988 - fPhiTwist*(fTAlph*p.x() + p.y()))* 994 - fPhiTwist*(fTAlph*p.x() + p.y()))* std::cos(phi)
989 + 2*(-(fdeltaX*phi) + fdeltaY*fTAlph* 995 + 2*(-(fdeltaX*phi) + fdeltaY*fTAlph*phi + fPhiTwist*(p.x()
990 - fTAlph*p.y()))*std::sin(phi))/(2.*( 996 - fTAlph*p.y()))*std::sin(phi))/(2.*(fPhiTwist
991 + fPhiTwist*fTAlph*fTAlph)) ; 997 + fPhiTwist*fTAlph*fTAlph)) ;
992 } 998 }
993 999
994 1000
995 G4ThreeVector G4TwistBoxSide::ProjectPoint(con 1001 G4ThreeVector G4TwistBoxSide::ProjectPoint(const G4ThreeVector& p,
996 1002 G4bool isglobal)
997 { 1003 {
998 // Get Rho at p.z() on Hyperbolic Surface. 1004 // Get Rho at p.z() on Hyperbolic Surface.
999 G4ThreeVector tmpp; 1005 G4ThreeVector tmpp;
1000 if (isglobal) 1006 if (isglobal)
1001 { 1007 {
1002 tmpp = fRot.inverse()*p - fTrans; 1008 tmpp = fRot.inverse()*p - fTrans;
1003 } 1009 }
1004 else 1010 else
1005 { 1011 {
1006 tmpp = p; 1012 tmpp = p;
1007 } 1013 }
1008 1014
1009 G4double phi ; 1015 G4double phi ;
1010 G4double u ; 1016 G4double u ;
1011 1017
1012 GetPhiUAtX( tmpp, phi, u ) ; 1018 GetPhiUAtX( tmpp, phi, u ) ;
1013 // calculate (phi, u) for a point p close 1019 // calculate (phi, u) for a point p close the surface
1014 1020
1015 G4ThreeVector xx = SurfacePoint(phi,u) ; 1021 G4ThreeVector xx = SurfacePoint(phi,u) ;
1016 // transform back to cartesian coordinate 1022 // transform back to cartesian coordinates
1017 1023
1018 if (isglobal) 1024 if (isglobal)
1019 { 1025 {
1020 return (fRot * xx + fTrans); 1026 return (fRot * xx + fTrans);
1021 } 1027 }
1022 else 1028 else
1023 { 1029 {
1024 return xx; 1030 return xx;
1025 } 1031 }
1026 } 1032 }
1027 1033
1028 void G4TwistBoxSide::GetFacets( G4int k, G4in 1034 void G4TwistBoxSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
1029 G4int faces[] 1035 G4int faces[][4], G4int iside )
1030 { 1036 {
1031 G4double phi ; 1037 G4double phi ;
1032 G4double b ; 1038 G4double b ;
1033 1039
1034 G4double z, u ; // the two parameters f 1040 G4double z, u ; // the two parameters for the surface equation
1035 G4ThreeVector p ; // a point on the surfac 1041 G4ThreeVector p ; // a point on the surface, given by (z,u)
1036 1042
1037 G4int nnode ; 1043 G4int nnode ;
1038 G4int nface ; 1044 G4int nface ;
1039 1045
1040 // calculate the (n-1)*(k-1) vertices 1046 // calculate the (n-1)*(k-1) vertices
1041 1047
1042 G4int i,j ; 1048 G4int i,j ;
1043 1049
1044 for ( i = 0 ; i<n ; ++i ) 1050 for ( i = 0 ; i<n ; ++i )
1045 { 1051 {
1046 z = -fDz+i*(2.*fDz)/(n-1) ; 1052 z = -fDz+i*(2.*fDz)/(n-1) ;
1047 phi = z*fPhiTwist/(2*fDz) ; 1053 phi = z*fPhiTwist/(2*fDz) ;
1048 b = GetValueB(phi) ; 1054 b = GetValueB(phi) ;
1049 1055
1050 for ( j = 0 ; j<k ; ++j ) 1056 for ( j = 0 ; j<k ; ++j )
1051 { 1057 {
1052 nnode = GetNode(i,j,k,n,iside) ; 1058 nnode = GetNode(i,j,k,n,iside) ;
1053 u = -b/2 +j*b/(k-1) ; 1059 u = -b/2 +j*b/(k-1) ;
1054 p = SurfacePoint(phi,u,true) ; 1060 p = SurfacePoint(phi,u,true) ;
1055 // surface point in global coordinate 1061 // surface point in global coordinate system
1056 1062
1057 xyz[nnode][0] = p.x() ; 1063 xyz[nnode][0] = p.x() ;
1058 xyz[nnode][1] = p.y() ; 1064 xyz[nnode][1] = p.y() ;
1059 xyz[nnode][2] = p.z() ; 1065 xyz[nnode][2] = p.z() ;
1060 1066
1061 if ( i<n-1 && j<k-1 ) // contercloc 1067 if ( i<n-1 && j<k-1 ) // conterclock wise filling
1062 { 1068 {
1063 nface = GetFace(i,j,k,n,iside) ; 1069 nface = GetFace(i,j,k,n,iside) ;
1064 faces[nface][0] = GetEdgeVisibility(i 1070 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * (GetNode(i ,j ,k,n,iside)+1) ; // fortran numbering
1065 faces[nface][1] = GetEdgeVisibility(i 1071 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * (GetNode(i ,j+1,k,n,iside)+1) ;
1066 faces[nface][2] = GetEdgeVisibility(i 1072 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * (GetNode(i+1,j+1,k,n,iside)+1) ;
1067 faces[nface][3] = GetEdgeVisibility(i 1073 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * (GetNode(i+1,j ,k,n,iside)+1) ;
1068 } 1074 }
1069 } 1075 }
1070 } 1076 }
1071 } 1077 }
1072 1078