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