source: trunk/source/geometry/solids/specific/src/G4TwistTubsSide.cc @ 1058

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

file release beta
File size: 34.9 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: G4TwistTubsSide.cc,v 1.5 2006/06/29 18:49:18 gunter Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31// --------------------------------------------------------------------
32// GEANT 4 class source file
33//
34//
35// G4TwistTubsSide.cc
36//
37// Author:
38//   01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
39//
40// History:
41//   13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
42//                 from original version in Jupiter-2.5.02 application.
43//   29-Apr-2004 - O.Link. Bug fixed in GetAreaCode
44// --------------------------------------------------------------------
45
46#include "G4TwistTubsSide.hh"
47
48//=====================================================================
49//* constructors ------------------------------------------------------
50
51G4TwistTubsSide::G4TwistTubsSide(const G4String         &name,
52                                   const G4RotationMatrix &rot,
53                                   const G4ThreeVector    &tlate,
54                                         G4int             handedness,
55                                   const G4double          kappa,
56                                   const EAxis             axis0,
57                                   const EAxis             axis1,
58                                         G4double          axis0min,
59                                         G4double          axis1min,
60                                         G4double          axis0max,
61                                         G4double          axis1max)
62   : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
63                axis0min, axis1min, axis0max, axis1max),
64     fKappa(kappa)
65{
66   if (axis0 == kZAxis && axis1 == kXAxis) {
67      G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "InvalidSetup",
68                  FatalException, "Should swap axis0 and axis1!");
69   }
70   fIsValidNorm = false;
71   SetCorners();
72   SetBoundaries();
73}
74
75G4TwistTubsSide::G4TwistTubsSide(const G4String     &name,
76                                         G4double      EndInnerRadius[2],
77                                         G4double      EndOuterRadius[2],
78                                         G4double      DPhi,
79                                         G4double      EndPhi[2],
80                                         G4double      EndZ[2], 
81                                         G4double      InnerRadius,
82                                         G4double      OuterRadius,
83                                         G4double      Kappa,
84                                         G4int         handedness)
85  : G4VTwistSurface(name)
86{ 
87   fHandedness = handedness;   // +z = +ve, -z = -ve
88   fAxis[0]    = kXAxis; // in local coordinate system
89   fAxis[1]    = kZAxis;
90   fAxisMin[0] = InnerRadius;  // Inner-hype radius at z=0
91   fAxisMax[0] = OuterRadius;  // Outer-hype radius at z=0
92   fAxisMin[1] = EndZ[0];
93   fAxisMax[1] = EndZ[1];
94
95   fKappa = Kappa;
96   fRot.rotateZ( fHandedness > 0
97                 ? -0.5*DPhi
98                 :  0.5*DPhi );
99   fTrans.set(0, 0, 0);
100   fIsValidNorm = false;
101   
102   SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
103   SetBoundaries();
104}
105
106
107//=====================================================================
108//* Fake default constructor ------------------------------------------
109
110G4TwistTubsSide::G4TwistTubsSide( __void__& a )
111  : G4VTwistSurface(a)
112{
113}
114
115
116//=====================================================================
117//* destructor --------------------------------------------------------
118
119G4TwistTubsSide::~G4TwistTubsSide()
120{
121}
122
123//=====================================================================
124//* GetNormal ---------------------------------------------------------
125
126G4ThreeVector G4TwistTubsSide::GetNormal(const G4ThreeVector &tmpxx, 
127                                                G4bool isGlobal) 
128{
129   // GetNormal returns a normal vector at a surface (or very close
130   // to surface) point at tmpxx.
131   // If isGlobal=true, it returns the normal in global coordinate.
132   //
133   G4ThreeVector xx;
134   if (isGlobal) {
135      xx = ComputeLocalPoint(tmpxx);
136      if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
137         return ComputeGlobalDirection(fCurrentNormal.normal);
138      }
139   } else {
140      xx = tmpxx;
141      if (xx == fCurrentNormal.p) {
142         return fCurrentNormal.normal;
143      }
144   }
145   
146   G4ThreeVector er(1, fKappa * xx.z(), 0);
147   G4ThreeVector ez(0, fKappa * xx.x(), 1);
148   G4ThreeVector normal = fHandedness*(er.cross(ez));
149
150   if (isGlobal) {
151      fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
152   } else {
153      fCurrentNormal.normal = normal.unit();
154   }
155   return fCurrentNormal.normal;
156}
157
158//=====================================================================
159//* DistanceToSurface -------------------------------------------------
160
161G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp,
162                                          const G4ThreeVector &gv,
163                                                G4ThreeVector  gxx[],
164                                                G4double       distance[],
165                                                G4int          areacode[],
166                                                G4bool         isvalid[],
167                                                EValidate      validate)
168{
169   // Coordinate system:
170   //
171   //    The coordinate system is so chosen that the intersection of
172   //    the twisted surface with the z=0 plane coincides with the
173   //    x-axis.
174   //    Rotation matrix from this coordinate system (local system)
175   //    to global system is saved in fRot field.
176   //    So the (global) particle position and (global) velocity vectors,
177   //    p and v, should be rotated fRot.inverse() in order to convert
178   //    to local vectors.
179   //
180   // Equation of a twisted surface:
181   //
182   //    x(rho(z=0), z) = rho(z=0)
183   //    y(rho(z=0), z) = rho(z=0)*K*z
184   //    z(rho(z=0), z) = z
185   //    with
186   //       K = std::tan(fPhiTwist/2)/fZHalfLen
187   //
188   // Equation of a line:
189   //
190   //    gxx = p + t*v
191   //    with
192   //       p = fRot.inverse()*gp 
193   //       v = fRot.inverse()*gv
194   //
195   // Solution for intersection:
196   //
197   //    Required time for crossing is given by solving the
198   //    following quadratic equation:
199   //
200   //       a*t^2 + b*t + c = 0
201   //
202   //    where
203   //
204   //       a = K*v_x*v_z
205   //       b = K*(v_x*p_z + v_z*p_x) - v_y
206   //       c = K*p_x*p_z - p_y
207   //
208   //    Out of the possible two solutions you must choose
209   //    the one that gives a positive rho(z=0).
210   //
211   //
212     
213   fCurStatWithV.ResetfDone(validate, &gp, &gv);
214
215   if (fCurStatWithV.IsDone()) {
216      G4int i;
217      for (i=0; i<fCurStatWithV.GetNXX(); i++) {
218         gxx[i] = fCurStatWithV.GetXX(i);
219         distance[i] = fCurStatWithV.GetDistance(i);
220         areacode[i] = fCurStatWithV.GetAreacode(i);
221         isvalid[i]  = fCurStatWithV.IsValid(i);
222      }
223      return fCurStatWithV.GetNXX();
224   } else {
225      // initialize
226      G4int i;
227      for (i=0; i<2; 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   G4ThreeVector xx[2]; 
238
239
240   //
241   // special case!
242   // p is origin or
243   //
244
245   G4double absvz = std::fabs(v.z());
246
247   if ((absvz < DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x()) < DBL_MIN)) {
248      // no intersection
249
250      isvalid[0] = false;
251      fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
252                                isvalid[0], 0, validate, &gp, &gv);
253      return 0;
254   } 
255   
256   //
257   // special case end
258   //
259
260   
261   G4double a = fKappa * v.x() * v.z();
262   G4double b = fKappa * (v.x() * p.z() + v.z() * p.x()) - v.y();
263   G4double c = fKappa * p.x() * p.z() - p.y();
264   G4double D = b * b - 4 * a * c;             // discriminant
265
266   if (std::fabs(a) < DBL_MIN) {
267      if (std::fabs(b) > DBL_MIN) { 
268
269         // single solution
270
271         distance[0] = - c / b;
272         xx[0]       = p + distance[0]*v;
273         gxx[0]      = ComputeGlobalPoint(xx[0]);
274
275         if (validate == kValidateWithTol) {
276            areacode[0] = GetAreaCode(xx[0]);
277            if (!IsOutside(areacode[0])) {
278               if (distance[0] >= 0) isvalid[0] = true;
279            }
280         } else if (validate == kValidateWithoutTol) {
281            areacode[0] = GetAreaCode(xx[0], false);
282            if (IsInside(areacode[0])) {
283               if (distance[0] >= 0) isvalid[0] = true;
284            }
285         } else { // kDontValidate                       
286            // we must omit x(rho,z) = rho(z=0) < 0
287            if (xx[0].x() > 0) {
288               areacode[0] = sInside;
289               if (distance[0] >= 0) isvalid[0] = true;
290            } else {
291               distance[0] = kInfinity;
292               fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
293                                              areacode[0], isvalid[0],
294                                              0, validate, &gp, &gv);
295               return 0;
296            } 
297         }
298                 
299         fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
300                                        isvalid[0], 1, validate, &gp, &gv);
301         return 1;
302
303      } else {
304         // if a=b=0 , v.y=0 and (v.x=0 && p.x=0) or (v.z=0 && p.z=0) .
305         //    if v.x=0 && p.x=0, no intersection unless p is on z-axis
306         //    (in that case, v is paralell to surface).
307         //    if v.z=0 && p.z=0, no intersection unless p is on x-axis
308         //    (in that case, v is paralell to surface).
309         // return distance = infinity.
310
311         fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
312                                        isvalid[0], 0, validate, &gp, &gv);
313
314         return 0;
315      }
316     
317   } else if (D > DBL_MIN) {   
318
319      // double solutions
320
321      D = std::sqrt(D);
322      G4double      factor = 0.5/a;
323      G4double      tmpdist[2] = {kInfinity, kInfinity};
324      G4ThreeVector tmpxx[2];
325      G4int         tmpareacode[2] = {sOutside, sOutside};
326      G4bool        tmpisvalid[2]  = {false, false};
327      G4int i;
328
329      for (i=0; i<2; i++) {
330         G4double bminusD = - b - D;
331
332         // protection against round off error 
333         //G4double protection = 1.0e-6;
334         G4double protection = 0;
335         if ( b * D < 0 && std::fabs(bminusD / D) < protection ) {
336            G4double acovbb = (a*c)/(b*b);
337            tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
338         } else { 
339            tmpdist[i] = factor * bminusD;
340         }
341
342         D = -D;
343         tmpxx[i] = p + tmpdist[i]*v;
344         
345         if (validate == kValidateWithTol) {
346            tmpareacode[i] = GetAreaCode(tmpxx[i]);
347            if (!IsOutside(tmpareacode[i])) {
348               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
349               continue;
350            }
351         } else if (validate == kValidateWithoutTol) {
352            tmpareacode[i] = GetAreaCode(tmpxx[i], false);
353            if (IsInside(tmpareacode[i])) {
354               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
355               continue;
356            }
357         } else { // kDontValidate
358            // we must choose x(rho,z) = rho(z=0) > 0
359            if (tmpxx[i].x() > 0) {
360               tmpareacode[i] = sInside;
361               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
362               continue;
363            } else {
364               tmpdist[i] = kInfinity;
365               continue;
366            }                     
367         }
368      }
369     
370      if (tmpdist[0] <= tmpdist[1]) {
371         distance[0] = tmpdist[0];
372         distance[1] = tmpdist[1];
373         xx[0]       = tmpxx[0];
374         xx[1]       = tmpxx[1];
375         gxx[0]      = ComputeGlobalPoint(tmpxx[0]);
376         gxx[1]      = ComputeGlobalPoint(tmpxx[1]);
377         areacode[0] = tmpareacode[0];
378         areacode[1] = tmpareacode[1];
379         isvalid[0]  = tmpisvalid[0];
380         isvalid[1]  = tmpisvalid[1];
381      } else {
382         distance[0] = tmpdist[1];
383         distance[1] = tmpdist[0];
384         xx[0]       = tmpxx[1];
385         xx[1]       = tmpxx[0];
386         gxx[0]      = ComputeGlobalPoint(tmpxx[1]);
387         gxx[1]      = ComputeGlobalPoint(tmpxx[0]);
388         areacode[0] = tmpareacode[1];
389         areacode[1] = tmpareacode[0];
390         isvalid[0]  = tmpisvalid[1];
391         isvalid[1]  = tmpisvalid[0];
392      }
393         
394      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
395                                     isvalid[0], 2, validate, &gp, &gv);
396      fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
397                                     isvalid[1], 2, validate, &gp, &gv);
398
399      // protection against roundoff error
400
401      for (G4int k=0; k<2; k++) {
402         if (!isvalid[k]) continue;
403
404         G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
405                                              * xx[k].z() , xx[k].z());
406         G4double      deltaY  =  (xx[k] - xxonsurface).mag();
407
408         if ( deltaY > 0.5*kCarTolerance ) {
409
410           G4int maxcount = 10;
411           G4int l;
412           G4double      lastdeltaY = deltaY; 
413           G4ThreeVector last = deltaY; 
414           for (l=0; l<maxcount; l++) {
415             G4ThreeVector surfacenormal = GetNormal(xxonsurface); 
416             distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
417                                                surfacenormal, xx[k]);
418             deltaY      = (xx[k] - xxonsurface).mag();
419             if (deltaY > lastdeltaY) {
420               
421             }
422             gxx[k]      = ComputeGlobalPoint(xx[k]);
423
424               if (deltaY <= 0.5*kCarTolerance) {
425           
426                 break;
427               }
428               xxonsurface.set(xx[k].x(),
429                               fKappa * std::fabs(xx[k].x()) * xx[k].z(),
430                               xx[k].z());
431            }
432            if (l == maxcount) {
433               G4cerr << "ERROR - G4TwistTubsSide::DistanceToSurface(p,v)"
434                      << G4endl
435                      << "        maxloop count " << maxcount << G4endl;
436               G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
437                           "InvalidSetup",  FatalException,
438                           "Exceeded maxloop count!");
439            }
440         }
441
442      }
443      return 2;
444   } else {
445      // if D<0, no solution
446      // if D=0, just grazing the surfaces, return kInfinity
447
448      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
449                                     isvalid[0], 0, validate, &gp, &gv);
450
451      return 0;
452   }
453   G4Exception("G4TwistTubsSide::DistanceToSurface(p,v)",
454               "InvalidCondition", FatalException, "Illegal operation !");
455   return 1;
456}
457
458//=====================================================================
459//* DistanceToSurface -------------------------------------------------
460
461G4int G4TwistTubsSide::DistanceToSurface(const G4ThreeVector &gp,
462                                                G4ThreeVector  gxx[],
463                                                G4double       distance[],
464                                                G4int          areacode[])
465{ 
466   fCurStat.ResetfDone(kDontValidate, &gp);
467   G4int i = 0;
468   if (fCurStat.IsDone()) {
469      for (i=0; i<fCurStat.GetNXX(); i++) {
470         gxx[i] = fCurStat.GetXX(i);
471         distance[i] = fCurStat.GetDistance(i);
472         areacode[i] = fCurStat.GetAreacode(i);
473      }
474      return fCurStat.GetNXX();
475   } else {
476      // initialize
477      for (i=0; i<2; i++) {
478         distance[i] = kInfinity;
479         areacode[i] = sOutside;
480         gxx[i].set(kInfinity, kInfinity, kInfinity);
481      }
482   }
483   
484   static const G4double halftol = 0.5 * kCarTolerance; 
485
486   G4ThreeVector  p       = ComputeLocalPoint(gp);
487   G4ThreeVector  xx;
488   G4int          parity  = (fKappa >= 0 ? 1 : -1);
489 
490   //
491   // special case!
492   // If p is on surface, or
493   // p is on z-axis,
494   // return here immediatery.
495   //
496   
497   G4ThreeVector  lastgxx[2];
498   G4double       distfromlast[2];
499   for (i=0; i<2; i++) {
500      lastgxx[i] = fCurStatWithV.GetXX(i);
501      distfromlast[i] = (gp - lastgxx[i]).mag();
502   } 
503
504   if  ((gp - lastgxx[0]).mag() < halftol
505     || (gp - lastgxx[1]).mag() < halftol) { 
506      // last winner, or last poststep point is on the surface.
507      xx = p;
508      distance[0] = 0;
509      gxx[0] = gp;
510
511      G4bool isvalid = true;
512      fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
513                             isvalid, 1, kDontValidate, &gp);
514      return 1;
515   }
516         
517   if (p.getRho() == 0) { 
518      // p is on z-axis. Namely, p is on twisted surface (invalid area).
519      // We must return here, however, returning distance to x-minimum
520      // boundary is better than return 0-distance.
521      //
522      G4bool isvalid = true;
523      if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
524         distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
525         areacode[0] = sInside;
526      } else {
527         distance[0] = 0;
528         xx.set(0., 0., 0.);
529      }
530      gxx[0] = ComputeGlobalPoint(xx);
531      fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
532                                isvalid, 0, kDontValidate, &gp);
533      return 1;
534   } 
535
536   //
537   // special case end
538   //
539
540   // set corner points of quadrangle try area ...
541
542   G4ThreeVector A;  // foot of normal from p to boundary of sAxis0 & sAxisMin
543   G4ThreeVector C;  // foot of normal from p to boundary of sAxis0 & sAxisMax
544   G4ThreeVector B;       // point on boundary sAxis0 & sAxisMax at z = A.z()
545   G4ThreeVector D;       // point on boundary sAxis0 & sAxisMin at z = C.z()
546   G4double      distToA; // distance from p to A
547   G4double      distToC; // distance from p to C
548
549   distToA = DistanceToBoundary(sAxis0 & sAxisMin, A, p);
550   distToC = DistanceToBoundary(sAxis0 & sAxisMax, C, p);
551   
552   // is p.z between a.z and c.z?
553   // p.z must be bracketed a.z and c.z.
554   if (A.z() > C.z()) {
555      if (p.z() > A.z()) {
556         A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
557      } else if (p.z() < C.z()) {
558         C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
559      }
560   } else {
561      if (p.z() > C.z()) {
562         C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
563      } else if (p.z() < A.z()) {
564         A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
565      }
566   }
567
568   G4ThreeVector  d[2];     // direction vectors of boundary
569   G4ThreeVector  x0[2];    // foot of normal from line to p
570   G4int          btype[2]; // boundary type
571
572   for (i=0; i<2; i++) {
573      if (i == 0) {
574         GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
575         B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i]; 
576         // x0 + t*d , d is direction unit vector.
577      } else {
578         GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
579         D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i]; 
580      }
581   }
582
583   // In order to set correct diagonal, swap A and D, C and B if needed. 
584   G4ThreeVector pt(p.x(), p.y(), 0.);
585   G4double      rc = std::fabs(p.x());
586   G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.); 
587   G4int         pside = AmIOnLeftSide(pt, surfacevector); 
588   G4double      test  = (A.z() - C.z()) * parity * pside; 
589
590   if (test == 0) {
591      if (pside == 0) {
592         // p is on surface.
593         xx = p;
594         distance[0] = 0;
595         gxx[0] = gp;
596
597         G4bool isvalid = true;
598         fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
599                                isvalid, 1, kDontValidate, &gp);
600         return 1;
601      } else {
602         // A.z = C.z(). return distance to line.
603         d[0] = C - A;
604         distance[0] = DistanceToLine(p, A, d[0], xx);
605         areacode[0] = sInside;
606         gxx[0] = ComputeGlobalPoint(xx);
607         G4bool isvalid = true;
608         fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
609                                isvalid, 1, kDontValidate, &gp);
610         return 1;
611      } 
612
613   } else if (test < 0) {
614
615      // wrong diagonal. vector AC is crossing the surface! 
616      // swap A and D, C and B
617      G4ThreeVector tmp;
618      tmp = A;
619      A = D;
620      D = tmp;
621      tmp = C;
622      C = B;
623      B = tmp; 
624
625   } else {
626      // correct diagonal. nothing to do. 
627   }
628
629   // Now, we chose correct diaglnal.
630   // First try. divide quadrangle into double triangle by diagonal and
631   // calculate distance to both surfaces.
632
633   G4ThreeVector xxacb;   // foot of normal from plane ACB to p
634   G4ThreeVector nacb;    // normal of plane ACD
635   G4ThreeVector xxcad;   // foot of normal from plane CAD to p
636   G4ThreeVector ncad;    // normal of plane CAD
637   G4ThreeVector AB(A.x(), A.y(), 0);
638   G4ThreeVector DC(C.x(), C.y(), 0);
639
640   G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB, xxacb, nacb) * parity;
641   G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC, xxcad, ncad) * parity;
642
643   // if calculated distance = 0, return 
644
645   if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol) {
646      xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad); 
647      areacode[0] = sInside;
648      gxx[0] = ComputeGlobalPoint(xx);
649      distance[0] = 0;
650      G4bool isvalid = true;
651      fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
652                                isvalid, 1, kDontValidate, &gp);
653      return 1;
654   }
655   
656   if (distToACB * distToCAD > 0 && distToACB < 0) {
657      // both distToACB and distToCAD are negative.
658      // divide quadrangle into double triangle by diagonal
659      G4ThreeVector normal;
660      distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
661   } else {
662      if (distToACB * distToCAD > 0) {
663         // both distToACB and distToCAD are positive.
664         // Take smaller one.
665         if (distToACB <= distToCAD) {
666            distance[0] = distToACB;
667            xx   = xxacb;
668         } else {
669            distance[0] = distToCAD;
670            xx   = xxcad;
671         }
672      } else {
673         // distToACB * distToCAD is negative.
674         // take positive one
675         if (distToACB > 0) {
676            distance[0] = distToACB;
677            xx   = xxacb;
678         } else {
679            distance[0] = distToCAD;
680            xx   = xxcad;
681         }
682      }
683     
684   }
685   areacode[0] = sInside;
686   gxx[0]      = ComputeGlobalPoint(xx);
687   G4bool isvalid = true;
688   fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
689                             isvalid, 1, kDontValidate, &gp);
690   return 1;
691}
692
693//=====================================================================
694//* DistanceToPlane ---------------------------------------------------
695
696G4double G4TwistTubsSide::DistanceToPlane(const G4ThreeVector &p,
697                                           const G4ThreeVector &A,
698                                           const G4ThreeVector &B,
699                                           const G4ThreeVector &C,
700                                           const G4ThreeVector &D,
701                                           const G4int          parity,
702                                                 G4ThreeVector &xx,
703                                                 G4ThreeVector &n)
704{
705   static const G4double halftol = 0.5 * kCarTolerance;
706   
707   G4ThreeVector M = 0.5*(A + B);
708   G4ThreeVector N = 0.5*(C + D);
709   G4ThreeVector xxanm;  // foot of normal from p to plane ANM
710   G4ThreeVector nanm;   // normal of plane ANM
711   G4ThreeVector xxcmn;  // foot of normal from p to plane CMN
712   G4ThreeVector ncmn;   // normal of plane CMN
713
714   G4double distToanm = G4VTwistSurface::DistanceToPlane(p, A, (N - A), (M - A), xxanm, nanm) * parity;
715   G4double distTocmn = G4VTwistSurface::DistanceToPlane(p, C, (M - C), (N - C), xxcmn, ncmn) * parity;
716
717   // if p is behind of both surfaces, abort.
718   if (distToanm * distTocmn > 0 && distToanm < 0) {
719     G4Exception("G4TwistTubsSide::DistanceToPlane()",
720                 "InvalidCondition", FatalException,
721                 "Point p is behind the surfaces.");
722   }
723
724   // if p is on surface, return 0.
725   if (std::fabs(distToanm) <= halftol) {
726      xx = xxanm;
727      n  = nanm * parity;
728      return 0;
729   } else if (std::fabs(distTocmn) <= halftol) {
730      xx = xxcmn;
731      n  = ncmn * parity;
732      return 0;
733   }
734   
735   if (distToanm <= distTocmn) {
736      if (distToanm > 0) {
737         // both distanses are positive. take smaller one.
738         xx = xxanm;
739         n  = nanm * parity;
740         return distToanm;
741      } else {
742         // take -ve distance and call the function recursively.
743         return DistanceToPlane(p, A, M, N, D, parity, xx, n);
744      }
745   } else {
746      if (distTocmn > 0) {
747         // both distanses are positive. take smaller one.
748         xx = xxcmn;
749         n  = ncmn * parity;
750         return distTocmn;
751      } else {
752         // take -ve distance and call the function recursively.
753         return DistanceToPlane(p, C, N, M, B, parity, xx, n);
754      }
755   }
756}
757
758//=====================================================================
759//* GetAreaCode -------------------------------------------------------
760
761G4int G4TwistTubsSide::GetAreaCode(const G4ThreeVector &xx, 
762                                          G4bool withTol)
763{
764   // We must use the function in local coordinate system.
765   // See the description of DistanceToSurface(p,v).
766   
767   static const G4double ctol = 0.5 * kCarTolerance;
768   G4int areacode = sInside;
769   
770   if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
771      G4int xaxis = 0;
772      G4int zaxis = 1;
773     
774      if (withTol) {
775
776         G4bool isoutside   = false;
777
778         // test boundary of xaxis
779
780         if (xx.x() < fAxisMin[xaxis] + ctol) {
781            areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; 
782            if (xx.x() <= fAxisMin[xaxis] - ctol) isoutside = true;
783
784         } else if (xx.x() > fAxisMax[xaxis] - ctol) {
785            areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
786            if (xx.x() >= fAxisMax[xaxis] + ctol)  isoutside = true;
787         }
788
789         // test boundary of z-axis
790
791         if (xx.z() < fAxisMin[zaxis] + ctol) {
792            areacode |= (sAxis1 & (sAxisZ | sAxisMin)); 
793
794            if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
795            else                        areacode |= sBoundary;
796            if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
797
798         } else if (xx.z() > fAxisMax[zaxis] - ctol) {
799            areacode |= (sAxis1 & (sAxisZ | sAxisMax));
800
801            if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
802            else                        areacode |= sBoundary; 
803            if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
804         }
805
806         // if isoutside = true, clear inside bit.             
807         // if not on boundary, add axis information.             
808         
809         if (isoutside) {
810            G4int tmpareacode = areacode & (~sInside);
811            areacode = tmpareacode;
812         } else if ((areacode & sBoundary) != sBoundary) {
813            areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
814         }           
815         
816      } else {
817
818         // boundary of x-axis
819
820         if (xx.x() < fAxisMin[xaxis] ) {
821            areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
822         } else if (xx.x() > fAxisMax[xaxis]) {
823            areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
824         }
825         
826         // boundary of z-axis
827
828         if (xx.z() < fAxisMin[zaxis]) {
829            areacode |= (sAxis1 & (sAxisZ | sAxisMin));
830            if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
831            else                        areacode |= sBoundary; 
832           
833         } else if (xx.z() > fAxisMax[zaxis]) {
834            areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
835            if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
836            else                        areacode |= sBoundary; 
837         }
838
839         if ((areacode & sBoundary) != sBoundary) {
840            areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
841         }           
842      }
843      return areacode;
844   } else {
845      G4Exception("G4TwistTubsSide::GetAreaCode()",
846                  "NotImplemented", FatalException,
847                  "Feature NOT implemented !");
848   }
849   return areacode;
850}
851
852//=====================================================================
853//* SetCorners( arglist ) -------------------------------------------------
854
855void G4TwistTubsSide::SetCorners(
856                                  G4double      endInnerRad[2],
857                                  G4double      endOuterRad[2],
858                                  G4double      endPhi[2],
859                                  G4double      endZ[2])
860{
861   // Set Corner points in local coodinate.   
862
863   if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
864   
865      G4int zmin = 0 ;  // at -ve z
866      G4int zmax = 1 ;  // at +ve z
867
868      G4double x, y, z;
869     
870      // corner of Axis0min and Axis1min
871      x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
872      y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
873      z = endZ[zmin];
874      SetCorner(sC0Min1Min, x, y, z);
875     
876      // corner of Axis0max and Axis1min
877      x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
878      y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
879      z = endZ[zmin];
880      SetCorner(sC0Max1Min, x, y, z);
881     
882      // corner of Axis0max and Axis1max
883      x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
884      y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
885      z = endZ[zmax];
886      SetCorner(sC0Max1Max, x, y, z);
887     
888      // corner of Axis0min and Axis1max
889      x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
890      y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
891      z = endZ[zmax];
892      SetCorner(sC0Min1Max, x, y, z);
893
894   } else {
895      G4cerr << "ERROR - G4TwistTubsFlatSide::SetCorners()" << G4endl
896             << "        fAxis[0] = " << fAxis[0] << G4endl
897             << "        fAxis[1] = " << fAxis[1] << G4endl;
898      G4Exception("G4TwistTubsSide::SetCorners()",
899                  "NotImplemented", FatalException,
900                  "Feature NOT implemented !");
901   }
902}
903
904//=====================================================================
905//* SetCorners() ------------------------------------------------------
906
907void G4TwistTubsSide::SetCorners()
908{
909   G4Exception("G4TwistTubsSide::SetCorners()",
910               "NotImplemented", FatalException,
911               "Method NOT implemented !");
912}
913
914//=====================================================================
915//* SetBoundaries() ---------------------------------------------------
916
917void G4TwistTubsSide::SetBoundaries()
918{
919   // Set direction-unit vector of boundary-lines in local coodinate.
920   //   
921   G4ThreeVector direction;
922   
923   if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
924     
925      // sAxis0 & sAxisMin
926      direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
927      direction = direction.unit();
928      SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction, 
929                  GetCorner(sC0Min1Min), sAxisZ) ;
930     
931      // sAxis0 & sAxisMax
932      direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
933      direction = direction.unit();
934      SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction, 
935                  GetCorner(sC0Max1Min), sAxisZ);
936                 
937      // sAxis1 & sAxisMin
938      direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
939      direction = direction.unit();
940      SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction, 
941                  GetCorner(sC0Min1Min), sAxisX);
942                 
943      // sAxis1 & sAxisMax
944      direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
945      direction = direction.unit();
946      SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction, 
947                  GetCorner(sC0Min1Max), sAxisX);
948                 
949   } else {
950      G4cerr << "ERROR - G4TwistTubsFlatSide::SetBoundaries()" << G4endl
951             << "        fAxis[0] = " << fAxis[0] << G4endl
952             << "        fAxis[1] = " << fAxis[1] << G4endl;
953      G4Exception("G4TwistTubsSide::SetCorners()",
954                  "NotImplemented", FatalException,
955                  "Feature NOT implemented !");
956   }
957}
958
959//=====================================================================
960//* GetFacets() -------------------------------------------------------
961
962void G4TwistTubsSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
963                                 G4int faces[][4], G4int iside ) 
964{
965
966  G4double z ;     // the two parameters for the surface equation
967  G4double x,xmin,xmax ;
968
969  G4ThreeVector p ;  // a point on the surface, given by (z,u)
970
971  G4int nnode ;
972  G4int nface ;
973
974  // calculate the (n-1)*(m-1) vertices
975
976  G4int i,j ;
977
978  for ( i = 0 ; i<n ; i++ )
979  {
980
981    z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
982
983    for ( j = 0 ; j<m ; j++ ) {
984
985      nnode = GetNode(i,j,m,n,iside) ;
986
987      xmin = GetBoundaryMin(z) ; 
988      xmax = GetBoundaryMax(z) ;
989
990      if (fHandedness < 0) { 
991        x = xmin + j*(xmax-xmin)/(m-1) ;
992      } else {               
993        x = xmax - j*(xmax-xmin)/(m-1) ;
994      }
995
996      p = SurfacePoint(x,z,true) ;  // surface point in global coord.system
997
998      xyz[nnode][0] = p.x() ;
999      xyz[nnode][1] = p.y() ;
1000      xyz[nnode][2] = p.z() ;
1001
1002      if ( i<n-1 && j<m-1 ) {   // clock wise filling
1003       
1004        nface = GetFace(i,j,m,n,iside) ;
1005
1006        faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i  ,j  ,m,n,iside)+1) ; 
1007        faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,1) * ( GetNode(i+1,j  ,m,n,iside)+1) ;
1008        faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
1009        faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,1) * ( GetNode(i  ,j+1,m,n,iside)+1) ;
1010
1011      }
1012    }
1013  }
1014}
Note: See TracBrowser for help on using the repository browser.