//
// ********************************************************************
// * 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.68 2008/07/07 09:35:16 grichine Exp $
// GEANT4 tag $Name: geant4-09-02-ref-02 $
//
// class G4Sphere
//
// Implementation for G4Sphere class
//
// History:
//
// 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 <assert.h>

#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)
{
  fEpsilon = 1.0e-14;

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

  // Check radii

  if (pRmin<pRmax&&pRmin>=0)
  {
    fRmin=pRmin; fRmax=pRmax;
  }
  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

  if (pDPhi>=twopi)
  {
    fDPhi=twopi;
  }
  else if (pDPhi>0)
  {
    fDPhi=pDPhi;
  }
  else
  {
    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
           << "        Negative Z delta-Phi ! - "
           << pDPhi << G4endl;
    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
                "Invalid DPhi.");
  }

  // Convert fSPhi to 0-2PI

  if (pSPhi<0)
  {
    fSPhi=twopi-std::fmod(std::fabs(pSPhi),twopi);
  }
  else
  {
    fSPhi=std::fmod(pSPhi,twopi);
  }

  // Sphere is placed such that fSPhi+fDPhi>twopi !
  // fSPhi could be < 0 !!?
  //
  if (fSPhi+fDPhi>twopi) fSPhi-=twopi;

  // Check theta angles

  if (pSTheta<0 || pSTheta>pi)
  {
    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl;
    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
                "stheta outside 0-PI range.");
  }
  else
  {
    fSTheta=pSTheta;
  }

  if (pDTheta+pSTheta>=pi)
  {
    fDTheta=pi-pSTheta;
  }
  else if (pDTheta>0)
  {
    fDTheta=pDTheta;
  }
  else
  {
    G4cerr << "ERROR - G4Sphere()::G4Sphere(): " << GetName() << G4endl
           << "        Negative delta-Theta ! - "
           << pDTheta << G4endl;
    G4Exception("G4Sphere::G4Sphere()", "InvalidSetup", FatalException,
                "Invalid 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 ( fDPhi==twopi && fDTheta==pi)  // !pTransform.IsRotated() &&
  {
    // 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;

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

  //  if(rad2 >= 1.369e+19) DBG();
  //  G4double rad = std::sqrt(rad2);
  // Check radial surfaces
  // sets `in'

  if ( fRmin ) tolRMin = fRmin + kRadTolerance*0.5;
  else         tolRMin = 0 ;
  
  tolRMax = fRmax - kRadTolerance*0.5 ;
  //  const G4double  fractionTolerance = 1.0e-12;
  const G4double  flexRadMaxTolerance = // kRadTolerance;
    std::max(kRadTolerance, fEpsilon * fRmax);

  const G4double  Rmax_minus = fRmax - flexRadMaxTolerance*0.5;
  const G4double  flexRadMinTolerance = std::max(kRadTolerance, 
                     fEpsilon * fRmin);
  const G4double  Rmin_plus = (fRmin > 0) ? fRmin + flexRadMinTolerance*0.5 : 0 ;
    
if(rad2 <= Rmax_minus*Rmax_minus && rad2 >= Rmin_plus*Rmin_plus) in = kInside ; 

// if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin )  in = kInside ;
  // if ( rad <= tolRMax && rad >= tolRMin )  in = kInside ;
  else
  {
    tolRMax = fRmax + kRadTolerance*0.5 ;
    tolRMin = fRmin - kRadTolerance*0.5 ;

    if ( tolRMin < 0.0 ) tolRMin = 0.0 ;
    
     if ( rad2 <= tolRMax*tolRMax && rad2 >= tolRMin*tolRMin )  in = kSurface ;
    //  if ( rad <= tolRMax && rad >= tolRMin )  in = kSurface ;
    else                                                return in = kOutside ;
  }

  // Phi boundaries   : Do not check if it has no phi boundary!
  // (in != kOutside). It is new J.Apostolakis proposal of 30.10.03

  if ( ( fDPhi < twopi - kAngTolerance ) &&
       ( (p.x() != 0.0 ) || (p.y() != 0.0) ) )
  {
    pPhi = std::atan2(p.y(),p.x()) ;

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

  // Theta bondaries
  // (in!=kOutside)
  
  if ( (rho2 || p.z()) && fDTheta < pi - kAngTolerance*0.5 )
  {
    rho    = std::sqrt(rho2);
    pTheta = std::atan2(rho,p.z());

    if ( in == kInside )
    {
      if ( (pTheta < fSTheta + kAngTolerance*0.5)
        || (pTheta > fSTheta + fDTheta - kAngTolerance*0.5) )
      {
        if ( (pTheta >= fSTheta - kAngTolerance*0.5)
          && (pTheta <= fSTheta + fDTheta + kAngTolerance*0.5) )
        {
          in = kSurface ;
        }
        else
        {
          in = kOutside ;
        }
      }
    }
    else
    {
      if ( (pTheta < fSTheta - kAngTolerance*0.5)
        || (pTheta > fSTheta + fDTheta + kAngTolerance*0.5) )
      {
        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;
  G4double delta = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
  G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
  G4ThreeVector norm, sumnorm(0.,0.,0.);

  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 && (fDPhi < twopi || fDTheta < pi) )
  {
    pPhi = std::atan2(p.y(),p.x());

    if(pPhi  < fSPhi-dAngle)           pPhi     += twopi;
    else if(pPhi > fSPhi+fDPhi+dAngle) pPhi     -= twopi;
  }
  if ( fDPhi < twopi ) // && rho ) // old limitation against (0,0,z)
  {
    if ( rho )
    {
      distSPhi = std::fabs( pPhi - fSPhi ); 
      distEPhi = std::fabs(pPhi-fSPhi-fDPhi); 
    }
    else if( !fRmin )
    {
      distSPhi = 0.; 
      distEPhi = 0.; 
    }
    nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
  }        
  if ( fDTheta < pi ) // && rad ) // old limitation against (0,0,0)
  {
    if ( rho )
    {
      pTheta     = std::atan2(rho,p.z());
      distSTheta = std::fabs(pTheta-fSTheta); 
      distETheta = std::fabs(pTheta-fSTheta-fDTheta);
 
      nTs = G4ThreeVector(-std::cos(fSTheta)*p.x()/rho, //  *std::cos(pPhi),
                          -std::cos(fSTheta)*p.y()/rho, //  *std::sin(pPhi),
                           std::sin(fSTheta)                   );

      nTe = G4ThreeVector( std::cos(fSTheta+fDTheta)*p.x()/rho, // *std::cos(pPhi),
                           std::cos(fSTheta+fDTheta)*p.y()/rho, // *std::sin(pPhi),
                          -std::sin(fSTheta+fDTheta)                  );    
    }
    else if( !fRmin )
    {
      if ( fSTheta )  
      {              
        distSTheta = 0.;
	nTs = G4ThreeVector(0.,0.,-1.);
      }
      if ( fSTheta + fDTheta < pi ) // distETheta = 0.;
      {              
        distETheta = 0.;
	nTe = G4ThreeVector(0.,0.,1.);
      }
    }    
  }
  if( rad )  nR = G4ThreeVector(p.x()/rad,p.y()/rad,p.z()/rad);

  if( distRMax <= delta )
  {
    noSurfaces ++;
    sumnorm += nR;
  }
  if( fRmin && distRMin <= delta )
  {
    noSurfaces ++;
    sumnorm -= nR;
  }
  if( fDPhi < twopi )   
  {
    if (distSPhi <= dAngle)
    {
      noSurfaces ++;
      sumnorm += nPs;
    }
    if (distEPhi <= dAngle) 
    {
      noSurfaces ++;
      sumnorm += nPe;
    }
  }
  if ( fDTheta < pi )
  {
    if (distSTheta <= dAngle && fSTheta > 0.)
    {
      noSurfaces ++;
      if( rad <= delta && fDPhi >= twopi) sumnorm += nZ;
      else                                sumnorm += nTs;
    }
    if (distETheta <= dAngle && fSTheta+fDTheta < pi) 
    {
      noSurfaces ++;
      if( rad <= delta && fDPhi >= twopi) 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 (fDPhi<twopi&&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 (fDTheta<pi&&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(std::sin(fSPhi),-std::cos(fSPhi),0);
      break;
    case kNEPhi:
      norm=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
      break;
    case kNSTheta:
      norm=G4ThreeVector(-std::cos(fSTheta)*std::cos(pPhi),
                         -std::cos(fSTheta)*std::sin(pPhi),
                          std::sin(fSTheta)            );
      //  G4cout<<G4endl<<" case kNSTheta:"<<G4endl;
      //  G4cout<<"pPhi = "<<pPhi<<G4endl;
      //  G4cout<<"rad  = "<<rad<<G4endl;
      //  G4cout<<"pho  = "<<rho<<G4endl;
      //  G4cout<<"p:    "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl;
      //  G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl;
      break;
    case kNETheta:
      norm=G4ThreeVector( std::cos(fSTheta+fDTheta)*std::cos(pPhi),
                          std::cos(fSTheta+fDTheta)*std::sin(pPhi),
                         -std::sin(fSTheta+fDTheta)              );

      //  G4cout<<G4endl<<" case kNETheta:"<<G4endl;
      //  G4cout<<"pPhi = "<<pPhi<<G4endl;
      //  G4cout<<"rad  = "<<rad<<G4endl;
      //  G4cout<<"pho  = "<<rho<<G4endl;
      //  G4cout<<"p:    "<<p.x()<<"; "<<p.y()<<"; "<<p.z()<<G4endl;
      //  G4cout<<"norm: "<<norm.x()<<"; "<<norm.y()<<"; "<<norm.z()<<G4endl;
      break;
    default:
      DumpInfo();
      G4Exception("G4Sphere::ApproxSurfaceNormal()", "Notification", JustWarning,
                  "Undefined side for valid surface normal to solid.");
      break;    
  } // end case

  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 tolIRMin2, tolORMin2, tolORMax2, tolIRMax2 ;
  G4double tolSTheta=0., tolETheta=0. ;

  // Intersection point

  G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;

  // Phi intersection

  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi , Comp ; 

  // Phi flag and precalcs

  G4bool segPhi ;        
  G4double hDPhi, hDPhiOT, hDPhiIT, cPhi, sinCPhi=0., cosCPhi=0. ; 
  G4double cosHDPhiOT=0., cosHDPhiIT=0. ;
  G4double Dist, cosPsi ;

  G4bool segTheta ;                             // Theta flag and precals
  G4double tanSTheta, tanETheta ;
  G4double tanSTheta2, tanETheta2 ;
  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() ;

  // Radial Precalcs

  if (fRmin > kRadTolerance*0.5)
  {
    tolORMin2=(fRmin-kRadTolerance*0.5)*(fRmin-kRadTolerance*0.5);
  }
  else
  {
    tolORMin2 = 0 ;
  }
  tolIRMin2 = (fRmin+kRadTolerance*0.5)*(fRmin+kRadTolerance*0.5) ;
  tolORMax2 = (fRmax+kRadTolerance*0.5)*(fRmax+kRadTolerance*0.5) ;
  tolIRMax2 = (fRmax-kRadTolerance*0.5)*(fRmax-kRadTolerance*0.5) ;

  // Set phi divided flag and precalcs

  if (fDPhi < twopi)
  {
    segPhi = true ;
    hDPhi = 0.5*fDPhi ;    // half delta phi
    cPhi = fSPhi + hDPhi ;

    hDPhiOT = hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi 
    hDPhiIT = hDPhi-0.5*kAngTolerance;

    sinCPhi    = std::sin(cPhi) ;
    cosCPhi    = std::cos(cPhi) ;
    cosHDPhiOT = std::cos(hDPhiOT) ;
    cosHDPhiIT = std::cos(hDPhiIT) ;
  }
  else
  {
    segPhi = false ;
  }

  // Theta precalcs
    
  if (fDTheta < pi )
  {
    segTheta  = true ;
    tolSTheta = fSTheta - kAngTolerance*0.5 ;
    tolETheta = fSTheta + fDTheta + kAngTolerance*0.5 ;
  }
  else
  {
    segTheta = false ;
  }

  // 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 ;
  const G4double  flexRadMaxTolerance = // kRadTolerance;
    std::max(kRadTolerance, fEpsilon * fRmax);

  //  if (c > kRadTolerance*fRmax)
  if (c > flexRadMaxTolerance*fRmax)
  {
    // If outside toleranct boundary of outer G4Sphere
    // [should be std::sqrt(rad2)-fRmax > kRadTolerance*0.5]

    d2 = pDotV3d*pDotV3d - c ;

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

      if (s >= 0 )
      {
        xi   = p.x() + s*v.x() ;
        yi   = p.y() + s*v.y() ;
        rhoi = std::sqrt(xi*xi + yi*yi) ;

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

          if (cosPsi >= cosHDPhiOT)
          {
            if (segTheta)   // 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 (segTheta)    // 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 > tolIRMin2 && pDotV3d < 0 )

    if (rad2 > tolIRMax2 && ( d2 >= flexRadMaxTolerance*fRmax && pDotV3d < 0 ) )
    {
      if (segPhi)
      {
        // 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 (segTheta)
          {
            if ( (pTheta >= tolSTheta + kAngTolerance)
              && (pTheta <= tolETheta - kAngTolerance) )
            {
              return snxt=0;
            }
          }
          else    // strictly inside Theta in both cases
          {
            return snxt=0;
          }
        }
      }
      else
      {
        if ( segTheta )
        {
          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 for 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 >- kRadTolerance*0.5 && pDotV3d >= 0 && rad2 < tolIRMin2 )

    if ( c > -kRadTolerance*0.5 && rad2 < tolIRMin2 && 
         ( d2 < fRmin*kCarTolerance || pDotV3d >= 0 ) )
    {
      if (segPhi)
      {
        // 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 (segTheta)
          {
            if ( (pTheta >= tolSTheta + kAngTolerance)
              && (pTheta <= tolETheta - kAngTolerance) )
            {
              return snxt=0;
            }
          }
          else
          {
            return snxt = 0 ;
          }
        }
      }
      else
      {
        if (segTheta)
        {
          if ( (pTheta >= tolSTheta + kAngTolerance)
            && (pTheta <= tolETheta - kAngTolerance) )
          {
            return snxt = 0 ;
          }
        }
        else
        {
          return snxt=0;
        }
      }
    }
    else   // Not special tolerant case
    {
      //  d2 = pDotV3d*pDotV3d - c ;

      if (d2 >= 0)
      {
        s = -pDotV3d + std::sqrt(d2) ;
        if ( s >= kRadTolerance*0.5 )  // It was >= 0 ??
        {
          xi   = p.x() + s*v.x() ;
          yi   = p.y() + s*v.y() ;
          rhoi = std::sqrt(xi*xi+yi*yi) ;

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

            if (cosPsi >= cosHDPhiOT)
            {
              if (segTheta)  // 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 (segTheta)   // 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 ( segPhi )
  {
    // First phi surface (`S'tarting phi)

    sinSPhi = std::sin(fSPhi) ;
    cosSPhi = std::cos(fSPhi) ;

    // 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 < kCarTolerance*0.5)
      {
        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 ( segTheta )
            {
              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)

    ePhi    = fSPhi + fDPhi ;
    sinEPhi = std::sin(ePhi)     ;
    cosEPhi = std::cos(ePhi)     ;

    // Compnent in outwards normal dirn

    Comp    = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
        
    if (Comp < 0)
    {
      Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
      if ( Dist < kCarTolerance*0.5 )
      {
        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 ( segTheta )
            {
              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 ( segTheta )
  {

    // 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

    tanSTheta  = std::tan(fSTheta)         ;
    tanSTheta2 = tanSTheta*tanSTheta  ;
    tanETheta  = std::tan(fSTheta+fDTheta) ;
    tanETheta2 = tanETheta*tanETheta  ;
      
    if (fSTheta)
    {
      dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
    }
    else
    {
      dist2STheta = kInfinity ;
    }
    if ( fSTheta + fDTheta < pi )
    {
      dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
    }
      else
    {
      dist2ETheta=kInfinity;
    }      
    if ( pTheta < tolSTheta) // dist2STheta<-kRadTolerance*0.5 && dist2ETheta>0)
    {
      // 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 ;
        
      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 ( segPhi && 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 ( fSTheta + fDTheta < pi )
      {
        t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
        t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
        
        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*(fSTheta + fDTheta - halfpi) <= 0) )
            {
              if (segPhi && 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 ;
        
      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*(fSTheta + fDTheta - 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 + fDTheta - halfpi) <= 0) )
          {
            if (segPhi && 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 ;

        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 (segPhi && 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 > kAngTolerance) )
    {
      // 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<pi*.5)
        || (t2<0  && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>pi*.5)
        || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==pi*.5) )
      {
        if (segPhi && 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) ;
      b  = t2/t1 ;
      c  = dist2STheta/t1 ;
      d2 = b*b - c ;

      if (d2 >= 0)
      {
        d = std::sqrt(d2) ;
        s = -b + d ;
        if ( (s >= kCarTolerance*0.5) && (s < snxt) && (fSTheta < pi*0.5) )
        {
          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 ( segPhi && 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)
           && ((fSTheta + fDTheta) < 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    && (fSTheta+fDTheta) <pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
 || (t2>=0   && (fSTheta+fDTheta) >pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
 || (v.z()>0 && (fSTheta+fDTheta)==pi*0.5 && tolIRMin2<rad2 && rad2<tolIRMax2)
         )
      {
        if (segPhi && 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) ;
      b  = t2/t1 ;
      c  = dist2ETheta/t1 ;
      d2 = b*b - c ;

      if (d2 >= 0)
      {
        d = std::sqrt(d2) ;
        s = -b + d ;
        
        if ( (s >= kCarTolerance*0.5)
          && (s < snxt) && ((fSTheta + fDTheta) > pi*0.5) )
        {
          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 + fDTheta - halfpi) <= 0) )
          {
            if (segPhi && 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 ;

      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 (segPhi && 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 ;
        
      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*(fSTheta + fDTheta - halfpi) <= 0) )
          {
            if (segPhi && 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,rad,rho;
  G4double phiC,cosPhiC,sinPhiC,cosPsi,ePhi;
  G4double pTheta,dTheta1,dTheta2;
  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
  //    
  if (fRmin)
  {
    safeRMin=fRmin-rad;
    safeRMax=rad-fRmax;
    if (safeRMin>safeRMax)
    {
      safe=safeRMin;
    }
    else
    {
      safe=safeRMax;
    }
  }
  else
  {
    safe=rad-fRmax;
  }

  //
  // Distance to phi extent
  //
  if (fDPhi<twopi&&rho)
  {
    phiC=fSPhi+fDPhi*0.5;
    cosPhiC=std::cos(phiC);
    sinPhiC=std::sin(phiC);

    // Psi=angle from central phi to point
    //
    cosPsi=(p.x()*cosPhiC+p.y()*sinPhiC)/rho;
    if (cosPsi<std::cos(fDPhi*0.5))
    {
      // Point lies outside phi range
      //
      if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
      {
        safePhi=std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
      }
      else
      {
        ePhi=fSPhi+fDPhi;
        safePhi=std::fabs(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi));
      }
      if (safePhi>safe) safe=safePhi;
    }
  }
  //
  // Distance to Theta extent
  //    
  if ((rad!=0.0) && (fDTheta<pi))
  {
    pTheta=std::acos(p.z()/rad);
    if (pTheta<0) pTheta+=pi;
    dTheta1=fSTheta-pTheta;
    dTheta2=pTheta-(fSTheta+fDTheta);
    if (dTheta1>dTheta2)
    {
      if (dTheta1>=0)             // WHY ???????????
      {
        safeTheta=rad*std::sin(dTheta1);
        if (safe<=safeTheta)
        {
          safe=safeTheta;
        }
      }
    }
    else
    {
      if (dTheta2>=0)
      {
        safeTheta=rad*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;  

  G4double t1,t2;
  G4double b,c,d;

  // Variables for phi intersection:

  G4double sinSPhi,cosSPhi,ePhi,sinEPhi,cosEPhi;
  G4double cPhi,sinCPhi,cosCPhi;
  G4double pDistS,compS,pDistE,compE,sphi2,vphi;
    
  G4double rho2,rad2,pDotV2d,pDotV3d,pTheta;

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

  // G4double Comp; // Phi intersection

  G4bool segPhi;        // Phi flag and precalcs
  G4double hDPhi,hDPhiOT,hDPhiIT; 
  G4double cosHDPhiOT,cosHDPhiIT;

  G4bool segTheta;                             // Theta flag and precals
  G4double tanSTheta=0.,tanETheta=0., rhoSecTheta;
  G4double tanSTheta2=0.,tanETheta2=0.;
  G4double dist2STheta, dist2ETheta, distTheta;
  G4double d2,s;

  // General Precalcs

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

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

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

  // Set phi divided flag and precalcs

  if( fDPhi < twopi )
  {
    segPhi=true;
    hDPhi=0.5*fDPhi;    // half delta phi
    cPhi=fSPhi+hDPhi;;
    hDPhiOT=hDPhi+0.5*kAngTolerance; // Outer Tolerant half delta phi 
    hDPhiIT=hDPhi-0.5*kAngTolerance;
    sinCPhi=std::sin(cPhi);
    cosCPhi=std::cos(cPhi);
    cosHDPhiOT=std::cos(hDPhiOT);
    cosHDPhiIT=std::cos(hDPhiIT);
  }
  else
  {
    segPhi=false;
  }

  // Theta precalcs
    
  if ( fDTheta < pi )
  {
    segTheta  = true;
    tolSTheta = fSTheta - kAngTolerance*0.5;
    tolETheta = fSTheta + fDTheta + kAngTolerance*0.5;
  }
  else segTheta = false;

    
  // 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))
  //
  // const G4double  fractionTolerance = 1.0e-12;

  const G4double  flexRadMaxTolerance =            // kRadTolerance;
    std::max(kRadTolerance, fEpsilon * fRmax);

  const G4double  Rmax_plus = fRmax + flexRadMaxTolerance*0.5;

  const G4double  flexRadMinTolerance = std::max(kRadTolerance, 
                     fEpsilon * fRmin);

  const G4double  Rmin_minus= (fRmin > 0) ? fRmin-flexRadMinTolerance*0.5 : 0 ;

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

    if (c < flexRadMaxTolerance*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 >- flexRadMaxTolerance*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 >- flexRadMinTolerance*fRmin ) // 2.0 * (0.5*kRadTolerance) * fRmin
      {
        if( c < flexRadMinTolerance*fRmin && 
            d2 >= flexRadMinTolerance*fRmin && pDotV3d < 0 ) // leaving from Rmin
        {
          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 (segTheta)
  {
    // 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
    //
  
    /* ////////////////////////////////////////////////////////

    tanSTheta=std::tan(fSTheta);
    tanSTheta2=tanSTheta*tanSTheta;
    tanETheta=std::tan(fSTheta+fDTheta);
    tanETheta2=tanETheta*tanETheta;
      
    if (fSTheta)
    {
      dist2STheta=rho2-p.z()*p.z()*tanSTheta2;
    }
    else
    {
      dist2STheta = kInfinity;
    }
    if (fSTheta + fDTheta < pi)
    {
      dist2ETheta = rho2-p.z()*p.z()*tanETheta2;
    }
    else
    {
      dist2ETheta = kInfinity ;
    }
    if (pTheta > tolSTheta && pTheta < tolETheta)   // Inside theta  
    {
      // In tolerance of STheta and possible leaving out to small thetas N-

      if(pTheta < tolSTheta + kAngTolerance  && fSTheta > kAngTolerance)  
      {
        t2=pDotV2d-p.z()*v.z()*tanSTheta2 ; // =(VdotN+)*rhoSecSTheta

        if( fSTheta < pi*0.5 && t2 < 0)
        {
          if(calcNorm) *validNorm = false ;
          return snxt = 0 ;
        }
        else if(fSTheta > pi*0.5 && t2 >= 0)
        {
          if(calcNorm)
          {
            rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2)) ;
            *validNorm = true ;
            *n = G4ThreeVector(-p.x()/rhoSecTheta,   // N-
                               -p.y()/rhoSecTheta,
                               tanSTheta/std::sqrt(1+tanSTheta2) ) ;
          }
          return snxt = 0 ;
        }
        else if( fSTheta == pi*0.5 && v.z() > 0)
        {
          if(calcNorm)
          {
            *validNorm = true ;
            *n = G4ThreeVector(0,0,1) ;
          }
          return snxt = 0 ;
        }
      }

      // In tolerance of ETheta and possible leaving out to larger thetas N+

      if ( (pTheta  > tolETheta - kAngTolerance)
        && (( fSTheta + fDTheta) < pi - kAngTolerance) )  
      {
        t2=pDotV2d-p.z()*v.z()*tanETheta2 ;
        if((fSTheta+fDTheta)>pi*0.5 && t2<0)
        {
          if(calcNorm) *validNorm = false ;
          return snxt = 0 ;
        }
        else if( (fSTheta+fDTheta) < pi*0.5 && t2 >= 0 )
        {
          if(calcNorm)
          {
            rhoSecTheta = std::sqrt(rho2*(1+tanETheta2)) ;
            *validNorm = true ;
            *n = G4ThreeVector( p.x()/rhoSecTheta,  // N+
                                p.y()/rhoSecTheta,
                                -tanETheta/std::sqrt(1+tanETheta2)  ) ; 
          }
          return snxt = 0 ;
        }
        else if( ( fSTheta+fDTheta) == pi*0.5 && v.z() < 0 )
        {
          if(calcNorm)
          {
            *validNorm = true ;
            *n = G4ThreeVector(0,0,-1) ;
          }
          return snxt = 0 ;
        }
      }
      if( fSTheta > 0 )
      {       
        // First root of fSTheta cone, second if first root -ve

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

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

          if ( s < 0 )
          {
            s = -b + d ;    // Second root
          }
          if (s > flexRadMaxTolerance*0.5 )   // && s<sr)
          {
            // check against double cone solution
            zi=p.z()+s*v.z();
            if (fSTheta<pi*0.5 && zi<0)
            {
              s = kInfinity ;  // wrong cone
            }
            if (fSTheta>pi*0.5 && zi>0)
            {
              s = kInfinity ;  // wrong cone
            }
            stheta = s ;
            sidetheta = kSTheta ;
          }
        }
      }

      // Possible intersection with ETheta cone  
      
      if (fSTheta + fDTheta < pi)
      {
        t1 = 1-v.z()*v.z()*(1+tanETheta2);
        t2 = pDotV2d-p.z()*v.z()*tanETheta2;        
        b  = t2/t1;
        c  = dist2ETheta/t1;
        d2 = b*b-c ;

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

          if ( s < 0 )
          {
            s=-b+d;    // Second root
          }
          if (s > flexRadMaxTolerance*0.5 && s < stheta )
          {
            // check against double cone solution
            zi=p.z()+s*v.z();
            if (fSTheta+fDTheta<pi*0.5 && zi<0)
            {
              s = kInfinity ;  // wrong cone
            }
            if (fSTheta+fDTheta>pi*0.5 && zi>0)
            {
              s = kInfinity ;  // wrong cone
            }
          }
          if (s < stheta)
          {
            stheta = s ;
            sidetheta = kETheta ;
          }
        }
      }
    }  
    */  ////////////////////////////////////////////////////////////

    if(fSTheta) // intersection with first cons
    {

      tanSTheta = std::tan(fSTheta);

      if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0
      {
        if( v.z() > 0. ) 
        {
          if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 )
          {
            if(calcNorm)
            {
              *validNorm = true;
              *n = G4ThreeVector(0.,0.,1.);
            }
            return snxt = 0 ;
          }  
          // s = -p.z()/v.z();
          stheta    = -p.z()/v.z();
          sidetheta = kSTheta;
        }
      }
      else // kons is not plane 
      {
        tanSTheta2  = tanSTheta*tanSTheta;
        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(std::fabs(dist2STheta/(1+tanSTheta2)));
        distTheta = std::sqrt(rho2)-p.z()*tanSTheta;

        if( std::fabs(t1) < 0.5*kAngTolerance ) // 1st order equation, v parallel to kons
        {
          if( v.z() > 0. )
          {
            if(std::fabs(distTheta) < flexRadMaxTolerance*0.5) // 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.;               
              }
            }
            // s = -0.5*dist2STheta/t2;

            stheta    = -0.5*dist2STheta/t2;
            sidetheta = kSTheta;
          }  
        }
        else   // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve
        {
          if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface
          {
            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) < flexRadMaxTolerance*0.5 && t2 < 0.) ||
                              s < 0.  ||
                  ( s > 0. && p.z() + s*v.z() > 0.)                   ) 
              {
                s = -b + d ; // 2nd root
              }
              if( s >  flexRadMaxTolerance*0.5 && 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) < flexRadMaxTolerance*0.5 && t2 >= 0.) ||
                              s < 0.                                   ||
                  ( s > 0. && p.z() + s*v.z() < 0.)                      )
              {
                s = -b + d ; // 2nd root
              }
              if( s >  flexRadMaxTolerance*0.5 && p.z() + s*v.z() >= 0.)  
              {
                stheta    = s;
                sidetheta = kSTheta;
              }            
            }
          }
        }
      }
    }
    if (fSTheta + fDTheta < pi) // intersection with second cons
    {

      tanETheta = std::tan(fSTheta+fDTheta);

      if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0
      {
        if( v.z() < 0. ) 
        {
          if ( std::fabs( p.z() ) <= flexRadMaxTolerance*0.5 )
          {
            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 
      {
        tanETheta2  = tanETheta*tanETheta;
        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(std::fabs(dist2ETheta/(1+tanETheta2)));
        distTheta = std::sqrt(rho2)-p.z()*tanETheta;

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

            if( s < stheta)
            {
              stheta    = s;
              sidetheta = kETheta;
            }
          }  
        }
        else   // 2nd order equation, 1st root of fSTheta cone, 2nd if 1st root -ve
        {
          if( std::fabs(distTheta) < flexRadMaxTolerance*0.5) // && t2 >= 0.) surface
          {
            if( fSTheta+fDTheta < halfpi && t2 >= 0. ) // leave
            {
              if( calcNorm ) 
              {
                *validNorm = true;
                if (rho2)
                {
                    rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
                   
                    *n = G4ThreeVector( p.x()/rhoSecTheta,   
                                        p.y()/rhoSecTheta,
                                        -std::sin(fSTheta+fDTheta)  );
                }
                else *n = G4ThreeVector(0.,0.,-1.);
              }                           
              return snxt = 0.;
            }
            else if( fSTheta+fDTheta > 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( fSTheta+fDTheta < halfpi )
            {
              s = -b - d;         // First root

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

  } // end theta intersections

  // Phi Intersection
    
  if ( fDPhi < twopi)
  {
    sinSPhi=std::sin(fSPhi);
    cosSPhi=std::cos(fSPhi);
    ePhi=fSPhi+fDPhi;
    sinEPhi=std::sin(ePhi);
    cosEPhi=std::cos(ePhi);
    cPhi=fSPhi+fDPhi*0.5;
    sinCPhi=std::sin(cPhi);
    cosCPhi=std::cos(cPhi);

    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 intersecting with correct half-plane
          // (if not -> no intersect)

          if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
          {
            sphi=kInfinity;
          }
          else
          {
            sidephi = kSPhi ;
            if ( pDistS > -0.5*kCarTolerance) 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 intersecting with correct half-plane
 
            if ((yi*cosCPhi-xi*sinCPhi)>=0) // Leaving via ending phi
            {
              sidephi = kEPhi ;
              if ( pDistE <= -0.5*kCarTolerance )
              {
                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 ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
            {
              sphi = kInfinity ;
            }
            else // Leaving via Ending phi
            {
              sidephi = kEPhi ;
              if ( pDistE > -0.5*kCarTolerance ) 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 ( ( 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 ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
            {
              sphi = kInfinity ;
            }
           else  // Leaving via Starting phi
           {
              sidephi = kSPhi ;
             if ( pDistS > -0.5*kCarTolerance ) 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 ( ( 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 < vphi && vphi < fSPhi + fDPhi )
        {
          sphi=kInfinity;
        }
        else
        {
          sidephi = kSPhi ; // arbitrary 
          sphi    = 0     ;
        }
      }
      else  // travel along z - no phi intersaction
      {
        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(std::sin(fSPhi),-std::cos(fSPhi),0);
          *validNorm=true;
        }
        else *validNorm=false;
        break ;

      case kEPhi:
        if ( fDPhi <= pi )      // Normal to Phi+
        {
          *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),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();
          rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanSTheta2));
          *n = G4ThreeVector( xi/rhoSecTheta,   // N-
                              yi/rhoSecTheta,
                              -tanSTheta/std::sqrt(1+tanSTheta2));
          *validNorm=true;
        }
        else *validNorm=false;  // Concave STheta cone
        break;

      case kETheta:
        if( ( fSTheta + fDTheta ) == halfpi )
        {
          *n         = G4ThreeVector(0.,0.,-1.);
          *validNorm = true;
        }
        else if ( ( fSTheta + fDTheta ) < halfpi)
        {
          xi=p.x()+snxt*v.x();
          yi=p.y()+snxt*v.y();
          rhoSecTheta = std::sqrt((xi*xi+yi*yi)*(1+tanETheta2));
          *n = G4ThreeVector( xi/rhoSecTheta,   // N+
                              yi/rhoSecTheta,
                              -tanETheta/std::sqrt(1+tanETheta2) );
          *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;
}

/////////////////////////////////////////////////////////////////////////
//
// Calcluate 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,rad,rho;
  G4double phiC,cosPhiC,sinPhiC,ePhi;
  G4double pTheta,dTheta1,dTheta2;
  rho2=p.x()*p.x()+p.y()*p.y();
  rad=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=rad-fRmin;
    safeRMax=fRmax-rad;
    if (safeRMin<safeRMax)
    {
      safe=safeRMin;
    }
    else
    {
      safe=safeRMax;
    }
  }
  else
  {
    safe=fRmax-rad;
  }

  //
  // Distance to phi extent
  //
  if (fDPhi<twopi && rho)
  {
    phiC=fSPhi+fDPhi*0.5;
    cosPhiC=std::cos(phiC);
    sinPhiC=std::sin(phiC);
    if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
    {
      safePhi=-(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
    }
    else
    {
      ePhi=fSPhi+fDPhi;
      safePhi=(p.x()*std::sin(ePhi)-p.y()*std::cos(ePhi));
    }
    if (safePhi<safe) safe=safePhi;
  }

  //
  // Distance to Theta extent
  //    
  if (rad)
  {
    pTheta=std::acos(p.z()/rad);
    if (pTheta<0) pTheta+=pi;
    dTheta1=pTheta-fSTheta;
    dTheta2=(fSTheta+fDTheta)-pTheta;
    if (dTheta1<dTheta2)
    {
      safeTheta=rad*std::sin(dTheta1);
      if (safe>safeTheta)
      {
        safe=safeTheta;
      }
    }
    else
    {
      safeTheta=rad*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 (fDPhi==pi*2.0 && fSPhi==0)
  {
    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 (fDTheta==pi && fSTheta==0)
  {
    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)*std::cos(fSTheta);
  height2 = (fRmax-fRmin)*std::cos(fSTheta+fDTheta);
  slant1  = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta))
                      + height1*height1);
  slant2  = std::sqrt(sqr((fRmax - fRmin)*std::sin(fSTheta+fDTheta))
                      + height2*height2);
  rRand   = RandFlat::shoot(fRmin,fRmax);
  
  aOne = fRmax*fRmax*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta));
  aTwo = fRmin*fRmin*fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta));
  aThr = fDPhi*((fRmax + fRmin)*std::sin(fSTheta))*slant1;
  aFou = fDPhi*((fRmax + fRmin)*std::sin(fSTheta+fDTheta))*slant2;
  aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
  
  phi = RandFlat::shoot(fSPhi, fSPhi + fDPhi); 
  cosphi = std::cos(phi); 
  sinphi = std::sin(phi);
  theta = RandFlat::shoot(fSTheta,fSTheta+fDTheta);
  costheta = std::cos(theta);
  sintheta = std::sqrt(1.-sqr(costheta));

  if( ((fSPhi==0) && (fDPhi==2.*pi)) || (fDPhi==2.*pi) ) {aFiv = 0;}
  if(fSTheta == 0)  {aThr=0;}
  if(fDTheta + fSTheta == pi) {aFou = 0;}
  if(fSTheta == 0.5*pi) {aThr = pi*(fRmax*fRmax-fRmin*fRmin);}
  if(fSTheta + fDTheta == 0.5*pi) { 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 != 0.5*pi)
    {
      zRand = RandFlat::shoot(fRmin*std::cos(fSTheta),fRmax*std::cos(fSTheta));
      return G4ThreeVector(std::tan(fSTheta)*zRand*cosphi,
                           std::tan(fSTheta)*zRand*sinphi,zRand);
    }
    else
    {
      return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
    }    
  }
  else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
  {
    if(fSTheta + fDTheta != 0.5*pi)
    {
      zRand = RandFlat::shoot(fRmin*std::cos(fSTheta+fDTheta),
                              fRmax*std::cos(fSTheta+fDTheta));
      return G4ThreeVector  (std::tan(fSTheta+fDTheta)*zRand*cosphi,
                             std::tan(fSTheta+fDTheta)*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*std::cos(fSPhi),
                         rRand*sintheta*std::sin(fSPhi),rRand*costheta);
  }
  else
  {
    return G4ThreeVector(rRand*sintheta*std::cos(fSPhi+fDPhi),
                         rRand*sintheta*std::sin(fSPhi+fDPhi),rRand*costheta);
  }
}

/////////////////////////////////////////////////////////////////////////////
//
// 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!!!
}