2 // ********************************************************************
 
    3 // * License and Disclaimer                                           *
 
    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.                             *
 
   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.         *
 
   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 // ********************************************************************
 
   27 // $Id: G4VTwistSurface.icc 66356 2012-12-18 09:02:32Z gcosmo $
 
   30 // --------------------------------------------------------------------
 
   31 // G4VTwistSurface class inline methods
 
   34 //   01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
 
   37 //   13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
 
   38 //                 from original version in Jupiter-2.5.02 application.
 
   39 // --------------------------------------------------------------------
 
   41 //=====================================================================
 
   42 //* DistanceToPlaneWithV ----------------------------------------------
 
   45 G4double G4VTwistSurface::DistanceToPlaneWithV(const G4ThreeVector &p,
 
   46                                           const G4ThreeVector &v,
 
   47                                           const G4ThreeVector &x0,
 
   48                                           const G4ThreeVector &n0,
 
   52    G4double t = kInfinity;
 
   53    if (q) { t = (n0 * (x0 - p)) / q; }
 
   58 //=====================================================================
 
   59 //* DistanceToPlane ---------------------------------------------------
 
   62 G4double G4VTwistSurface::DistanceToPlane(const G4ThreeVector &p,
 
   63                                      const G4ThreeVector &x0,
 
   64                                      const G4ThreeVector &n0,
 
   68    // Calculate distance to plane in local coordinate,
 
   69    // then return distance and global intersection points.
 
   71    // p          - location of flying particle   
 
   72    // x0         - reference point of surface
 
   73    // xx         - a foot of perpendicular line from p to the plane 
 
   74    // t          - distance from xx to p
 
   75    // n          - a unit normal of this plane from plane to p. 
 
   84    //         t = n * (p - x0) / std::fabs(n)
 
   87    G4ThreeVector n = n0.unit();
 
   93 //=====================================================================
 
   94 //* DistanceToPlane ---------------------------------------------------
 
   97 G4double G4VTwistSurface::DistanceToPlane(const G4ThreeVector &p,
 
   98                                      const G4ThreeVector &x0,
 
   99                                      const G4ThreeVector &t1,
 
  100                                      const G4ThreeVector &t2,
 
  105    // Calculate distance to plane in local coordinate,
 
  106    // then return distance and global intersection points.
 
  107    // t1         - 1st. vector lying on the plane 
 
  108    // t2         - 2nd. vector lying on the plane
 
  110    n = (t1.cross(t2)).unit();
 
  111    return DistanceToPlane(p, x0, n, xx); 
 
  114 //=====================================================================
 
  115 //* DistanceToLine ----------------------------------------------------
 
  118 G4double G4VTwistSurface::DistanceToLine(const G4ThreeVector &p,
 
  119                                     const G4ThreeVector &x0,
 
  120                                     const G4ThreeVector &d,
 
  124    // Calculate distance to line,
 
  125    // then return distance and global intersection points.
 
  127    // p          - location of flying particle   
 
  128    // x0         - reference point of line
 
  129    // d          - direction vector of line
 
  130    // xx         - a foot of perpendicular line from p to the plane 
 
  131    // t          - distance from xx to p
 
  135    //    distance^2 = |(xx - p)|^2
 
  139    //   (d/dt)distance^2 = (d/dt)|((x0 + t*d) - p)|^2
 
  140    //                    = 2*t*|d|^2 + 2*d*(x0 - p)
 
  141    //                    = 0  // smallest distance
 
  143    //      t = - d*(x0 - p) / |d|^2
 
  147    G4ThreeVector dir = d.unit();
 
  148    t  = - dir * (x0 - p);      // |dir|^2 = 1.
 
  151    G4ThreeVector dist = xx - p;
 
  155 //=====================================================================
 
  156 //* IsAxis0 -----------------------------------------------------------
 
  159 G4bool G4VTwistSurface::IsAxis0(G4int areacode) const 
 
  161    if (areacode & sAxis0) return true;
 
  165 //=====================================================================
 
  166 //* IsAxis1 -----------------------------------------------------------
 
  169 G4bool G4VTwistSurface::IsAxis1(G4int areacode) const 
 
  171    if (areacode & sAxis1) return true;
 
  175 //=====================================================================
 
  176 //* IsOutside ---------------------------------------------------------
 
  179 G4bool G4VTwistSurface::IsOutside(G4int areacode) const 
 
  181    if (areacode & sInside) return false;
 
  185 //=====================================================================
 
  186 //* IsInside ----------------------------------------------------------
 
  189 G4bool G4VTwistSurface::IsInside(G4int areacode, G4bool testbitmode) const 
 
  191    if (areacode & sInside) {
 
  195        if (!((areacode & sBoundary) || (areacode & sCorner))) return true;
 
  201 //=====================================================================
 
  202 //* IsBoundary --------------------------------------------------------
 
  205 G4bool G4VTwistSurface::IsBoundary(G4int areacode, G4bool testbitmode) const 
 
  207    if ((areacode & sBoundary) == sBoundary) {
 
  211          if ((areacode & sInside) == sInside) return true;
 
  217 //=====================================================================
 
  218 //* IsCorner ----------------------------------------------------------
 
  221 G4bool G4VTwistSurface::IsCorner(G4int areacode, G4bool testbitmode) const 
 
  223    if ((areacode & sCorner) == sCorner) {
 
  227          if ((areacode & sInside) == sInside) return true;
 
  233 //=====================================================================
 
  234 //* GetAxisType -------------------------------------------------------
 
  237 G4int G4VTwistSurface::GetAxisType(G4int areacode, G4int whichaxis) const
 
  239    G4int axiscode = areacode & sAxisMask & whichaxis;
 
  241    if (axiscode == (sAxisX & sAxis0) ||
 
  242        axiscode == (sAxisX & sAxis1)) {
 
  244    } else if (axiscode == (sAxisY & sAxis0) ||
 
  245               axiscode == (sAxisY & sAxis1)) {
 
  247    } else if (axiscode == (sAxisZ & sAxis0) ||
 
  248               axiscode == (sAxisZ & sAxis1)) {
 
  250    } else if (axiscode == (sAxisRho & sAxis0) ||
 
  251               axiscode == (sAxisRho & sAxis1)) {
 
  253    } else if (axiscode == (sAxisPhi & sAxis0) ||
 
  254               axiscode == (sAxisPhi & sAxis1)) {
 
  257       std::ostringstream message;
 
  258       message << "Configuration not supported." << G4endl
 
  259               << "        areacode = " << areacode;
 
  260       G4Exception("G4VTwistSurface::GetAxisType()","GeomSolids0001",
 
  261                   FatalException, message);
 
  266 //=====================================================================
 
  267 //* ComputeGlobalPoint ------------------------------------------------
 
  270 G4ThreeVector G4VTwistSurface::ComputeGlobalPoint(const G4ThreeVector &lp) const
 
  272    return fRot * G4ThreeVector(lp) + fTrans;
 
  275 //=====================================================================
 
  276 //* ComputeGlobalPoint ------------------------------------------------
 
  279 G4ThreeVector G4VTwistSurface::ComputeLocalPoint(const G4ThreeVector &gp) const
 
  281    return fRot.inverse() * ( G4ThreeVector(gp) - fTrans ) ;
 
  284 //=====================================================================
 
  285 //* ComputeGlobalDirection --------------------------------------------
 
  288 G4ThreeVector G4VTwistSurface::ComputeGlobalDirection(const G4ThreeVector &lp) const
 
  290    return fRot * G4ThreeVector(lp); 
 
  293 //=====================================================================
 
  294 //* ComputeLocalDirection ---------------------------------------------
 
  297 G4ThreeVector G4VTwistSurface::ComputeLocalDirection(const G4ThreeVector &gp) const
 
  299    return fRot.inverse() * G4ThreeVector(gp);
 
  302 //=====================================================================
 
  303 //* SetNeighbours -----------------------------------------------------
 
  306 void G4VTwistSurface::SetNeighbours(G4VTwistSurface* axis0min, G4VTwistSurface* axis1min, 
 
  307                                G4VTwistSurface* axis0max, G4VTwistSurface* axis1max)
 
  309    fNeighbours[0] = axis0min;
 
  310    fNeighbours[1] = axis1min;
 
  311    fNeighbours[2] = axis0max;
 
  312    fNeighbours[3] = axis1max;
 
  315 //=====================================================================
 
  316 //* GetNeighbours -----------------------------------------------------
 
  319 G4int G4VTwistSurface::GetNeighbours(G4int areacode, G4VTwistSurface** surfaces) 
 
  322   int sAxis0Min = sAxis0 & sAxisMin ;
 
  323   int sAxis1Min = sAxis1 & sAxisMin ;
 
  324   int sAxis0Max = sAxis0 & sAxisMax ;
 
  325   int sAxis1Max = sAxis1 & sAxisMax ;
 
  329   if ( (areacode & sAxis0Min ) == sAxis0Min ) {
 
  330     surfaces[i] = fNeighbours[0] ;
 
  334   if ( ( areacode & sAxis1Min ) == sAxis1Min ) {
 
  335      surfaces[i] = fNeighbours[1] ;
 
  337     if ( i == 2 ) return i ;
 
  340   if ( ( areacode & sAxis0Max ) == sAxis0Max ) {
 
  341     surfaces[i] = fNeighbours[2] ;
 
  343     if ( i == 2 ) return i ;
 
  346   if ( ( areacode & sAxis1Max ) == sAxis1Max ) {
 
  347     surfaces[i] = fNeighbours[3] ;
 
  349     if ( i == 2 ) return i ;
 
  355 //=====================================================================
 
  356 //* GetCorner ---------------------------------------------------------
 
  359 G4ThreeVector G4VTwistSurface::GetCorner(G4int areacode) const
 
  361    if (!(areacode & sCorner)){
 
  362       std::ostringstream message;
 
  363       message << "Area code must represent corner." << G4endl
 
  364               << "        areacode = " << areacode;
 
  365       G4Exception("G4VTwistSurface::GetCorner()","GeomSolids0002",
 
  366                   FatalException, message);
 
  369    if ((areacode & sC0Min1Min) == sC0Min1Min) { 
 
  371    } else if ((areacode & sC0Max1Min) == sC0Max1Min) { 
 
  373    } else if ((areacode & sC0Max1Max) == sC0Max1Max) {
 
  375    } else if ((areacode & sC0Min1Max) == sC0Min1Max) { 
 
  378       std::ostringstream message;
 
  379       message << "Configuration not supported." << G4endl
 
  380               << "        areacode = " << areacode;
 
  381       G4Exception("G4VTwistSurface::GetCorner()", "GeomSolids0001",
 
  382                   FatalException, message);