//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: G4Sphere.cc,v 1.84 2009/08/07 15:56:23 gcosmo Exp $
// GEANT4 tag $Name: geant4-09-04-beta-01 $
//
// class G4Sphere
//
// Implementation for G4Sphere class
//
// History:
//
// 14.09.09 T.Nikitina: fix for phi section in DistanceToOut(p,v,..),as for G4Tubs,G4Cons 
// 26.03.09 G.Cosmo   : optimisations and uniform use of local radial tolerance
// 12.06.08 V.Grichine: fix for theta intersections in DistanceToOut(p,v,...)
// 22.07.05 O.Link    : Added check for intersection with double cone
// 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
// 16.09.04 V.Grichine: bug fixed in SurfaceNormal(p), theta normals
// 16.07.04 V.Grichine: bug fixed in DistanceToOut(p,v), Rmin go outside
// 02.06.04 V.Grichine: bug fixed in DistanceToIn(p,v), on Rmax,Rmin go inside
// 30.10.03 J.Apostolakis: new algorithm in Inside for SPhi-sections
// 29.10.03 J.Apostolakis: fix in Inside for SPhi-0.5*kAngTol < phi<SPhi, SPhi<0
// 19.06.02 V.Grichine: bug fixed in Inside(p), && -> && fDTheta - kAngTolerance
// 30.01.02 V.Grichine: bug fixed in Inside(p), && -> || at l.451
// 06.03.00 V.Grichine: modifications in Distance ToOut(p,v,...)
// 18.11.99 V.Grichine: side = kNull in Distance ToOut(p,v,...)
// 25.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), phi intersections
// 12.11.98 V.Grichine: bug fixed in DistanceToIn(p,v), theta intersections
// 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...)
// 17.09.96 V.Grichine: final modifications to commit
// 28.03.94 P.Kent: old C++ code converted to tolerant geometry
// --------------------------------------------------------------------

#include "G4Sphere.hh"

#include "G4VoxelLimits.hh"
#include "G4AffineTransform.hh"
#include "G4GeometryTolerance.hh"

#include "G4VPVParameterisation.hh"

#include "Randomize.hh"

#include "meshdefs.hh"

#include "G4VGraphicsScene.hh"
#include "G4VisExtent.hh"
#include "G4Polyhedron.hh"
#include "G4NURBS.hh"
#include "G4NURBSbox.hh"

using namespace CLHEP;

// Private enum: Not for external use - used by distanceToOut

enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kSTheta,kETheta};

// used by normal

enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNSTheta,kNETheta};

////////////////////////////////////////////////////////////////////////
//
// constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
//             - note if pDPhi>2PI then reset to 2PI

G4Sphere::G4Sphere( const G4String& pName,
                          G4double pRmin, G4double pRmax,
                          G4double pSPhi, G4double pDPhi,
                          G4double pSTheta, G4double pDTheta )
  : G4CSGSolid(pName), fFullPhiSphere(true), fFullThetaSphere(true)
{
  fEpsilon = 2.0e-11;  // relative radial tolerance constant

  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();

  // Check radii and set radial tolerances

  G4double kRadTolerance = G4GeometryTolerance::GetInstance()
                         ->GetRadialTolerance();
  if ( (pRmin < pRmax) && (pRmax >= 10*kRadTolerance) && (pRmin >= 0) )
  {
    fRmin=pRmin; fRmax=pRmax;
    fRminTolerance = (pRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
    fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
  }
  else
  {
    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
           << "        Invalide values for radii ! - "
           << "        pRmin = " << pRmin << ", pRmax = " << pRmax << G4endl;
    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
                "Invalid radii");
  }

  // Check angles

  CheckPhiAngles(pSPhi, pDPhi);
  CheckThetaAngles(pSTheta, pDTheta);
}

///////////////////////////////////////////////////////////////////////
//
// Fake default constructor - sets only member data and allocates memory
//                            for usage restricted to object persistency.
//
G4Sphere::G4Sphere( __void__& a )
  : G4CSGSolid(a)
{
}

/////////////////////////////////////////////////////////////////////
//
// Destructor

G4Sphere::~G4Sphere()
{
}

//////////////////////////////////////////////////////////////////////////
//
// Dispatch to parameterisation for replication mechanism dimension
// computation & modification.

void G4Sphere::ComputeDimensions(       G4VPVParameterisation* p,
                                  const G4int n,
                                  const G4VPhysicalVolume* pRep)
{
  p->ComputeDimensions(*this,n,pRep);
}

////////////////////////////////////////////////////////////////////////////
//
// Calculate extent under transform and specified limit

G4bool G4Sphere::CalculateExtent( const EAxis pAxis,
                                  const G4VoxelLimits& pVoxelLimit,
                                  const G4AffineTransform& pTransform,
                                        G4double& pMin, G4double& pMax ) const
{
  if ( fFullSphere )
  {
    // Special case handling for solid spheres-shells
    // (rotation doesn't influence).
    // Compute x/y/z mins and maxs for bounding box respecting limits,
    // with early returns if outside limits. Then switch() on pAxis,
    // and compute exact x and y limit for x/y case
      
    G4double xoffset,xMin,xMax;
    G4double yoffset,yMin,yMax;
    G4double zoffset,zMin,zMax;

    G4double diff1,diff2,maxDiff,newMin,newMax;
    G4double xoff1,xoff2,yoff1,yoff2;

    xoffset=pTransform.NetTranslation().x();
    xMin=xoffset-fRmax;
    xMax=xoffset+fRmax;
    if (pVoxelLimit.IsXLimited())
    {
      if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
        || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
      {
        return false;
      }
      else
      {
        if (xMin<pVoxelLimit.GetMinXExtent())
        {
          xMin=pVoxelLimit.GetMinXExtent();
        }
        if (xMax>pVoxelLimit.GetMaxXExtent())
        {
          xMax=pVoxelLimit.GetMaxXExtent();
        }
      }
    }

    yoffset=pTransform.NetTranslation().y();
    yMin=yoffset-fRmax;
    yMax=yoffset+fRmax;
    if (pVoxelLimit.IsYLimited())
    {
      if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
        || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
      {
        return false;
      }
      else
      {
        if (yMin<pVoxelLimit.GetMinYExtent())
        {
          yMin=pVoxelLimit.GetMinYExtent();
        }
        if (yMax>pVoxelLimit.GetMaxYExtent())
        {
          yMax=pVoxelLimit.GetMaxYExtent();
        }
      }
    }

    zoffset=pTransform.NetTranslation().z();
    zMin=zoffset-fRmax;
    zMax=zoffset+fRmax;
    if (pVoxelLimit.IsZLimited())
    {
      if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
        || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
      {
        return false;
      }
      else
      {
        if (zMin<pVoxelLimit.GetMinZExtent())
        {
          zMin=pVoxelLimit.GetMinZExtent();
        }
        if (zMax>pVoxelLimit.GetMaxZExtent())
        {
          zMax=pVoxelLimit.GetMaxZExtent();
        }
      }
    }

    // Known to cut sphere

    switch (pAxis)
    {
      case kXAxis:
        yoff1=yoffset-yMin;
        yoff2=yMax-yoffset;
        if ((yoff1>=0) && (yoff2>=0))
        {
          // Y limits cross max/min x => no change
          //
          pMin=xMin;
          pMax=xMax;
        }
        else
        {
          // Y limits don't cross max/min x => compute max delta x,
          // hence new mins/maxs
          //
          diff1=std::sqrt(fRmax*fRmax-yoff1*yoff1);
          diff2=std::sqrt(fRmax*fRmax-yoff2*yoff2);
          maxDiff=(diff1>diff2) ? diff1:diff2;
          newMin=xoffset-maxDiff;
          newMax=xoffset+maxDiff;
          pMin=(newMin<xMin) ? xMin : newMin;
          pMax=(newMax>xMax) ? xMax : newMax;
        }
        break;
      case kYAxis:
        xoff1=xoffset-xMin;
        xoff2=xMax-xoffset;
        if ((xoff1>=0) && (xoff2>=0))
        {
          // X limits cross max/min y => no change
          //
          pMin=yMin;
          pMax=yMax;
        }
        else
        {
          // X limits don't cross max/min y => compute max delta y,
          // hence new mins/maxs
          //
          diff1=std::sqrt(fRmax*fRmax-xoff1*xoff1);
          diff2=std::sqrt(fRmax*fRmax-xoff2*xoff2);
          maxDiff=(diff1>diff2) ? diff1:diff2;
          newMin=yoffset-maxDiff;
          newMax=yoffset+maxDiff;
          pMin=(newMin<yMin) ? yMin : newMin;
          pMax=(newMax>yMax) ? yMax : newMax;
        }
        break;
      case kZAxis:
        pMin=zMin;
        pMax=zMax;
        break;
      default:
        break;
    }
    pMin-=kCarTolerance;
    pMax+=kCarTolerance;

    return true;  
  }
  else       // Transformed cutted sphere
  {
    G4int i,j,noEntries,noBetweenSections;
    G4bool existsAfterClip=false;

    // Calculate rotated vertex coordinates

    G4ThreeVectorList* vertices;
    G4int  noPolygonVertices ;
    vertices=CreateRotatedVertices(pTransform,noPolygonVertices);

    pMin=+kInfinity;
    pMax=-kInfinity;

    noEntries=vertices->size();  // noPolygonVertices*noPhiCrossSections
    noBetweenSections=noEntries-noPolygonVertices;

    G4ThreeVectorList ThetaPolygon ;
    for (i=0;i<noEntries;i+=noPolygonVertices)
    {
      for(j=0;j<(noPolygonVertices/2)-1;j++)
      {
        ThetaPolygon.push_back((*vertices)[i+j]) ;      
        ThetaPolygon.push_back((*vertices)[i+j+1]) ;      
        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ;      
        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ;      
        CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
        ThetaPolygon.clear() ;
      }
    }
    for (i=0;i<noBetweenSections;i+=noPolygonVertices)
    {
      for(j=0;j<noPolygonVertices-1;j++)
      {
        ThetaPolygon.push_back((*vertices)[i+j]) ;      
        ThetaPolygon.push_back((*vertices)[i+j+1]) ;      
        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ;      
        ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ;      
        CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
        ThetaPolygon.clear() ;
      }
      ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ;      
      ThetaPolygon.push_back((*vertices)[i]) ;  
      ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ;      
      ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ;      
      CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
      ThetaPolygon.clear() ;
    }
      
    if ((pMin!=kInfinity) || (pMax!=-kInfinity))
    {
      existsAfterClip=true;

      // Add 2*tolerance to avoid precision troubles
      //
      pMin-=kCarTolerance;
      pMax+=kCarTolerance;
    }
    else
    {
      // Check for case where completely enveloping clipping volume
      // If point inside then we are confident that the solid completely
      // envelopes the clipping volume. Hence set min/max extents according
      // to clipping volume extents along the specified axis.

      G4ThreeVector clipCentre(
          (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
          (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
          (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
        
      if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
      {
        existsAfterClip=true;
        pMin=pVoxelLimit.GetMinExtent(pAxis);
        pMax=pVoxelLimit.GetMaxExtent(pAxis);
      }
    }
    delete vertices;
    return existsAfterClip;
  }
}

///////////////////////////////////////////////////////////////////////////
//
// Return whether point inside/outside/on surface
// Split into radius, phi, theta checks
// Each check modifies 'in', or returns as approprate

EInside G4Sphere::Inside( const G4ThreeVector& p ) const
{
  G4double rho,rho2,rad2,tolRMin,tolRMax;
  G4double pPhi,pTheta;
  EInside in = kOutside;
  static const G4double halfAngTolerance = kAngTolerance*0.5;
  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
  const G4double halfRminTolerance = fRminTolerance*0.5;
  const G4double Rmax_minus = fRmax - halfRmaxTolerance;
  const G4double Rmin_plus  = (fRmin > 0) ? fRmin+halfRminTolerance : 0;

  rho2 = p.x()*p.x() + p.y()*p.y() ;
  rad2 = rho2 + p.z()*p.z() ;

  // Check radial surfaces. Sets 'in'

  tolRMin = Rmin_plus;
  tolRMax = Rmax_minus;

  if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
  {
    in = kInside;
  }
  else
  {
    tolRMax = fRmax + halfRmaxTolerance;                  // outside case
    tolRMin = std::max(fRmin-halfRminTolerance, 0.);      // outside case
    if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
    {
      in = kSurface;
    }
    else
    {
      return in = kOutside;
    }
  }

  // Phi boundaries   : Do not check if it has no phi boundary!

  if ( !fFullPhiSphere && rho2 )  // [fDPhi < twopi] and [p.x or p.y]
  {
    pPhi = std::atan2(p.y(),p.x()) ;

    if      ( pPhi < fSPhi - halfAngTolerance  ) { pPhi += twopi; }
    else if ( pPhi > ePhi + halfAngTolerance )   { pPhi -= twopi; }
    
    if ( (pPhi < fSPhi - halfAngTolerance)
      || (pPhi > ePhi + halfAngTolerance) )      { return in = kOutside; }
    
    else if (in == kInside)  // else it's kSurface anyway already
    {
      if ( (pPhi < fSPhi + halfAngTolerance)
        || (pPhi > ePhi - halfAngTolerance) )    { in = kSurface; }     
    }
  }

  // Theta bondaries
  
  if ( (rho2 || p.z()) && (!fFullThetaSphere) )
  {
    rho    = std::sqrt(rho2);
    pTheta = std::atan2(rho,p.z());

    if ( in == kInside )
    {
      if ( (pTheta < fSTheta + halfAngTolerance)
        || (pTheta > eTheta - halfAngTolerance) )
      {
        if ( (pTheta >= fSTheta - halfAngTolerance)
          && (pTheta <= eTheta + halfAngTolerance) )
        {
          in = kSurface;
        }
        else
        {
          in = kOutside;
        }
      }
    }
    else
    {
      if ( (pTheta < fSTheta - halfAngTolerance)
        || (pTheta > eTheta + halfAngTolerance) )
      {
        in = kOutside;
      }
    }
  }
  return in;
}

/////////////////////////////////////////////////////////////////////
//
// Return unit normal of surface closest to p
// - note if point on z axis, ignore phi divided sides
// - unsafe if point close to z axis a rmin=0 - no explicit checks

G4ThreeVector G4Sphere::SurfaceNormal( const G4ThreeVector& p ) const
{
  G4int noSurfaces = 0;  
  G4double rho, rho2, rad, pTheta, pPhi=0.;
  G4double distRMin = kInfinity;
  G4double distSPhi = kInfinity, distEPhi = kInfinity;
  G4double distSTheta = kInfinity, distETheta = kInfinity;
  G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
  G4ThreeVector norm, sumnorm(0.,0.,0.);

  static const G4double halfCarTolerance = 0.5*kCarTolerance;
  static const G4double halfAngTolerance = 0.5*kAngTolerance;

  rho2 = p.x()*p.x()+p.y()*p.y();
  rad  = std::sqrt(rho2+p.z()*p.z());
  rho  = std::sqrt(rho2);

  G4double    distRMax = std::fabs(rad-fRmax);
  if (fRmin)  distRMin = std::fabs(rad-fRmin);
    
  if ( rho && !fFullSphere )
  {
    pPhi = std::atan2(p.y(),p.x());

    if (pPhi < fSPhi-halfAngTolerance)     { pPhi += twopi; }
    else if (pPhi > ePhi+halfAngTolerance) { pPhi -= twopi; }
  }
  if ( !fFullPhiSphere )
  {
    if ( rho )
    {
      distSPhi = std::fabs( pPhi-fSPhi ); 
      distEPhi = std::fabs( pPhi-ePhi ); 
    }
    else if( !fRmin )
    {
      distSPhi = 0.; 
      distEPhi = 0.; 
    }
    nPs = G4ThreeVector(sinSPhi,-cosSPhi,0);
    nPe = G4ThreeVector(-sinEPhi,cosEPhi,0);
  }        
  if ( !fFullThetaSphere )
  {
    if ( rho )
    {
      pTheta     = std::atan2(rho,p.z());
      distSTheta = std::fabs(pTheta-fSTheta); 
      distETheta = std::fabs(pTheta-eTheta);
 
      nTs = G4ThreeVector(-cosSTheta*p.x()/rho,
                          -cosSTheta*p.y()/rho,
                           sinSTheta          );

      nTe = G4ThreeVector( cosETheta*p.x()/rho,
                           cosETheta*p.y()/rho,
                          -sinETheta          );    
    }
    else if( !fRmin )
    {
      if ( fSTheta )  
      {              
        distSTheta = 0.;
        nTs = G4ThreeVector(0.,0.,-1.);
      }
      if ( eTheta < pi )
      {              
        distETheta = 0.;
        nTe = G4ThreeVector(0.,0.,1.);
      }
    }    
  }
  if( rad )  { nR = G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad); }

  if( distRMax <= halfCarTolerance )
  {
    noSurfaces ++;
    sumnorm += nR;
  }
  if( fRmin && (distRMin <= halfCarTolerance) )
  {
    noSurfaces ++;
    sumnorm -= nR;
  }
  if( !fFullPhiSphere )   
  {
    if (distSPhi <= halfAngTolerance)
    {
      noSurfaces ++;
      sumnorm += nPs;
    }
    if (distEPhi <= halfAngTolerance) 
    {
      noSurfaces ++;
      sumnorm += nPe;
    }
  }
  if ( !fFullThetaSphere )
  {
    if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
    {
      noSurfaces ++;
      if ((rad <= halfCarTolerance) && fFullPhiSphere)  { sumnorm += nZ;  }
      else                                              { sumnorm += nTs; }
    }
    if ((distETheta <= halfAngTolerance) && (eTheta < pi)) 
    {
      noSurfaces ++;
      if ((rad <= halfCarTolerance) && fFullPhiSphere)  { sumnorm -= nZ;  }
      else                                              { sumnorm += nTe; }
      if(sumnorm.z() == 0.)  { sumnorm += nZ; }
    }
  }
  if ( noSurfaces == 0 )
  {
#ifdef G4CSGDEBUG
    G4Exception("G4Sphere::SurfaceNormal(p)", "Notification", JustWarning, 
                "Point p is not on surface !?" ); 
#endif
     norm = ApproxSurfaceNormal(p);
  }
  else if ( noSurfaces == 1 )  { norm = sumnorm; }
  else                         { norm = sumnorm.unit(); }
  return norm;
}


/////////////////////////////////////////////////////////////////////////////////////////////
//
// Algorithm for SurfaceNormal() following the original specification
// for points not on the surface

G4ThreeVector G4Sphere::ApproxSurfaceNormal( const G4ThreeVector& p ) const
{
  ENorm side;
  G4ThreeVector norm;
  G4double rho,rho2,rad,pPhi,pTheta;
  G4double distRMin,distRMax,distSPhi,distEPhi,
           distSTheta,distETheta,distMin;

  rho2=p.x()*p.x()+p.y()*p.y();
  rad=std::sqrt(rho2+p.z()*p.z());
  rho=std::sqrt(rho2);

  //
  // Distance to r shells
  //

  distRMax=std::fabs(rad-fRmax);
  if (fRmin)
  {
    distRMin=std::fabs(rad-fRmin);
      
    if (distRMin<distRMax)
    {
      distMin=distRMin;
      side=kNRMin;
    }
    else
    {   
      distMin=distRMax;
      side=kNRMax;
    }
  }
  else
  {
    distMin=distRMax;
    side=kNRMax;
  }

  //
  // Distance to phi planes
  //
  // Protected against (0,0,z) 
    
  pPhi = std::atan2(p.y(),p.x());
  if (pPhi<0) { pPhi += twopi; }

  if (!fFullPhiSphere && rho)
  {
    if (fSPhi<0)
    {
      distSPhi=std::fabs(pPhi-(fSPhi+twopi))*rho;
    }
    else
    {
      distSPhi=std::fabs(pPhi-fSPhi)*rho;
    }

    distEPhi=std::fabs(pPhi-fSPhi-fDPhi)*rho;

    // Find new minimum
    //
    if (distSPhi<distEPhi)
    {
      if (distSPhi<distMin)
      {
        distMin=distSPhi;
        side=kNSPhi;
      }
    }
    else
    {
      if (distEPhi<distMin)
      {
        distMin=distEPhi;
        side=kNEPhi;
      }
    }
  }

  //
  // Distance to theta planes
  //

  if (!fFullThetaSphere && rad)
  {
    pTheta=std::atan2(rho,p.z());
    distSTheta=std::fabs(pTheta-fSTheta)*rad;
    distETheta=std::fabs(pTheta-fSTheta-fDTheta)*rad;

    // Find new minimum
    //
    if (distSTheta<distETheta)
    {
      if (distSTheta<distMin)
      {
        distMin = distSTheta ;
        side = kNSTheta ;
      }
    }
    else
    {
      if (distETheta<distMin)
      {
        distMin = distETheta ;
        side = kNETheta ;
      }
    }
  }

  switch (side)
  {
    case kNRMin:      // Inner radius
      norm=G4ThreeVector(-p.x()/rad,-p.y()/rad,-p.z()/rad);
      break;
    case kNRMax:      // Outer radius
      norm=G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad);
      break;
    case kNSPhi:
      norm=G4ThreeVector(sinSPhi,-cosSPhi,0);
      break;
    case kNEPhi:
      norm=G4ThreeVector(-sinEPhi,cosEPhi,0);
      break;
    case kNSTheta:
      norm=G4ThreeVector(-cosSTheta*std::cos(pPhi),
                         -cosSTheta*std::sin(pPhi),
                          sinSTheta            );
      break;
    case kNETheta:
      norm=G4ThreeVector( cosETheta*std::cos(pPhi),
                          cosETheta*std::sin(pPhi),
                         -sinETheta              );
      break;
    default:
      DumpInfo();
      G4Exception("G4Sphere::ApproxSurfaceNormal()","Notification",JustWarning,
                  "Undefined side for valid surface normal to solid.");
      break;    
  }

  return norm;
}

///////////////////////////////////////////////////////////////////////////////
//
// Calculate distance to shape from outside, along normalised vector
// - return kInfinity if no intersection, or intersection distance <= tolerance
//
// -> If point is outside outer radius, compute intersection with rmax
//        - if no intersection return
//        - if  valid phi,theta return intersection Dist
//
// -> If shell, compute intersection with inner radius, taking largest +ve root
//        - if valid phi,theta, save intersection
//
// -> If phi segmented, compute intersection with phi half planes
//        - if valid intersection(r,theta), return smallest intersection of
//          inner shell & phi intersection
//
// -> If theta segmented, compute intersection with theta cones
//        - if valid intersection(r,phi), return smallest intersection of
//          inner shell & theta intersection
//
//
// NOTE:
// - `if valid' (above) implies tolerant checking of intersection points
//
// OPT:
// Move tolIO/ORmin/RMax2 precalcs to where they are needed -
// not required for most cases.
// Avoid atan2 for non theta cut G4Sphere.

G4double G4Sphere::DistanceToIn( const G4ThreeVector& p,
                                 const G4ThreeVector& v  ) const
{
  G4double snxt = kInfinity ;      // snxt = default return value
  G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
  G4double tolSTheta=0., tolETheta=0. ;
  const G4double dRmax = 100.*fRmax;

  static const G4double halfCarTolerance = kCarTolerance*0.5;
  static const G4double halfAngTolerance = kAngTolerance*0.5;
  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
  const G4double halfRminTolerance = fRminTolerance*0.5;
  const G4double tolORMin2 = (fRmin>halfRminTolerance)
               ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
  const G4double tolIRMin2 =
               (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
  const G4double tolORMax2 =
               (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
  const G4double tolIRMax2 =
               (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);

  // Intersection point
  //
  G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;

  // Phi intersection
  //
  G4double Comp ; 

  // Phi precalcs
  //
  G4double Dist, cosPsi ;

  // Theta precalcs
  //
  G4double dist2STheta, dist2ETheta ;
  G4double t1, t2, b, c, d2, d, s = kInfinity ;

  // General Precalcs
  //
  rho2 = p.x()*p.x() + p.y()*p.y() ;
  rad2 = rho2 + p.z()*p.z() ;
  pTheta = std::atan2(std::sqrt(rho2),p.z()) ;

  pDotV2d = p.x()*v.x() + p.y()*v.y() ;
  pDotV3d = pDotV2d + p.z()*v.z() ;

  // Theta precalcs
  //
  if (!fFullThetaSphere)
  {
    tolSTheta = fSTheta - halfAngTolerance ;
    tolETheta = eTheta + halfAngTolerance ;
  }

  // Outer spherical shell intersection
  // - Only if outside tolerant fRmax
  // - Check for if inside and outer G4Sphere heading through solid (-> 0)
  // - No intersect -> no intersection with G4Sphere
  //
  // Shell eqn: x^2+y^2+z^2=RSPH^2
  //
  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
  //
  // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
  // =>      rad2        +2s(pDotV3d)       +s^2                =R^2
  //
  // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))

  c = rad2 - fRmax*fRmax ;

  if (c > fRmaxTolerance*fRmax)
  {
    // If outside tolerant boundary of outer G4Sphere
    // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance]

    d2 = pDotV3d*pDotV3d - c ;

    if ( d2 >= 0 )
    {
      s = -pDotV3d - std::sqrt(d2) ;

      if (s >= 0 )
      {
        if ( s>dRmax ) // Avoid rounding errors due to precision issues seen on
        {              // 64 bits systems. Split long distances and recompute
          G4double fTerm = s-std::fmod(s,dRmax);
          s = fTerm + DistanceToIn(p+fTerm*v,v);
        } 
        xi   = p.x() + s*v.x() ;
        yi   = p.y() + s*v.y() ;
        rhoi = std::sqrt(xi*xi + yi*yi) ;

        if (!fFullPhiSphere && rhoi)    // Check phi intersection
        {
          cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;

          if (cosPsi >= cosHDPhiOT)
          {
            if (!fFullThetaSphere)   // Check theta intersection
            {
              zi = p.z() + s*v.z() ;

              // rhoi & zi can never both be 0
              // (=>intersect at origin =>fRmax=0)
              //
              iTheta = std::atan2(rhoi,zi) ;
              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
              {
                return snxt = s ;
              }
            }
            else
            {
              return snxt=s;
            }
          }
        }
        else
        {
          if (!fFullThetaSphere)    // Check theta intersection
          {
            zi = p.z() + s*v.z() ;

            // rhoi & zi can never both be 0
            // (=>intersect at origin => fRmax=0 !)
            //
            iTheta = std::atan2(rhoi,zi) ;
            if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
            {
              return snxt=s;
            }
          }
          else
          {
            return snxt = s ;
          }
        }          
      }
    }
    else    // No intersection with G4Sphere
    {
      return snxt=kInfinity;
    }
  }
  else
  {
    // Inside outer radius
    // check not inside, and heading through G4Sphere (-> 0 to in)

    d2 = pDotV3d*pDotV3d - c ;

    if ( (rad2 > tolIRMax2)
      && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
    {
      if (!fFullPhiSphere)
      {
        // Use inner phi tolerant boundary -> if on tolerant
        // phi boundaries, phi intersect code handles leaving/entering checks

        cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;

        if (cosPsi>=cosHDPhiIT)
        { 
          // inside radii, delta r -ve, inside phi

          if ( !fFullThetaSphere )
          {
            if ( (pTheta >= tolSTheta + kAngTolerance)
              && (pTheta <= tolETheta - kAngTolerance) )
            {
              return snxt=0;
            }
          }
          else    // strictly inside Theta in both cases
          {
            return snxt=0;
          }
        }
      }
      else
      {
        if ( !fFullThetaSphere )
        {
          if ( (pTheta >= tolSTheta + kAngTolerance)
            && (pTheta <= tolETheta - kAngTolerance) )
          {
            return snxt=0;
          }
        }
        else   // strictly inside Theta in both cases
        {
          return snxt=0;
        }
      }
    }
  }

  // Inner spherical shell intersection
  // - Always farthest root, because would have passed through outer
  //   surface first.
  // - Tolerant check if travelling through solid

  if (fRmin)
  {
    c  = rad2 - fRmin*fRmin ;
    d2 = pDotV3d*pDotV3d - c ;

    // Within tolerance inner radius of inner G4Sphere
    // Check for immediate entry/already inside and travelling outwards

    if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
      && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
    {
      if ( !fFullPhiSphere )
      {
        // Use inner phi tolerant boundary -> if on tolerant
        // phi boundaries, phi intersect code handles leaving/entering checks

        cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ;
        if (cosPsi >= cosHDPhiIT)
        { 
          // inside radii, delta r -ve, inside phi
          //
          if ( !fFullThetaSphere )
          {
            if ( (pTheta >= tolSTheta + kAngTolerance)
              && (pTheta <= tolETheta - kAngTolerance) )
            {
              return snxt=0;
            }
          }
          else
          {
            return snxt = 0 ;
          }
        }
      }
      else
      {
        if ( !fFullThetaSphere )
        {
          if ( (pTheta >= tolSTheta + kAngTolerance)
            && (pTheta <= tolETheta - kAngTolerance) )
          {
            return snxt = 0 ;
          }
        }
        else
        {
          return snxt=0;
        }
      }
    }
    else   // Not special tolerant case
    {
      if (d2 >= 0)
      {
        s = -pDotV3d + std::sqrt(d2) ;
        if ( s >= halfRminTolerance )  // It was >= 0 ??
        {
          xi   = p.x() + s*v.x() ;
          yi   = p.y() + s*v.y() ;
          rhoi = std::sqrt(xi*xi+yi*yi) ;

          if ( !fFullPhiSphere && rhoi )   // Check phi intersection
          {
            cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;

            if (cosPsi >= cosHDPhiOT)
            {
              if ( !fFullThetaSphere )  // Check theta intersection
              {
                zi = p.z() + s*v.z() ;

                // rhoi & zi can never both be 0
                // (=>intersect at origin =>fRmax=0)
                //
                iTheta = std::atan2(rhoi,zi) ;
                if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
                {
                  snxt = s ;
                }
              }
              else
              {
                snxt=s;
              }
            }
          }
          else
          {
            if ( !fFullThetaSphere )   // Check theta intersection
            {
              zi = p.z() + s*v.z() ;

              // rhoi & zi can never both be 0
              // (=>intersect at origin => fRmax=0 !)
              //
              iTheta = std::atan2(rhoi,zi) ;
              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
              {
                snxt = s;
              }
            }
            else
            {
              snxt = s;
            }
          }
        }
      }
    }
  }

  // Phi segment intersection
  //
  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
  //
  // o NOTE: Large duplication of code between sphi & ephi checks
  //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
  //            intersection check <=0 -> >=0
  //         -> Should use some form of loop Construct
  //
  if ( !fFullPhiSphere )
  {
    // First phi surface ('S'tarting phi)
    // Comp = Component in outwards normal dirn
    //
    Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
                    
    if ( Comp < 0 )
    {
      Dist = p.y()*cosSPhi - p.x()*sinSPhi ;

      if (Dist < halfCarTolerance)
      {
        s = Dist/Comp ;

        if (s < snxt)
        {
          if ( s > 0 )
          {
            xi    = p.x() + s*v.x() ;
            yi    = p.y() + s*v.y() ;
            zi    = p.z() + s*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;
          }
          else
          {
            s     = 0     ;
            xi    = p.x() ;
            yi    = p.y() ;
            zi    = p.z() ;
            rhoi2 = rho2  ;
            radi2 = rad2  ;
          }
          if ( (radi2 <= tolORMax2)
            && (radi2 >= tolORMin2)
            && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
          {
            // Check theta intersection
            // rhoi & zi can never both be 0
            // (=>intersect at origin =>fRmax=0)
            //
            if ( !fFullThetaSphere )
            {
              iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
              {
                // r and theta intersections good
                // - check intersecting with correct half-plane

                if ((yi*cosCPhi-xi*sinCPhi) <= 0)
                {
                  snxt = s ;
                }
              }
            }
            else
            {
              snxt = s ;
            }
          }
        }
      }
    }

    // Second phi surface ('E'nding phi)
    // Component in outwards normal dirn

    Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
        
    if (Comp < 0)
    {
      Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
      if ( Dist < halfCarTolerance )
      {
        s = Dist/Comp ;

        if ( s < snxt )
        {
          if (s > 0)
          {
            xi    = p.x() + s*v.x() ;
            yi    = p.y() + s*v.y() ;
            zi    = p.z() + s*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;
          }
          else
          {
            s     = 0     ;
            xi    = p.x() ;
            yi    = p.y() ;
            zi    = p.z() ;
            rhoi2 = rho2  ;
            radi2 = rad2  ;
          }
          if ( (radi2 <= tolORMax2)
            && (radi2 >= tolORMin2)
            && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
          {
            // Check theta intersection
            // rhoi & zi can never both be 0
            // (=>intersect at origin =>fRmax=0)
            //
            if ( !fFullThetaSphere )
            {
              iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
              if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
              {
                // r and theta intersections good
                // - check intersecting with correct half-plane

                if ((yi*cosCPhi-xi*sinCPhi) >= 0)
                {
                  snxt = s ;
                }
              }
            }
            else
            {
              snxt = s ;
            }
          }
        }
      }
    }
  }

  // Theta segment intersection

  if ( !fFullThetaSphere )
  {

    // Intersection with theta surfaces
    // Known failure cases:
    // o  Inside tolerance of stheta surface, skim
    //    ~parallel to cone and Hit & enter etheta surface [& visa versa]
    //
    //    To solve: Check 2nd root of etheta surface in addition to stheta
    //
    // o  start/end theta is exactly pi/2 
    // Intersections with cones
    //
    // Cone equation: x^2+y^2=z^2tan^2(t)
    //
    // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
    //
    // => (px^2+py^2-pz^2tan^2(t))+2s(pxvx+pyvy-pzvztan^2(t))
    //       + s^2(vx^2+vy^2-vz^2tan^2(t)) = 0
    //
    // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0

    if (fSTheta)
    {
      dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
    }
    else
    {
      dist2STheta = kInfinity ;
    }
    if ( eTheta < pi )
    {
      dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
    }
    else
    {
      dist2ETheta=kInfinity;
    }      
    if ( pTheta < tolSTheta )
    {
      // Inside (theta<stheta-tol) s theta cone
      // First root of stheta cone, second if first root -ve

      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
      if (t1)
      {   
        b  = t2/t1 ;
        c  = dist2STheta/t1 ;
        d2 = b*b - c ;

        if ( d2 >= 0 )
        {
          d = std::sqrt(d2) ;
          s = -b - d ;    // First root
          zi    = p.z() + s*v.z();

          if ( (s < 0) || (zi*(fSTheta - halfpi) > 0) )
          {
            s = -b+d;    // Second root
          }
          if ((s >= 0) && (s < snxt))
          {
            xi    = p.x() + s*v.x();
            yi    = p.y() + s*v.y();
            zi    = p.z() + s*v.z();
            rhoi2 = xi*xi + yi*yi;
            radi2 = rhoi2 + zi*zi;
            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(fSTheta - halfpi) <= 0) )
            {
              if ( !fFullPhiSphere && rhoi2 )  // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if (cosPsi >= cosHDPhiOT)
                {
                  snxt = s ;
                }
              }
              else
              {
                snxt = s ;
              }
            }
          }
        }
      }

      // Possible intersection with ETheta cone. 
      // Second >= 0 root should be considered
        
      if ( eTheta < pi )
      {
        t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
        t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
        if (t1)
        { 
          b  = t2/t1 ;
          c  = dist2ETheta/t1 ;
          d2 = b*b - c ;

          if (d2 >= 0)
          {
            d = std::sqrt(d2) ;
            s = -b + d ;    // Second root

            if ( (s >= 0) && (s < snxt) )
            {
              xi    = p.x() + s*v.x() ;
              yi    = p.y() + s*v.y() ;
              zi    = p.z() + s*v.z() ;
              rhoi2 = xi*xi + yi*yi   ;
              radi2 = rhoi2 + zi*zi   ;

              if ( (radi2 <= tolORMax2)
                && (radi2 >= tolORMin2)
                && (zi*(eTheta - halfpi) <= 0) )
              {
                if (!fFullPhiSphere && rhoi2)   // Check phi intersection
                {
                  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                  if (cosPsi >= cosHDPhiOT)
                  {
                    snxt = s ;
                  }
                }
                else
                {
                  snxt = s ;
                }
              }
            }
          }
        }
      }
    }  
    else if ( pTheta > tolETheta ) 
    { 
      // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0)
      // Inside (theta > etheta+tol) e-theta cone
      // First root of etheta cone, second if first root 'imaginary'

      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
      if (t1)
      {  
        b  = t2/t1 ;
        c  = dist2ETheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d = std::sqrt(d2) ;
          s = -b - d ;    // First root
          zi    = p.z() + s*v.z();

          if ( (s < 0) || (zi*(eTheta - halfpi) > 0) )
          {
            s = -b + d ;           // second root
          }
          if ( (s >= 0) && (s < snxt) )
          {
            xi    = p.x() + s*v.x() ;
            yi    = p.y() + s*v.y() ;
            zi    = p.z() + s*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2) 
              && (zi*(eTheta - halfpi) <= 0) )
            {
              if (!fFullPhiSphere && rhoi2)  // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if (cosPsi >= cosHDPhiOT)
                {
                  snxt = s ;
                }
              }
              else
              {
                snxt = s ;
              }
            }
          }
        }
      }

      // Possible intersection with STheta cone. 
      // Second >= 0 root should be considered
        
      if ( fSTheta )
      {
        t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
        t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
        if (t1)
        { 
          b  = t2/t1 ;
          c  = dist2STheta/t1 ;
          d2 = b*b - c ;

          if (d2 >= 0)
          {
            d = std::sqrt(d2) ;
            s = -b + d ;    // Second root

            if ( (s >= 0) && (s < snxt) )
            {
              xi    = p.x() + s*v.x() ;
              yi    = p.y() + s*v.y() ;
              zi    = p.z() + s*v.z() ;
              rhoi2 = xi*xi + yi*yi   ;
              radi2 = rhoi2 + zi*zi   ;

              if ( (radi2 <= tolORMax2)
                && (radi2 >= tolORMin2)
                && (zi*(fSTheta - halfpi) <= 0) )
              {
                if (!fFullPhiSphere && rhoi2)   // Check phi intersection
                {
                  cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                  if (cosPsi >= cosHDPhiOT)
                  {
                    snxt = s ;
                  }
                }
                else
                {
                  snxt = s ;
                }
              }
            }
          }
        }
      }  
    }     
    else if ( (pTheta < tolSTheta + kAngTolerance)
           && (fSTheta > halfAngTolerance) )
    {
      // In tolerance of stheta
      // If entering through solid [r,phi] => 0 to in
      // else try 2nd root

      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
      if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
        || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
        || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
      {
        if (!fFullPhiSphere && rho2)  // Check phi intersection
        {
          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
          if (cosPsi >= cosHDPhiIT)
          {
            return 0 ;
          }
        }
        else
        {
          return 0 ;
        }
      }

      // Not entering immediately/travelling through

      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
      if (t1)
      { 
        b  = t2/t1 ;
        c  = dist2STheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d = std::sqrt(d2) ;
          s = -b + d ;
          if ( (s >= halfCarTolerance) && (s < snxt) && (fSTheta < halfpi) )
          {  // ^^^^^^^^^^^^^^^^^^^^^  shouldn't it be >=0 instead ?
            xi    = p.x() + s*v.x() ;
            yi    = p.y() + s*v.y() ;
            zi    = p.z() + s*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(fSTheta - halfpi) <= 0) )
            {
              if ( !fFullPhiSphere && rhoi2 )    // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if ( cosPsi >= cosHDPhiOT )
                {
                  snxt = s ;
                }
              }
              else
              {
                snxt = s ;
              }
            }
          }
        }
      }
    }   
    else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
    {

      // In tolerance of etheta
      // If entering through solid [r,phi] => 0 to in
      // else try 2nd root

      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;

      if (   ((t2<0) && (eTheta < halfpi)
          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
        ||   ((t2>=0) && (eTheta > halfpi)
          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
        ||   ((v.z()>0) && (eTheta == halfpi)
          && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))  )
      {
        if (!fFullPhiSphere && rho2)   // Check phi intersection
        {
          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
          if (cosPsi >= cosHDPhiIT)
          {
            return 0 ;
          }
        }
        else
        {
          return 0 ;
        }
      }

      // Not entering immediately/travelling through

      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
      if (t1)
      { 
        b  = t2/t1 ;
        c  = dist2ETheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d = std::sqrt(d2) ;
          s = -b + d ;
        
          if ( (s >= halfCarTolerance)
            && (s < snxt) && (eTheta > halfpi) )
          {
            xi    = p.x() + s*v.x() ;
            yi    = p.y() + s*v.y() ;
            zi    = p.z() + s*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(eTheta - halfpi) <= 0) )
            {
              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if (cosPsi >= cosHDPhiOT)
                {
                  snxt = s ;
                }
              }
              else
              {
                snxt = s ;
              }
            }
          }
        } 
      }   
    }  
    else
    {
      // stheta+tol<theta<etheta-tol
      // For BOTH stheta & etheta check 2nd root for validity [r,phi]

      t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
      t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
      if (t1)
      { 
        b  = t2/t1;
        c  = dist2STheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d = std::sqrt(d2) ;
          s = -b + d ;    // second root

          if ((s >= 0) && (s < snxt))
          {
            xi    = p.x() + s*v.x() ;
            yi    = p.y() + s*v.y() ;
            zi    = p.z() + s*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(fSTheta - halfpi) <= 0) )
            {
              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if (cosPsi >= cosHDPhiOT)
                {
                  snxt = s ;
                }
              }
              else
              {
                snxt = s ;
              }
            }
          }
        }
      }        
      t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
      t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
      if (t1)
      {   
        b  = t2/t1 ;
        c  = dist2ETheta/t1 ;
        d2 = b*b - c ;

        if (d2 >= 0)
        {
          d = std::sqrt(d2) ;
          s = -b + d;    // second root

          if ((s >= 0) && (s < snxt))
          {
            xi    = p.x() + s*v.x() ;
            yi    = p.y() + s*v.y() ;
            zi    = p.z() + s*v.z() ;
            rhoi2 = xi*xi + yi*yi   ;
            radi2 = rhoi2 + zi*zi   ;

            if ( (radi2 <= tolORMax2)
              && (radi2 >= tolORMin2)
              && (zi*(eTheta - halfpi) <= 0) )
            {
              if (!fFullPhiSphere && rhoi2)   // Check phi intersection
              {
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
                if ( cosPsi >= cosHDPhiOT )
                {
                  snxt=s;
                }
              }
              else
              {
                snxt = s ;
              }
            }
          }
        }
      }
    }  
  }
  return snxt;
}

//////////////////////////////////////////////////////////////////////
//
// Calculate distance (<= actual) to closest surface of shape from outside
// - Calculate distance to radial planes
// - Only to phi planes if outside phi extent
// - Only to theta planes if outside theta extent
// - Return 0 if point inside

G4double G4Sphere::DistanceToIn( const G4ThreeVector& p ) const
{
  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
  G4double rho2,rds,rho;
  G4double cosPsi;
  G4double pTheta,dTheta1,dTheta2;
  rho2=p.x()*p.x()+p.y()*p.y();
  rds=std::sqrt(rho2+p.z()*p.z());
  rho=std::sqrt(rho2);

  //
  // Distance to r shells
  //    
  if (fRmin)
  {
    safeRMin=fRmin-rds;
    safeRMax=rds-fRmax;
    if (safeRMin>safeRMax)
    {
      safe=safeRMin;
    }
    else
    {
      safe=safeRMax;
    }
  }
  else
  {
    safe=rds-fRmax;
  }

  //
  // Distance to phi extent
  //
  if (!fFullPhiSphere && rho)
  {
    // Psi=angle from central phi to point
    //
    cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
    if (cosPsi<std::cos(hDPhi))
    {
      // Point lies outside phi range
      //
      if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
      {
        safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
      }
      else
      {
        safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
      }
      if (safePhi>safe)  { safe=safePhi; }
    }
  }
  //
  // Distance to Theta extent
  //    
  if ((rds!=0.0) && (!fFullThetaSphere))
  {
    pTheta=std::acos(p.z()/rds);
    if (pTheta<0)  { pTheta+=pi; }
    dTheta1=fSTheta-pTheta;
    dTheta2=pTheta-eTheta;
    if (dTheta1>dTheta2)
    {
      if (dTheta1>=0)             // WHY ???????????
      {
        safeTheta=rds*std::sin(dTheta1);
        if (safe<=safeTheta)
        {
          safe=safeTheta;
        }
      }
    }
    else
    {
      if (dTheta2>=0)
      {
        safeTheta=rds*std::sin(dTheta2);
        if (safe<=safeTheta)
        {
          safe=safeTheta;
        }
      }
    }
  }

  if (safe<0)  { safe=0; }
  return safe;
}

/////////////////////////////////////////////////////////////////////
//
// Calculate distance to surface of shape from 'inside', allowing for tolerance
// - Only Calc rmax intersection if no valid rmin intersection

G4double G4Sphere::DistanceToOut( const G4ThreeVector& p,
                                  const G4ThreeVector& v,
                                  const G4bool calcNorm,
                                        G4bool *validNorm,
                                        G4ThreeVector *n   ) const
{
  G4double snxt = kInfinity;     // snxt is default return value
  G4double sphi= kInfinity,stheta= kInfinity;
  ESide side=kNull,sidephi=kNull,sidetheta=kNull;  

  static const G4double halfCarTolerance = kCarTolerance*0.5;
  static const G4double halfAngTolerance = kAngTolerance*0.5;
  const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
  const G4double halfRminTolerance = fRminTolerance*0.5;
  const G4double Rmax_plus  = fRmax + halfRmaxTolerance;
  const G4double Rmin_minus = (fRmin) ? fRmin-halfRminTolerance : 0;
  G4double t1,t2;
  G4double b,c,d;

  // Variables for phi intersection:

  G4double pDistS,compS,pDistE,compE,sphi2,vphi;
    
  G4double rho2,rad2,pDotV2d,pDotV3d,pTheta;

  G4double tolSTheta=0.,tolETheta=0.;
  G4double xi,yi,zi;      // Intersection point

  // Theta precals
  //
  G4double rhoSecTheta;
  G4double dist2STheta, dist2ETheta, distTheta;
  G4double d2,s;

  // General Precalcs
  //
  rho2 = p.x()*p.x()+p.y()*p.y();
  rad2 = rho2+p.z()*p.z();

  pTheta = std::atan2(std::sqrt(rho2),p.z());

  pDotV2d = p.x()*v.x()+p.y()*v.y();
  pDotV3d = pDotV2d+p.z()*v.z();

  // Theta precalcs
    
  if ( !fFullThetaSphere )
  {
    tolSTheta = fSTheta - halfAngTolerance;
    tolETheta = eTheta + halfAngTolerance;
  }
    
  // Radial Intersections from G4Sphere::DistanceToIn
  //
  // Outer spherical shell intersection
  // - Only if outside tolerant fRmax
  // - Check for if inside and outer G4Sphere heading through solid (-> 0)
  // - No intersect -> no intersection with G4Sphere
  //
  // Shell eqn: x^2+y^2+z^2=RSPH^2
  //
  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
  //
  // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
  // =>      rad2        +2s(pDotV3d)       +s^2                =R^2
  //
  // => s=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2))

  if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
  {
    c = rad2 - fRmax*fRmax;

    if (c < fRmaxTolerance*fRmax) 
    {
      // Within tolerant Outer radius 
      // 
      // The test is
      //     rad  - fRmax < 0.5*kRadTolerance
      // =>  rad  < fRmax + 0.5*kRadTol
      // =>  rad2 < (fRmax + 0.5*kRadTol)^2
      // =>  rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol
      // =>  rad2 - fRmax^2    <~    fRmax*kRadTol 

      d2 = pDotV3d*pDotV3d - c;

      if( (c >- fRmaxTolerance*fRmax)       // on tolerant surface
       && ((pDotV3d >=0) || (d2 < 0)) )     // leaving outside from Rmax 
                                            // not re-entering
      {
        if(calcNorm)
        {
          *validNorm = true ;
          *n         = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
        }
        return snxt = 0;
      }
      else 
      {
        snxt = -pDotV3d+std::sqrt(d2);    // second root since inside Rmax
        side =  kRMax ; 
      }
    }

    // Inner spherical shell intersection:
    // Always first >=0 root, because would have passed
    // from outside of Rmin surface .

    if (fRmin)
    {
      c  = rad2 - fRmin*fRmin;
      d2 = pDotV3d*pDotV3d - c;

      if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin
      {
        if ( (c < fRminTolerance*fRmin)              // leaving from Rmin
          && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
        {
          if(calcNorm)  { *validNorm = false; }  // Rmin surface is concave         
          return snxt = 0 ;
        }
        else
        {  
          if ( d2 >= 0. )
          {
            s = -pDotV3d-std::sqrt(d2);

            if ( s >= 0. )     // Always intersect Rmin first
            {
              snxt = s ;
              side = kRMin ;
            }
          }
        }
      }
    }
  }

  // Theta segment intersection

  if ( !fFullThetaSphere )
  {
    // Intersection with theta surfaces
    //
    // Known failure cases:
    // o  Inside tolerance of stheta surface, skim
    //    ~parallel to cone and Hit & enter etheta surface [& visa versa]
    //
    //    To solve: Check 2nd root of etheta surface in addition to stheta
    //
    // o  start/end theta is exactly pi/2 
    //
    // Intersections with cones
    //
    // Cone equation: x^2+y^2=z^2tan^2(t)
    //
    // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t)
    //
    // => (px^2+py^2-pz^2tan^2(t))+2s(pxvx+pyvy-pzvztan^2(t))
    //       + s^2(vx^2+vy^2-vz^2tan^2(t)) = 0
    //
    // => s^2(1-vz^2(1+tan^2(t))+2s(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0
    //
  
    if(fSTheta) // intersection with first cons
    {
      if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
      {
        if( v.z() > 0. ) 
        {
          if ( std::fabs( p.z() ) <= halfRmaxTolerance )
          {
            if(calcNorm)
            {
              *validNorm = true;
              *n = G4ThreeVector(0.,0.,1.);
            }
            return snxt = 0 ;
          }  
          stheta    = -p.z()/v.z();
          sidetheta = kSTheta;
        }
      }
      else // kons is not plane 
      {
        t1          = 1-v.z()*v.z()*(1+tanSTheta2);
        t2          = pDotV2d-p.z()*v.z()*tanSTheta2;  // ~vDotN if p on cons
        dist2STheta = rho2-p.z()*p.z()*tanSTheta2;     // t3

        distTheta = std::sqrt(rho2)-p.z()*tanSTheta;

        if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
        {                                      // v parallel to kons
          if( v.z() > 0. )
          {
            if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
            {
              if( (fSTheta < halfpi) && (p.z() > 0.) )
              {
                if( calcNorm )  { *validNorm = false; }
                return snxt = 0.;
              }
              else if( (fSTheta > halfpi) && (p.z() <= 0) )
              {
                if( calcNorm ) 
                {
                  *validNorm = true;
                  if (rho2)
                  {
                    rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
                   
                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                        p.y()/rhoSecTheta,
                                        std::sin(fSTheta)  );
                  }
                  else *n = G4ThreeVector(0.,0.,1.);
                }
                return snxt = 0.;               
              }
            }
            stheta    = -0.5*dist2STheta/t2;
            sidetheta = kSTheta;
          }  
        }      // 2nd order equation, 1st root of fSTheta cone,
        else   // 2nd if 1st root -ve
        {
          if( std::fabs(distTheta) < halfRmaxTolerance )
          {
            if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave
            {
              if( calcNorm ) 
              {
                *validNorm = true;
                if (rho2)
                {
                  rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
                   
                  *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                      p.y()/rhoSecTheta,
                                      std::sin(fSTheta)  );
                }
                else  { *n = G4ThreeVector(0.,0.,1.); }
              }
              return snxt = 0.;
            }
            else if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) ) // leave
            {
              if( calcNorm )  { *validNorm = false; }
              return snxt = 0.;
            }                               
          }
          b  = t2/t1;
          c  = dist2STheta/t1;
          d2 = b*b - c ;

          if ( d2 >= 0. )
          {
            d = std::sqrt(d2);

            if( fSTheta > halfpi )
            {
              s = -b - d;         // First root

              if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
               ||  (s < 0.)  || ( (s > 0.) && (p.z() + s*v.z() > 0.) )     ) 
              {
                s = -b + d ; // 2nd root
              }
              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() <= 0.) )  
              {
                stheta    = s;
                sidetheta = kSTheta;
              }
            }
            else // sTheta < pi/2, concave surface, no normal
            {
              s = -b - d;         // First root

              if ( ( (std::fabs(s) < halfRmaxTolerance) && (t2 >= 0.) )
                || (s < 0.) || ( (s > 0.) && (p.z() + s*v.z() < 0.) )   )
              {
                s = -b + d ; // 2nd root
              }
              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() >= 0.) )  
              {
                stheta    = s;
                sidetheta = kSTheta;
              }            
            }
          }
        }
      }
    }
    if (eTheta < pi) // intersection with second cons
    {
      if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
      {
        if( v.z() < 0. ) 
        {
          if ( std::fabs( p.z() ) <= halfRmaxTolerance )
          {
            if(calcNorm)
            {
              *validNorm = true;
              *n = G4ThreeVector(0.,0.,-1.);
            }
            return snxt = 0 ;
          }  
          s = -p.z()/v.z();

          if( s < stheta )
          {
            stheta    = s;
            sidetheta = kETheta;
          }
        }
      }
      else // kons is not plane 
      {
        t1          = 1-v.z()*v.z()*(1+tanETheta2);
        t2          = pDotV2d-p.z()*v.z()*tanETheta2;  // ~vDotN if p on cons
        dist2ETheta = rho2-p.z()*p.z()*tanETheta2;     // t3

        distTheta = std::sqrt(rho2)-p.z()*tanETheta;

        if( std::fabs(t1) < halfAngTolerance ) // 1st order equation,
        {                                      // v parallel to kons
          if( v.z() < 0. )
          {
            if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface
            {
              if( (eTheta > halfpi) && (p.z() < 0.) )
              {
                if( calcNorm )  { *validNorm = false; }
                return snxt = 0.;
              }
              else if ( (eTheta < halfpi) && (p.z() >= 0) )
              {
                if( calcNorm ) 
                {
                  *validNorm = true;
                  if (rho2)
                  {
                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                        p.y()/rhoSecTheta,
                                        -sinETheta  );
                  }
                  else  { *n = G4ThreeVector(0.,0.,-1.); }
                }
                return snxt = 0.;               
              }
            }
            s = -0.5*dist2ETheta/t2;

            if( s < stheta )
            {
              stheta    = s;
              sidetheta = kETheta;
            }
          }  
        }      // 2nd order equation, 1st root of fSTheta cone
        else   // 2nd if 1st root -ve
        {
          if ( std::fabs(distTheta) < halfRmaxTolerance )
          {
            if( (eTheta < halfpi) && (t2 >= 0.) ) // leave
            {
              if( calcNorm ) 
              {
                *validNorm = true;
                if (rho2)
                {
                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                        p.y()/rhoSecTheta,
                                        -sinETheta  );
                }
                else *n = G4ThreeVector(0.,0.,-1.);
              }                           
              return snxt = 0.;
            }
            else if ( (eTheta > halfpi)
                   && (t2 < 0.) && (p.z() <=0.) ) // leave
            {
              if( calcNorm )  { *validNorm = false; }
              return snxt = 0.;
            }                               
          }
          b  = t2/t1;
          c  = dist2ETheta/t1;
          d2 = b*b - c ;

          if ( d2 >= 0. )
          {
            d = std::sqrt(d2);

            if( eTheta < halfpi )
            {
              s = -b - d;         // First root

              if( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
               || (s < 0.) ) 
              {
                s = -b + d ; // 2nd root
              }
              if( s > halfRmaxTolerance )  
              {
                if( s < stheta )
                {
                  stheta    = s;
                  sidetheta = kETheta;
                }
              }
            }
            else // sTheta+fDTheta > pi/2, concave surface, no normal
            {
              s = -b - d;         // First root

              if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 >= 0.))
                || (s < 0.) || ( (s > 0.) && (p.z() + s*v.z() > 0.) ) )
              {
                s = -b + d ; // 2nd root
              }
              if( (s > halfRmaxTolerance) && (p.z() + s*v.z() <= 0.) )  
              {
                if( s < stheta )
                {
                  stheta    = s;
                  sidetheta = kETheta;
                }
              }            
            }
          }
        }
      }
    }

  } // end theta intersections

  // Phi Intersection
    
  if ( !fFullPhiSphere )
  {
    if ( p.x() || p.y() ) // Check if on z axis (rho not needed later)
    {
      // pDist -ve when inside

      pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
      pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;

      // Comp -ve when in direction of outwards normal

      compS   = -sinSPhi*v.x()+cosSPhi*v.y() ;
      compE   =  sinEPhi*v.x()-cosEPhi*v.y() ;
      sidephi = kNull ;

      if ( (pDistS <= 0) && (pDistE <= 0) )
      {
        // Inside both phi *full* planes

        if ( compS < 0 )
        {
          sphi = pDistS/compS ;
          xi   = p.x()+sphi*v.x() ;
          yi   = p.y()+sphi*v.y() ;

          // Check intersection with correct half-plane (if not -> no intersect)
          //
          if( (std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance) )
          {
            vphi = std::atan2(v.y(),v.x());
            sidephi = kSPhi;
            if ( ( (fSPhi-halfAngTolerance) <= vphi)
              && ( (ePhi+halfAngTolerance)  >= vphi) )
            {
              sphi = kInfinity;
            }
          }
          else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
          {
            sphi=kInfinity;
          }
          else
          {
            sidephi = kSPhi ;
            if ( pDistS > -halfCarTolerance)  { sphi = 0; } // Leave by sphi 
          }
        }
        else  { sphi = kInfinity; }

        if ( compE < 0 )
        {
          sphi2=pDistE/compE ;
          if (sphi2 < sphi) // Only check further if < starting phi intersection
          {
            xi = p.x()+sphi2*v.x() ;
            yi = p.y()+sphi2*v.y() ;

            // Check intersection with correct half-plane
            //
            if ((std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance))
            {
              // Leaving via ending phi
              //
              vphi = std::atan2(v.y(),v.x()) ;
               
              if( !((fSPhi-halfAngTolerance <= vphi)
                  &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
              { 
                sidephi = kEPhi;
                if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
                else                                { sphi = 0.0;   }
              }
            }
            else if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
            {
              sidephi = kEPhi ;
              if ( pDistE <= -halfCarTolerance )
              {
                sphi=sphi2;
              }
              else 
              {
                sphi = 0 ;
              }
            }
          }
        }        
      }
      else if ((pDistS >= 0) && (pDistE >= 0)) // Outside both *full* phi planes
      {
        if ( pDistS <= pDistE )
        {
          sidephi = kSPhi ;
        }
        else
        {
          sidephi = kEPhi ;
        }
        if ( fDPhi > pi )
        {
          if ( (compS < 0) && (compE < 0) )  { sphi = 0; }
          else                               { sphi = kInfinity; }
        }
        else
        {
          // if towards both >=0 then once inside (after error)
          // will remain inside

          if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
          else                                { sphi = 0; }
        }    
      }
      else if ( (pDistS > 0) && (pDistE < 0) )
      {
        // Outside full starting plane, inside full ending plane

        if ( fDPhi > pi )
        {
          if ( compE < 0 )
          {
            sphi = pDistE/compE ;
            xi   = p.x() + sphi*v.x() ;
            yi   = p.y() + sphi*v.y() ;

            // Check intersection in correct half-plane
            // (if not -> not leaving phi extent)
            //
            if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
            {
              vphi = std::atan2(v.y(),v.x());
              sidephi = kSPhi;
              if ( ( (fSPhi-halfAngTolerance) <= vphi)
                && ( (ePhi+halfAngTolerance)  >= vphi) )
              {
                sphi = kInfinity;
              }
            }
            else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
            {
              sphi = kInfinity ;
            }
            else // Leaving via Ending phi
            {
              sidephi = kEPhi ;
              if ( pDistE > -halfCarTolerance )  { sphi = 0.; }
            }
          }
          else
          {
            sphi = kInfinity ;
          }
        }
        else
        {
          if ( compS >= 0 )
          {
            if ( compE < 0 )
            {            
              sphi = pDistE/compE ;
              xi   = p.x() + sphi*v.x() ;
              yi   = p.y() + sphi*v.y() ;

              // Check intersection in correct half-plane
              // (if not -> remain in extent)
              //
              if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
              {
                vphi = std::atan2(v.y(),v.x());
                sidephi = kSPhi;
                if ( ( (fSPhi-halfAngTolerance) <= vphi)
                  && ( (ePhi+halfAngTolerance)  >= vphi) )
                {
                  sphi = kInfinity;
                }
              }
              else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
              {
                sphi=kInfinity;
              }
              else // otherwise leaving via Ending phi
              {
                sidephi = kEPhi ;
              }
            }
            else sphi=kInfinity;
          }
          else // leaving immediately by starting phi
          {
            sidephi = kSPhi ;
            sphi    = 0 ;
          }
        }
      }
      else
      {
        // Must be pDistS < 0 && pDistE > 0
        // Inside full starting plane, outside full ending plane

        if ( fDPhi > pi )
        {
          if ( compS < 0 )
          {
            sphi=pDistS/compS;
            xi=p.x()+sphi*v.x();
            yi=p.y()+sphi*v.y();
            
            // Check intersection in correct half-plane
            // (if not -> not leaving phi extent)
            //
            if( (std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance) )
            {
              vphi = std::atan2(v.y(),v.x()) ;
              sidephi = kSPhi;
              if ( ( (fSPhi-halfAngTolerance) <= vphi)
                && ( (ePhi+halfAngTolerance)  >= vphi) )
              {
              sphi = kInfinity;
              }
            }
            else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
            {
              sphi = kInfinity ;
            }
            else  // Leaving via Starting phi
            {
              sidephi = kSPhi ;
              if ( pDistS > -halfCarTolerance )  { sphi = 0; }
            }
          }
          else
          {
            sphi = kInfinity ;
          }
        }
        else
        {
          if ( compE >= 0 )
          {
            if ( compS < 0 )
            {
              sphi = pDistS/compS ;
              xi   = p.x()+sphi*v.x() ;
              yi   = p.y()+sphi*v.y() ;
              
              // Check intersection in correct half-plane
              // (if not -> remain in extent)
              //
              if((std::abs(xi)<=kCarTolerance) && (std::abs(yi)<=kCarTolerance))
              {
                vphi = std::atan2(v.y(),v.x()) ;
                sidephi = kSPhi;
                if ( ( (fSPhi-halfAngTolerance) <= vphi)
                  && ( (ePhi+halfAngTolerance)  >= vphi) )
                {
                  sphi = kInfinity;
                }
              }
              else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
              {
                sphi = kInfinity ;
              }
              else // otherwise leaving via Starting phi
              {
                sidephi = kSPhi ;
              }
            }
            else
            {
              sphi = kInfinity ;
            }
          }
          else // leaving immediately by ending
          {
            sidephi = kEPhi ;
            sphi    = 0     ;
          }
        }
      }      
    }
    else
    {
      // On z axis + travel not || to z axis -> if phi of vector direction
      // within phi of shape, Step limited by rmax, else Step =0

      if ( v.x() || v.y() )
      {
        vphi = std::atan2(v.y(),v.x()) ;
        if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
        {
          sphi = kInfinity;
        }
        else
        {
          sidephi = kSPhi ; // arbitrary 
          sphi    = 0     ;
        }
      }
      else  // travel along z - no phi intersection
      {
        sphi = kInfinity ;
      }
    }
    if ( sphi < snxt )  // Order intersecttions
    {
      snxt = sphi ;
      side = sidephi ;
    }
  }
  if (stheta < snxt ) // Order intersections
  {
    snxt = stheta ;
    side = sidetheta ;
  }

  if (calcNorm)    // Output switch operator
  {
    switch( side )
    {
      case kRMax:
        xi=p.x()+snxt*v.x();
        yi=p.y()+snxt*v.y();
        zi=p.z()+snxt*v.z();
        *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
        *validNorm=true;
        break;

      case kRMin:
        *validNorm=false;  // Rmin is concave
        break;

      case kSPhi:
        if ( fDPhi <= pi )     // Normal to Phi-
        {
          *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
          *validNorm=true;
        }
        else  { *validNorm=false; }
        break ;

      case kEPhi:
        if ( fDPhi <= pi )      // Normal to Phi+
        {
          *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
          *validNorm=true;
        }
        else  { *validNorm=false; }
        break;

      case kSTheta:
        if( fSTheta == halfpi )
        {
          *n=G4ThreeVector(0.,0.,1.);
          *validNorm=true;
        }
        else if ( fSTheta > halfpi )
        {
          xi = p.x() + snxt*v.x();
          yi = p.y() + snxt*v.y();
          rho2=xi*xi+yi*yi;
          if (rho2)
          { 
            rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
            *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
                               -tanSTheta/std::sqrt(1+tanSTheta2));
          }
          else
          {
            *n = G4ThreeVector(0.,0.,1.);
          }
          *validNorm=true;
        }
        else  { *validNorm=false; }  // Concave STheta cone
        break;

      case kETheta:
        if( eTheta == halfpi )
        {
          *n         = G4ThreeVector(0.,0.,-1.);
          *validNorm = true;
        }
        else if ( eTheta < halfpi )
        {
          xi=p.x()+snxt*v.x();
          yi=p.y()+snxt*v.y();
          rho2=xi*xi+yi*yi;
          if (rho2)
          { 
            rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
            *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
                               -tanETheta/std::sqrt(1+tanETheta2) );
          }
          else
          {
            *n = G4ThreeVector(0.,0.,-1.);
          }
          *validNorm=true;
        }
        else  { *validNorm=false; }   // Concave ETheta cone
        break;

      default:
        G4cout.precision(16);
        G4cout << G4endl;
        DumpInfo();
        G4cout << "Position:"  << G4endl << G4endl;
        G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
        G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
        G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
        G4cout << "Direction:" << G4endl << G4endl;
        G4cout << "v.x() = "   << v.x() << G4endl;
        G4cout << "v.y() = "   << v.y() << G4endl;
        G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
        G4cout << "Proposed distance :" << G4endl << G4endl;
        G4cout << "snxt = "    << snxt/mm << " mm" << G4endl << G4endl;
        G4Exception("G4Sphere::DistanceToOut(p,v,..)",
                    "Notification", JustWarning,
                    "Undefined side for valid surface normal to solid.");
        break;
    }
  }
  if (snxt == kInfinity)
  {
    G4cout.precision(24);
    G4cout << G4endl;
    DumpInfo();
    G4cout << "Position:"  << G4endl << G4endl;
    G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl;
    G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl;
    G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl;
    G4cout << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm << " mm" 
           << G4endl << G4endl;
    G4cout << "Direction:" << G4endl << G4endl;
    G4cout << "v.x() = "   << v.x() << G4endl;
    G4cout << "v.y() = "   << v.y() << G4endl;
    G4cout << "v.z() = "   << v.z() << G4endl << G4endl;
    G4cout << "Proposed distance :" << G4endl << G4endl;
    G4cout << "snxt = "    << snxt/mm << " mm" << G4endl << G4endl;
    G4Exception("G4Sphere::DistanceToOut(p,v,..)",
                "Notification", JustWarning,
                "Logic error: snxt = kInfinity  ???");
  }

  return snxt;
}

/////////////////////////////////////////////////////////////////////////
//
// Calculate distance (<=actual) to closest surface of shape from inside

G4double G4Sphere::DistanceToOut( const G4ThreeVector& p ) const
{
  G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
  G4double rho2,rds,rho;
  G4double pTheta,dTheta1,dTheta2;
  rho2=p.x()*p.x()+p.y()*p.y();
  rds=std::sqrt(rho2+p.z()*p.z());
  rho=std::sqrt(rho2);

#ifdef G4CSGDEBUG
  if( Inside(p) == kOutside )
  {
     G4cout.precision(16) ;
     G4cout << G4endl ;
     DumpInfo();
     G4cout << "Position:"  << G4endl << G4endl ;
     G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
     G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
     G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
     G4Exception("G4Sphere::DistanceToOut(p)",
                 "Notification", JustWarning, "Point p is outside !?" );
  }
#endif

  //
  // Distance to r shells
  //    
  if (fRmin)
  {
    safeRMin=rds-fRmin;
    safeRMax=fRmax-rds;
    if (safeRMin<safeRMax)
    {
      safe=safeRMin;
    }
    else
    {
      safe=safeRMax;
    }
  }
  else
  {
    safe=fRmax-rds;
  }

  //
  // Distance to phi extent
  //
  if (!fFullPhiSphere && rho)
  {
    if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
    {
      safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
    }
    else
    {
      safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
    }
    if (safePhi<safe)  { safe=safePhi; }
  }

  //
  // Distance to Theta extent
  //    
  if (rds)
  {
    pTheta=std::acos(p.z()/rds);
    if (pTheta<0)  { pTheta+=pi; }
    dTheta1=pTheta-fSTheta;
    dTheta2=eTheta-pTheta;
    if (dTheta1<dTheta2)  { safeTheta=rds*std::sin(dTheta1); }
    else                  { safeTheta=rds*std::sin(dTheta2); }
    if (safe>safeTheta)   { safe=safeTheta; }
  }

  if (safe<0)  { safe=0; }
  return safe;
}

//////////////////////////////////////////////////////////////////////////
//
// Create a List containing the transformed vertices
// Ordering [0-3] -fDz cross section
//          [4-7] +fDz cross section such that [0] is below [4],
//                                             [1] below [5] etc.
// Note:
//  Caller has deletion resposibility
//  Potential improvement: For last slice, use actual ending angle
//                         to avoid rounding error problems.

G4ThreeVectorList*
G4Sphere::CreateRotatedVertices( const G4AffineTransform& pTransform,
                                       G4int& noPolygonVertices ) const
{
  G4ThreeVectorList *vertices;
  G4ThreeVector vertex;
  G4double meshAnglePhi,meshRMax,crossAnglePhi,
           coscrossAnglePhi,sincrossAnglePhi,sAnglePhi;
  G4double meshTheta,crossTheta,startTheta;
  G4double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ;
  G4int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections;

  // Phi cross sections
    
  noPhiCrossSections = G4int(fDPhi/kMeshAngleDefault)+1;
    
  if (noPhiCrossSections<kMinMeshSections)
  {
    noPhiCrossSections=kMinMeshSections;
  }
  else if (noPhiCrossSections>kMaxMeshSections)
  {
    noPhiCrossSections=kMaxMeshSections;
  }
  meshAnglePhi=fDPhi/(noPhiCrossSections-1);
    
  // If complete in phi, set start angle such that mesh will be at fRMax
  // on the x axis. Will give better extent calculations when not rotated.
    
  if (fFullPhiSphere)
  {
    sAnglePhi = -meshAnglePhi*0.5;
  }
    else
  {
    sAnglePhi=fSPhi;
  }    

  // Theta cross sections
    
  noThetaSections = G4int(fDTheta/kMeshAngleDefault)+1;
    
  if (noThetaSections<kMinMeshSections)
  {
    noThetaSections=kMinMeshSections;
  }
  else if (noThetaSections>kMaxMeshSections)
  {
    noThetaSections=kMaxMeshSections;
  }
  meshTheta=fDTheta/(noThetaSections-1);
    
  // If complete in Theta, set start angle such that mesh will be at fRMax
  // on the z axis. Will give better extent calculations when not rotated.
    
  if (fFullThetaSphere)
  {
    startTheta = -meshTheta*0.5;
  }
  else
  {
    startTheta=fSTheta;
  }    

  meshRMax = (meshAnglePhi >= meshTheta) ?
             fRmax/std::cos(meshAnglePhi*0.5) : fRmax/std::cos(meshTheta*0.5);
  G4double* cosCrossTheta = new G4double[noThetaSections];
  G4double* sinCrossTheta = new G4double[noThetaSections];    
  vertices=new G4ThreeVectorList();
  vertices->reserve(noPhiCrossSections*(noThetaSections*2));
  if (vertices && cosCrossTheta && sinCrossTheta)
  {
    for (crossSectionPhi=0;
         crossSectionPhi<noPhiCrossSections; crossSectionPhi++)
    {
      crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
      coscrossAnglePhi=std::cos(crossAnglePhi);
      sincrossAnglePhi=std::sin(crossAnglePhi);
      for (crossSectionTheta=0;
           crossSectionTheta<noThetaSections;crossSectionTheta++)
      {
        // Compute coordinates of cross section at section crossSectionPhi
        //
        crossTheta=startTheta+crossSectionTheta*meshTheta;
        cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
        sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);

        rMinX=fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
        rMinY=fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
        rMinZ=fRmin*cosCrossTheta[crossSectionTheta];
        
        vertex=G4ThreeVector(rMinX,rMinY,rMinZ);
        vertices->push_back(pTransform.TransformPoint(vertex));
        
      }    // Theta forward 
    
      for (crossSectionTheta=noThetaSections-1;
           crossSectionTheta>=0; crossSectionTheta--)
      {
        rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
        rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
        rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta];
        
        vertex=G4ThreeVector(rMaxX,rMaxY,rMaxZ);
        vertices->push_back(pTransform.TransformPoint(vertex));

      }   // Theta back 
    }   // Phi
    noPolygonVertices = noThetaSections*2 ;
  }
  else
  {
    DumpInfo();
    G4Exception("G4Sphere::CreateRotatedVertices()",
                "FatalError", FatalException,
                "Error in allocation of vertices. Out of memory !");
  }

  delete [] cosCrossTheta;
  delete [] sinCrossTheta;

  return vertices;
}

//////////////////////////////////////////////////////////////////////////
//
// G4EntityType

G4GeometryType G4Sphere::GetEntityType() const
{
  return G4String("G4Sphere");
}

//////////////////////////////////////////////////////////////////////////
//
// Stream object contents to an output stream

std::ostream& G4Sphere::StreamInfo( std::ostream& os ) const
{
  os << "-----------------------------------------------------------\n"
     << "    *** Dump for solid - " << GetName() << " ***\n"
     << "    ===================================================\n"
     << " Solid type: G4Sphere\n"
     << " Parameters: \n"
     << "    inner radius: " << fRmin/mm << " mm \n"
     << "    outer radius: " << fRmax/mm << " mm \n"
     << "    starting phi of segment  : " << fSPhi/degree << " degrees \n"
     << "    delta phi of segment     : " << fDPhi/degree << " degrees \n"
     << "    starting theta of segment: " << fSTheta/degree << " degrees \n"
     << "    delta theta of segment   : " << fDTheta/degree << " degrees \n"
     << "-----------------------------------------------------------\n";

  return os;
}

////////////////////////////////////////////////////////////////////////////////
//
// GetPointOnSurface

G4ThreeVector G4Sphere::GetPointOnSurface() const
{
  G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
  G4double height1, height2, slant1, slant2, costheta, sintheta,theta,rRand;

  height1 = (fRmax-fRmin)*cosSTheta;
  height2 = (fRmax-fRmin)*cosETheta;
  slant1  = std::sqrt(sqr((fRmax - fRmin)*sinSTheta) + height1*height1);
  slant2  = std::sqrt(sqr((fRmax - fRmin)*sinETheta) + height2*height2);
  rRand   = RandFlat::shoot(fRmin,fRmax);
  
  aOne = fRmax*fRmax*fDPhi*(cosSTheta-cosETheta);
  aTwo = fRmin*fRmin*fDPhi*(cosSTheta-cosETheta);
  aThr = fDPhi*((fRmax + fRmin)*sinSTheta)*slant1;
  aFou = fDPhi*((fRmax + fRmin)*sinETheta)*slant2;
  aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
  
  phi = RandFlat::shoot(fSPhi, ePhi); 
  cosphi = std::cos(phi); 
  sinphi = std::sin(phi);
  theta = RandFlat::shoot(fSTheta,eTheta);
  costheta = std::cos(theta);
  sintheta = std::sqrt(1.-sqr(costheta));

  if(fFullPhiSphere) { aFiv = 0; }
  if(fSTheta == 0)   { aThr=0; }
  if(eTheta == pi) { aFou = 0; }
  if(fSTheta == halfpi) { aThr = pi*(fRmax*fRmax-fRmin*fRmin); }
  if(eTheta == halfpi)  { aFou = pi*(fRmax*fRmax-fRmin*fRmin); }

  chose = RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
  if( (chose>=0.) && (chose<aOne) )
  {
    return G4ThreeVector(fRmax*sintheta*cosphi,
                         fRmax*sintheta*sinphi, fRmax*costheta);
  }
  else if( (chose>=aOne) && (chose<aOne+aTwo) )
  {
    return G4ThreeVector(fRmin*sintheta*cosphi,
                         fRmin*sintheta*sinphi, fRmin*costheta);
  }
  else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
  {
    if (fSTheta != halfpi)
    {
      zRand = RandFlat::shoot(fRmin*cosSTheta,fRmax*cosSTheta);
      return G4ThreeVector(tanSTheta*zRand*cosphi,
                           tanSTheta*zRand*sinphi,zRand);
    }
    else
    {
      return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
    }    
  }
  else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
  {
    if(eTheta != halfpi)
    {
      zRand = RandFlat::shoot(fRmin*cosETheta, fRmax*cosETheta);
      return G4ThreeVector  (tanETheta*zRand*cosphi,
                             tanETheta*zRand*sinphi,zRand);
    }
    else
    {
      return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
    }
  }
  else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
  {
    return G4ThreeVector(rRand*sintheta*cosSPhi,
                         rRand*sintheta*sinSPhi,rRand*costheta);
  }
  else
  {
    return G4ThreeVector(rRand*sintheta*cosEPhi,
                         rRand*sintheta*sinEPhi,rRand*costheta);
  }
}

////////////////////////////////////////////////////////////////////////////////
//
// GetSurfaceArea

G4double G4Sphere::GetSurfaceArea()
{
  if(fSurfaceArea != 0.) {;}
  else
  {   
    G4double Rsq=fRmax*fRmax;
    G4double rsq=fRmin*fRmin;
         
    fSurfaceArea = fDPhi*(rsq+Rsq)*(cosSTheta - cosETheta);
    if(!fFullPhiSphere)
    {
      fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
    }
    if(fSTheta >0)
    {
      G4double acos1=std::acos( std::pow(sinSTheta,2) * std::cos(fDPhi)
                              + std::pow(cosSTheta,2));
      if(fDPhi>pi)
      { 
        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
      }
      else
      {
        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
      }
    }
    if(eTheta < pi)
    {
      G4double acos2=std::acos( std::pow(sinETheta,2) * std::cos(fDPhi)
                              + std::pow(cosETheta,2));
      if(fDPhi>pi)
      { 
        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
      }
      else
      {
        fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
      }
    }
  }
  return fSurfaceArea;
}

/////////////////////////////////////////////////////////////////////////////
//
// Methods for visualisation

G4VisExtent G4Sphere::GetExtent() const
{
  return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax );
}


void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
{
  scene.AddSolid (*this);
}

G4Polyhedron* G4Sphere::CreatePolyhedron () const
{
  return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta);
}

G4NURBS* G4Sphere::CreateNURBS () const
{
  return new G4NURBSbox (fRmax, fRmax, fRmax);       // Box for now!!!
}