source: trunk/source/geometry/solids/specific/src/G4TwistTrapParallelSide.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

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