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