source: trunk/source/geometry/solids/specific/src/G4TwistBoxSide.cc @ 1202

Last change on this file since 1202 was 1058, checked in by garnier, 15 years ago

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