source: trunk/source/geometry/solids/BREPS/src/G4CylindricalSurface.cc @ 1337

//
// ********************************************************************
// * 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: G4CylindricalSurface.cc,v 1.8 2006/06/29 18:42:08 gunter Exp $
// GEANT4 tag $Name: geant4-09-04-beta-01 $
//
// ----------------------------------------------------------------------
// GEANT 4 class source file
//
// G4CylindricalSurface.cc
//
// ----------------------------------------------------------------------

#include "G4CylindricalSurface.hh"
#include "G4Sort.hh"
#include "G4Globals.hh"

G4CylindricalSurface::G4CylindricalSurface() : G4Surface()
{  
  // default constructor
  // default axis is ( 1.0, 0.0, 0.0 ), default radius is 1.0
  axis = G4Vector3D( 1.0, 0.0, 0.0 );
  radius = 1.0;
}


G4CylindricalSurface::G4CylindricalSurface( const G4Vector3D& o, 
                                            const G4Vector3D& a,
                                            G4double r ) //: G4Surface( o )
{ 
  // Normal constructor
  // require axis to be a unit vector
  G4double amag = a.mag();
  if ( amag != 0.0 )
    axis = a * (1/ amag);  // this makes the axis a unit vector
  else 
  {
    G4cerr << "Error in G4CylindricalSurface::G4CylindricalSurface--axis "
           <<"has zero length\n"
           << "\tDefault axis ( 1.0, 0.0, 0.0 ) is used.\n";
    
    axis = G4Vector3D( 1.0, 0.0, 0.0 );
  }

  //  Require radius to be non-negative
  if ( r >= 0.0 )
    radius = r;
  else 
  {
    G4cerr << "Error in G4CylindricalSurface::G4CylindricalSurface"
           << "--asked for negative radius\n"
           << "\tDefault radius of 1.0 is used.\n";
  
    radius = 1.0;
  }
        
  origin =o;
}


G4CylindricalSurface::~G4CylindricalSurface()
{
}

/*
G4CylindricalSurface::G4CylindricalSurface( const G4CylindricalSurface& c )
  : G4Surface( c.origin )
{
  axis = c.axis;  radius = c.radius;
}
*/

const char* G4CylindricalSurface::NameOf() const
{
  return "G4CylindricalSurface";
}

void G4CylindricalSurface::PrintOn( std::ostream& os ) const
{
  // printing function using C++ std::ostream class
  os << "G4CylindricalSurface surface with origin: " << origin << "\t"
     << "radius: " << radius << "\tand axis " << axis << "\n";
}

//G4int G4Surface::Intersect(const G4Ray& ry)
G4int G4CylindricalSurface::Intersect(const G4Ray& ry)
{

  //  L. Broglia : copy of G4FCylindricalSurface::Intersect
  
  //  Distance along a Ray (straight line with G4ThreeVec) to leave or enter
  //  a G4CylindricalSurface.  The input variable which_way should be set 
  //  to +1 to indicate leaving a G4CylindricalSurface, -1 to indicate 
  //  entering a G4CylindricalSurface.
  //  p is the point of intersection of the Ray with the G4CylindricalSurface.
  //  If the G4Vector3D of the Ray is opposite to that of the Normal to
  //  the G4CylindricalSurface at the intersection point, it will not leave 
  //  the G4CylindricalSurface.
  //  Similarly, if the G4Vector3D of the Ray is along that of the Normal 
  //  to the G4CylindricalSurface at the intersection point, it will not enter
  //  the G4CylindricalSurface.
  //  This method is called by all finite shapes sub-classed to 
  //  G4CylindricalSurface.
  //  Use the virtual function table to check if the intersection point
  //  is within the boundary of the finite shape.
  //  A negative result means no intersection.
  //  If no valid intersection point is found, set the distance
  //  and intersection point to large numbers.

  // G4int which_way = -1; 
  //Originally a parameter.Read explanation above. 
  
  G4int which_way=1;
  
  if(!Inside(ry.GetStart()))
    which_way = -1;        
  
  distance = FLT_MAXX;
  G4Vector3D lv ( FLT_MAXX, FLT_MAXX, FLT_MAXX );
  
  closest_hit = lv;

  //  Origin and G4Vector3D unit vector of Ray.
  G4Vector3D x    = ry.GetStart();
  G4Vector3D dhat = ry.GetDir();

  //  Axis unit vector of the G4CylindricalSurface.
  G4Vector3D ahat = GetAxis();
  G4int isoln   = 0, 
        maxsoln = 2;
  
  //  array of solutions in distance along the Ray
  G4double s[2];
  s[0] = -1.0; 
  s[1] = -1.0 ;

  //  calculate the two solutions (quadratic equation)
  G4Vector3D d    = x - GetOrigin();
  G4double   radiu = GetRadius();

  //quit with no intersection if the radius of the G4CylindricalSurface is zero
  //    if ( radiu <= 0.0 )
  //            return 0;
        
  G4double dsq  = d * d;
  G4double da   = d * ahat;
  G4double dasq = da * da;
  G4double rsq  = radiu * radiu;
  G4double qsq  = dsq - dasq;
  G4double dira = dhat * ahat;
  G4double a    = 1.0 - dira * dira;
  
  if ( a <= 0.0 )
    return 0;
        
  G4double b       = 2. * ( d * dhat - da * dira );
  G4double c       = rsq - qsq;
  G4double radical = b * b + 4. * a * c; 
        
  if ( radical < 0.0 ) 
    return 0;

  G4double root = std::sqrt( radical );
  s[0] = ( - b + root ) / ( 2. * a );
  s[1] = ( - b - root ) / ( 2. * a );

  //  order the possible solutions by increasing distance along the Ray
  //  (G4Sorting routines are in support/G4Sort.h)
  sort_double( s, isoln, maxsoln-1 );

  //  now loop over each positive solution, keeping the first one (smallest
  //  distance along the Ray) which is within the boundary of the sub-shape
  //  and which also has the correct G4Vector3D with respect to the Normal to
  //  the G4CylindricalSurface at the intersection point
  for ( isoln = 0; isoln < maxsoln; isoln++ ) 
  {
    if ( s[isoln] >= kCarTolerance*0.5 ) 
    {
      if ( s[isoln] >= FLT_MAXX )  // quit if too large
        return 0;
      
      distance = s[isoln];
      closest_hit = ry.GetPoint( distance );
      G4double  tmp =  dhat * (Normal( closest_hit ));
      
      if ((tmp * which_way) >= 0.0 )
        if ( WithinBoundary( closest_hit ) == 1 )
          distance =  distance*distance;
      
      return 1;
    }
  }
  
  //  get here only if there was no solution within the boundary, Reset
  //  distance and intersection point to large numbers
  distance = FLT_MAXX;
  closest_hit = lv;


  return 0;
}




G4double G4CylindricalSurface::HowNear( const G4Vector3D& x ) const
{
  //  Distance from the point x to the infinite G4CylindricalSurface.
  //  The distance will be positive if the point is Inside the 
  //  G4CylindricalSurface, negative if the point is outside.
  //  Note that this may not be correct for a bounded cylindrical object
  //  subclassed to G4CylindricalSurface.
  G4Vector3D d = x - origin;
  G4double dA = d * axis;
  G4double rad = std::sqrt( d.mag2() - dA*dA );
  G4double hownear = std::fabs( radius - rad );
 
  return hownear;
}


/*
G4double G4CylindricalSurface::distanceAlongRay( G4int which_way, const G4Ray* ry,
                                   G4Vector3D& p ) const 
{  //  Distance along a Ray (straight line with G4Vector3D) to leave or enter
   //  a G4CylindricalSurface.  The input variable which_way should be set to +1
   //  to indicate leaving a G4CylindricalSurface, -1 to indicate entering a
   //  G4CylindricalSurface.
   //  p is the point of intersection of the Ray with the G4CylindricalSurface.
   //  If the G4Vector3D of the Ray is opposite to that of the Normal to
   //  the G4CylindricalSurface at the intersection point, it will not leave the
   //  G4CylindricalSurface.
   //  Similarly, if the G4Vector3D of the Ray is along that of the Normal 
   //  to the G4CylindricalSurface at the intersection point, it will not enter
   //  the G4CylindricalSurface.
   //  This method is called by all finite shapes sub-classed to
   //  G4CylindricalSurface.
   //  Use the virtual function table to check if the intersection point
   //  is within the boundary of the finite shape.
   //  A negative result means no intersection.
   //  If no valid intersection point is found, set the distance
   //  and intersection point to large numbers.
        G4double Dist = FLT_MAXX;
        G4Vector3D lv ( FLT_MAXX, FLT_MAXX, FLT_MAXX );
        p = lv;
//  Origin and G4Vector3D unit vector of Ray.
        G4Vector3D x = ry->Position();
        G4Vector3D dhat = ry->Direction( 0.0 );
//  Axis unit vector of the G4CylindricalSurface.
        G4Vector3D ahat = GetAxis();
        G4int isoln = 0, maxsoln = 2;
//  array of solutions in distance along the Ray
//      G4double s[2] = { -1.0, -1.0 };
        G4double s[2];s[0] = -1.0; s[1]= -1.0 ;
//  calculate the two solutions (quadratic equation)
        G4Vector3D d = x - GetOrigin();
        G4double radius = GetRadius();
//  quit with no intersection if the radius of the G4CylindricalSurface is zero
        if ( radius <= 0.0 )
                return Dist;
        G4double dsq = d * d;
        G4double da = d * ahat;
        G4double dasq = da * da;
        G4double rsq = radius * radius;
        G4double qsq = dsq - dasq;
        G4double dira = dhat * ahat;
        G4double a = 1.0 - dira * dira;
        if ( a <= 0.0 )
                return Dist;
        G4double b = 2. * ( d * dhat - da * dira );
        G4double c = rsq - qsq;
        G4double radical = b * b + 4. * a * c; 
        if ( radical < 0.0 ) 
                return Dist;
        G4double root = std::sqrt( radical );
        s[0] = ( - b + root ) / ( 2. * a );
        s[1] = ( - b - root ) / ( 2. * a );
//  order the possible solutions by increasing distance along the Ray
//  (G4Sorting routines are in support/G4Sort.h)
        G4Sort_double( s, isoln, maxsoln-1 );
//  now loop over each positive solution, keeping the first one (smallest
//  distance along the Ray) which is within the boundary of the sub-shape
//  and which also has the correct G4Vector3D with respect to the Normal to
//  the G4CylindricalSurface at the intersection point
        for ( isoln = 0; isoln < maxsoln; isoln++ ) {
                if ( s[isoln] >= 0.0 ) {
                        if ( s[isoln] >= FLT_MAXX )  // quit if too large
                                return Dist;
                        Dist = s[isoln];
                        p = ry->Position( Dist );
                        if ( ( ( dhat * Normal( p ) * which_way ) >= 0.0 )
                          && ( WithinBoundary( p ) == 1 ) )
                                return Dist;
                }
        }
//  get here only if there was no solution within the boundary, Reset
//  distance and intersection point to large numbers
        p = lv;
        return FLT_MAXX;
}

*/


/*


G4double G4CylindricalSurface::distanceAlongHelix( G4int which_way,
                                   const Helix* hx,
                                   G4Vector3D& p ) const 
{  //  Distance along a Helix to leave or enter a G4CylindricalSurface.
   //  The input variable which_way should be set to +1 to
   //  indicate leaving a G4CylindricalSurface, -1 to indicate entering a
   //  G4CylindricalSurface.
   //  p is the point of intersection of the Helix with the G4CylindricalSurface.
   //  If the G4Vector3D of the Helix is opposite to that of the Normal to
   //  the G4CylindricalSurface at the intersection point, it will not leave the
   //  G4CylindricalSurface.
   //  Similarly, if the G4Vector3D of the Helix is along that of the Normal 
   //  to the G4CylindricalSurface at the intersection point, it will not enter
   //  the G4CylindricalSurface.
   //  This method is called by all finite shapes sub-classed to
   //  G4CylindricalSurface.
   //  Use the virtual function table to check if the intersection point
   //  is within the boundary of the finite shape.
   //  If no valid intersection point is found, set the distance
   //  and intersection point to large numbers.
   //  Possible negative distance solutions are discarded.
        G4double Dist = FLT_MAXX;
        G4Vector3D lv ( FLT_MAXX, FLT_MAXX, FLT_MAXX );
        G4Vector3D zerovec;  // zero G4Vector3D
        p = lv;
        G4int isoln = 0, maxsoln = 4;
//  Array of solutions in turning angle
//      G4double s[4] = { -1.0, -1.0, -1.0, -1.0 };
        G4double s[4];s[0]=-1.0;s[1]= -1.0;s[2]= -1.0;s[3]= -1.0;


//  Flag set to 1 if exact solution is found
        G4int exact = 0;
//  Helix parameters
        G4double rh = hx->GetRadius();          // radius of Helix
        G4Vector3D ah = hx->GetAxis();          // axis of Helix
        G4Vector3D oh = hx->position();  // origin of Helix
        G4Vector3D dh = hx->direction( 0.0 ); // initial G4Vector3D of Helix
        G4Vector3D prp = hx->getPerp(); // perpendicular vector
        G4double prpmag = prp.Magnitude();
        G4double rhp = rh / prpmag;
//  G4CylindricalSurface parameters
        G4double rc = GetRadius();  // radius of G4CylindricalSurface
        if ( rc == 0.0 )   // quit if zero radius
                return Dist;
        G4Vector3D oc = GetOrigin();  // origin of G4CylindricalSurface
        G4Vector3D ac = GetAxis();  // axis of G4CylindricalSurface
//
//  Calculate quantities of use later on.
        G4Vector3D alpha = rhp * prp;
        G4Vector3D beta = rhp * dh;
        G4Vector3D gamma = oh - oc;
//  Declare variables used later on in several places.
        G4double rcd2 = 0.0, alpha2 = 0.0;
        G4double A = 0.0, B = 0.0, C = 0.0, F = 0.0, G = 0.0, H = 0.0;
        G4double CoverB = 0.0, radical = 0.0, root = 0.0, s1 = 0.0, s2 = 0.0;
        G4Vector3D ghat;
//
//  Set flag for special cases
        G4int special_case = 0;  // 0 means general case
//
//  Test to see if axes of Helix and G4CylindricalSurface are parallel, in which
//  case there are exact solutions.
        if ( ( std::fabs( ah.AngleBetween(ac) ) < FLT_EPSILO )
          || ( std::fabs( ah.AngleBetween(ac) - pi ) < FLT_EPSILO ) ) {
                special_case = 1;
//  If, in addition, gamma is a zero vector or is parallel to the
//  G4CylindricalSurface axis, this simplifies the previous case.
                if ( gamma == zerovec ) {
                        special_case = 3;
                        ghat = gamma;
                }
                else {
                        ghat = gamma / gamma.Magnitude();
                        if ( ( std::fabs( ghat.AngleBetween(ac) ) < FLT_EPSILO )
                          || ( std::fabs( ghat.AngleBetween(ac) - pi ) <
                                                               FLT_EPSILO ) )
                                special_case = 3;
                }
//  Test to see if, in addition to the axes of the Helix and G4CylindricalSurface
//  being parallel, the axis of the G4CylindricalSurface is perpendicular to the
//  initial G4Vector3D of the Helix.
                if ( std::fabs( ( ac * dh ) ) < FLT_EPSILO ) {
//  And, if, in addition to all this, the difference in origins of the Helix
//  and G4CylindricalSurface is  perpendicular to the initial G4Vector3D of the
//  Helix, there is a separate special case.
                        if ( std::fabs( ( ghat * dh ) ) < FLT_EPSILO )
                                special_case = 4;
                }
        }  // end of section with axes of Helix and G4CylindricalSurface parallel
//
//  Another peculiar case occurs if the axis of the G4CylindricalSurface and the
//  initial G4Vector3D of the Helix line up and their origins are the same.
//  This will require a higher order approximation than the general case.
        if ( ( ( std::fabs( dh.AngleBetween(ac) ) < FLT_EPSILO )
            || ( std::fabs( dh.AngleBetween(ac) - pi ) < FLT_EPSILO ) )
            && ( gamma == zerovec ) )
                special_case = 2;
//
//  Now all the special cases have been tagged, so solutions are found
//  for each case.  Exact solutions are indicated by setting exact = 1.
//  [For some reason switch doesn't work here, so do series of if's.]
        if ( special_case == 0  ) {  // approximate quadratic solutions
                        A = beta * beta - ( beta * ac ) * ( beta * ac )
                          + gamma * alpha - ( gamma * ac ) * ( alpha * ac );
                        B = 2.0 * gamma * beta 
                          - 2.0 * ( gamma * ac ) * ( beta * ac );
                        C = gamma * gamma 
                          - ( gamma * ac ) * ( gamma * ac ) - rc * rc;
                        if ( std::fabs( A ) < FLT_EPSILO ) {  // no quadratic term
                                if ( B == 0.0 )  // no intersection, quit
                                        return Dist;
                                else             // B != 0
                                        s[0] = -C / B;
                        }
                        else {  // A != 0, general quadratic solution
                                radical = B * B - 4.0 * A * C;
                                if ( radical < 0.0 )  // no solution, quit 
                                        return Dist;
                                root = std::sqrt( radical );
                                s[0] = ( -B + root ) / ( 2.0 * A );
                                s[1] = ( -B - root ) / ( 2.0 * A );
                                if ( rh < 0.0 ) {
                                        s[0] = -s[0];
                                        s[1] = -s[1];
                                }
                                s[2] = s[0] + 2.0 * pi;
                                s[3] = s[1] + 2.0 * pi;
                        }
        }
//
        else if ( special_case == 1 ) {  // exact solutions
                        exact = 1;
                        H = 2.0 * ( alpha * alpha + gamma * alpha );
                        F = gamma * gamma 
                          - ( ( gamma * ac ) * ( gamma * ac ) )
                          - rc * rc  + H;
                        G = 2.0 * rhp * 
                            ( gamma * dh - ( gamma * ac ) * ( ac * dh ) );
                        A = G * G  + H * H;
                        B = -2.0 * F * H;
                        C = F * F - G * G;
                        if ( std::fabs( A ) < FLT_EPSILO ) {  // no quadratic term
                                if ( B == 0.0 )  // no intersection, quit
                                        return Dist;
                                else {  // B != 0
                                        CoverB = -C / B;
                                        if ( std::fabs( CoverB ) > 1.0 )
                                                return Dist;
                                        s[0] = std::acos( CoverB );
                                }
                        }
                        else {  // A != 0, general quadratic solution
//  Try a different method of calculation using F, G, and H to avoid
//  precision problems.
//                              radical = B * B - 4.0 * A * C;
//                              if ( radical < 0.0 ) {
                                  if ( std::fabs( H ) > FLT_EPSILO ) {
                                    G4double r1 = G / H;
                                    G4double r2 = F / H;
                                    G4double radsq = 1.0 + r1*r1 - r2*r2;
                                    if ( radsq < 0.0 )
                                      return Dist;
                                    root = G * std::sqrt( radsq );
                                    G4double denominator = H * ( 1.0 + r1*r1 );
                                    s1 = ( F + root ) / denominator;
                                    s2 = ( F - root ) / denominator;
                                  }
                                  else
                                    return Dist;
//                              }  //  end radical < 0 condition
//                              else {
//                                root = std::sqrt( radical );
//                                s1 = ( -B + root ) / ( 2.0 * A );
//                                s2 = ( -B - root ) / ( 2.0 * A );
//                              }
                                if ( std::fabs( s1 ) <= 1.0 ) {
                                        s[0] = std::acos( s1 );
                                        s[2] = 2.0 * pi - s[0];
                                }
                                if ( std::fabs( s2 ) <= 1.0 ) {
                                        s[1] = std::acos( s2 );
                                        s[3] = 2.0 * pi - s[1];
                                }
//  Must take only solutions which satisfy original unsquared equation:
//  Gsin(s) - Hcos(s) + F = 0.  Take best solution of pair and set false
//  solutions to -1.  Only do this if the result is significantly different
//  from zero.
                                G4double temp1 = 0.0, temp2 = 0.0;
                                G4double rsign = 1.0;
                                if ( rh < 0.0 ) rsign = -1.0;
                                if ( s[0] > 0.0 ) {
                                   temp1 =  G * rsign * std::sin( s[0] )
                                          - H * std::cos( s[0] ) + F;
                                   temp2 =  G * rsign * std::sin( s[2] )
                                          - H * std::cos( s[2] ) + F;
                                   if ( std::fabs( temp1 ) > std::fabs( temp2 ) )
                                        if ( std::fabs( temp1 ) > FLT_EPSILO ) 
                                          s[0] = -1.0;
                                   else
                                        if ( std::fabs( temp2 ) > FLT_EPSILO ) 
                                          s[2] = -1.0;
                                }
                                if ( s[1] > 0.0 ) {
                                   temp1 =  G * rsign * std::sin( s[1] )
                                          - H * std::cos( s[1] ) + F;
                                   temp2 =  G * rsign * std::sin( s[3] )
                                          - H * std::cos( s[3] ) + F;
                                   if ( std::fabs( temp1 ) > std::fabs( temp2 ) )
                                        if ( std::fabs( temp1 ) > FLT_EPSILO ) 
                                          s[1] = -1.0;
                                   else
                                        if ( std::fabs( temp2 ) > FLT_EPSILO ) 
                                          s[3] = -1.0;
                                }
                        }
        }
//
        else if ( special_case == 2 ) {  // approximate solution
                        G4Vector3D e = ah.cross( ac );
                        G4double re = std::fabs( rhp ) * e.Magnitude();
                        s[0] = std::sqrt( 2.0 * rc / re );
        }
//
        else if ( special_case == 3 ) {  // exact solutions
                        exact = 1;
                        alpha2 = alpha * alpha;
                        rcd2 = rhp * rhp * ( 1.0 - ( (ac*dh) * (ac*dh) ) );
                        A = alpha2 - rcd2;
                        B = - 2.0 * alpha2;
                        C = alpha2 + rcd2 - rc*rc;
                        if ( std::fabs( A ) < FLT_EPSILO ) {  // no quadratic term
                                if ( B == 0.0 )  // no intersection, quit
                                        return Dist;
                                else {  // B != 0
                                        CoverB = -C / B;
                                        if ( std::fabs( CoverB ) > 1.0 )
                                                return Dist;
                                        s[0] = std::acos( CoverB );
                                }
                        }
                        else {  // A != 0, general quadratic solution
                                radical = B * B - 4.0 * A * C;
                                if ( radical < 0.0 ) 
                                        return Dist;
                                root = std::sqrt( radical );
                                s1 = ( -B + root ) / ( 2.0 * A );
                                s2 = ( -B - root ) / ( 2.0 * A );
                                if ( std::fabs( s1 ) <= 1.0 )
                                        s[0] = std::acos( s1 );
                                if ( std::fabs( s2 ) <= 1.0 )
                                        s[1] = std::acos( s2 );
                        }
        }
//
        else if ( special_case == 4 ) {  // exact solution
                        exact = 1;
                        F = gamma * gamma 
                          - ( ( gamma * ac ) * ( gamma * ac ) )
                          - rc * rc;
                        G = 2.0 * ( rhp * rhp  +  gamma * alpha );
                        if ( G == 0.0 )  // no intersection, quit
                                return Dist;
                        G4double cs = 1.0 + ( F / G );
                        if ( std::fabs( cs ) > 1.0 )  // no intersection, quit
                                return Dist;
                        s[0] = std::acos( cs );
        }
//
        else                            // shouldn't get here
                        return Dist;
//
// **************************************************************************
//
//  Order the possible solutions by increasing turning angle
//  (G4Sorting routines are in support/G4Sort.h).
        G4Sort_double( s, isoln, maxsoln-1 );
//
//  Now loop over each positive solution, keeping the first one (smallest
//  distance along the Helix) which is within the boundary of the sub-shape.
        for ( isoln = 0; isoln < maxsoln; isoln++ ) {
                if ( s[isoln] >= 0.0 ) {
//  Calculate distance along Helix and position and G4Vector3D vectors.
                        Dist = s[isoln] * std::fabs( rhp );
                        p = hx->position( Dist );
                        G4Vector3D d = hx->direction( Dist );
                        if ( exact == 0 ) {  //  only for approximate solns
//  Now do approximation to get remaining distance to correct this solution
//  iterate it until the accuracy is below the user-set surface precision.
                           G4double delta = 0.0;  
                           G4double delta0 = FLT_MAXX;
                           G4int dummy = 1;
                           G4int iter = 0;
                           G4int in0 = Inside( hx->position ( 0.0 ) );
                           G4int in1 = Inside( p );
                           G4double sc = Scale();
                           while ( dummy ) {
                                iter++;
//  Terminate loop after 50 iterations and Reset distance to large number,
//  indicating no intersection with G4CylindricalSurface.
//  This generally occurs if the Helix curls too tightly to Intersect it.
                                if ( iter > 50 ) {
                                        Dist = FLT_MAXX;
                                        p = lv;
                                        break;
                                }
//  Find distance from the current point along the above-calculated
//  G4Vector3D using a Ray.
//  The G4Vector3D of the Ray and the Sign of the distance are determined
//  by whether the starting point of the Helix is Inside or outside of
//  the G4CylindricalSurface.
                                in1 = Inside( p );
                                if ( in1 ) {  //  current point Inside
                                   if ( in0 ) {  //  starting point Inside
                                      Ray* r = new Ray( p, d );
                                      delta = 
                                          distanceAlongRay( 1, r, p );
                                      delete r;
                                   }
                                   else {       //  starting point outside
                                      Ray* r = new Ray( p, -d );
                                      delta = 
                                         -distanceAlongRay( 1, r, p );
                                      delete r;
                                   }
                                }
                                else {        //  current point outside
                                   if ( in0 ) {  //  starting point Inside
                                      Ray* r = new Ray( p, -d );
                                      delta = 
                                         -distanceAlongRay( -1, r, p );
                                      delete r;
                                   }
                                   else {        //  starting point outside
                                      Ray* r = new Ray( p, d );
                                      delta = 
                                          distanceAlongRay( -1, r, p );
                                      delete r;
                                   }
                                }
//  Test if distance is less than the surface precision, if so Terminate loop.
                                if ( std::fabs( delta / sc ) <= SURFACE_PRECISION )
                                        break;
//  If delta has not changed sufficiently from the previous iteration, 
//  skip out of this loop.
                                if ( std::fabs( ( delta - delta0 ) / sc ) <=
                                                           SURFACE_PRECISION )
                                        break;
//  If delta has increased in absolute value from the previous iteration
//  either the Helix doesn't Intersect the G4CylindricalSurface or the approximate
//  solution is too far from the real solution.  Try groping for a solution.
//  If not found, Reset distance to large number, indicating no intersection with
//  the G4CylindricalSurface.
                                if ( std::fabs( delta ) > std::fabs( delta0 ) ) {
                                        Dist = std::fabs( rhp ) * 
                                               gropeAlongHelix( hx );
                                        if ( Dist < 0.0 ) {
                                                Dist = FLT_MAXX;
                                                p = lv;
                                        }
                                        else
                                                p = hx->position( Dist );
                                        break;
                                }
//  Set old delta to new one.
                                delta0 = delta;
//  Add distance to G4CylindricalSurface to distance along Helix.
                                Dist += delta;
//  Negative distance along Helix means Helix doesn't Intersect
//  G4CylindricalSurface.
//  Reset distance to large number, indicating no intersection with
//  G4CylindricalSurface.
                                if ( Dist < 0.0 ) {
                                        Dist = FLT_MAXX;
                                        p = lv;
                                        break;
                                }
//  Recalculate point along Helix and the G4Vector3D.
                                p = hx->position( Dist );
                                d = hx->direction( Dist );
                           }  //  end of while loop
                        }  //  end of exact == 0 condition
//  Now have best value of distance along Helix and position for this
//  solution, so test if it is within the boundary of the sub-shape
//  and require that it point in the correct G4Vector3D with respect to
//  the Normal to the G4CylindricalSurface.
                        if ( ( Dist < FLT_MAXX ) &&
                             ( ( hx->direction( Dist ) * Normal( p ) *
                                 which_way ) >= 0.0 ) &&
                             ( WithinBoundary( p ) == 1 ) )
                                return Dist;
                }  // end of if s[isoln] >= 0.0 condition
        }  // end of for loop over solutions
//  if one gets here, there is no solution, so set distance along Helix
//  and position to large numbers
        Dist = FLT_MAXX;
        p = lv;
        return Dist;
}
*/


G4Vector3D G4CylindricalSurface::Normal( const G4Vector3D& p ) const
{ 
  //  return the Normal unit vector to the G4CylindricalSurface 
  //  at a point p on (or nearly on) the G4CylindricalSurface
 
  G4Vector3D n = ( p - origin ) - ( ( p - origin ) * axis ) * axis;
  G4double nmag = n.mag();
 
  if ( nmag != 0.0 )
    n = n * (1/nmag);
  
  return n;
}


G4Vector3D G4CylindricalSurface::SurfaceNormal( const G4Point3D& p ) const
{
  //  return the Normal unit vector to the G4CylindricalSurface at a point 
  //  p on (or nearly on) the G4CylindricalSurface
  
  G4Vector3D n = ( p - origin ) - ( ( p - origin ) * axis ) * axis;
  G4double nmag = n.mag();
  
  if ( nmag != 0.0 )
    n = n * (1/nmag);
  
  return n;
}


G4int G4CylindricalSurface::Inside ( const G4Vector3D& x ) const
{ 
  //  Return 0 if point x is outside G4CylindricalSurface, 1 if Inside.
  //  Outside means that the distance to the G4CylindricalSurface would 
  //  be negative.
  //  Use the HowNear function to calculate this distance.
  if ( HowNear( x ) >= -0.5*kCarTolerance )
    return 1;
  else
    return 0; 
}


G4int G4CylindricalSurface::WithinBoundary( const G4Vector3D& x ) const
{ 
  //  return 1 if point x is on the G4CylindricalSurface, otherwise return zero
  //  base this on the surface precision factor set in support/globals.h
  if ( std::fabs( HowNear( x ) / Scale() ) <= SURFACE_PRECISION )
    return 1;
  else
    return 0;
}


G4double G4CylindricalSurface::Scale() const
{
  //  Returns the radius of a G4CylindricalSurface unless it is zero, in which
  //  case returns the arbitrary number 1.0.
  //  This is ok since derived finite-sized classes will overwrite this.
  //  Used for Scale-invariant tests of surface thickness.
  if ( radius == 0.0 )
    return 1.0;
  else
    return radius;
}


//void G4CylindricalSurface::rotate( G4double alpha, G4double beta,
//                     G4double gamma, G4ThreeMat& m, G4int inverse )
// //  rotate G4CylindricalSurface  first about global x-axis by angle alpha,
//                  second about global y-axis by angle beta,
//               and third about global z-axis by angle gamma
//  by creating and using G4ThreeMat objects in Surface::rotate
//  angles are assumed to be given in radians
//  if inverse is non-zero, the order of rotations is reversed
//  the axis is rotated here, the origin is rotated by calling
//  Surface::rotate
//      G4Surface::rotate( alpha, beta, gamma, m, inverse );
//      axis = m * axis;
//}

//void G4CylindricalSurface::rotate( G4double alpha, G4double beta,
//                     G4double gamma, G4int inverse )
//{  //  rotate G4CylindricalSurface  first about global x-axis by angle alpha,
//                  second about global y-axis by angle beta,
//               and third about global z-axis by angle gamma
//  by creating and using G4ThreeMat objects in Surface::rotate
//  angles are assumed to be given in radians
//  if inverse is non-zero, the order of rotations is reversed
//  the axis is rotated here, the origin is rotated by calling
//  Surface::rotate
//      G4ThreeMat m;
//      G4Surface::rotate( alpha, beta, gamma, m, inverse );
//      axis = m * axis;
//}


void G4CylindricalSurface::SetRadius( G4double r )
{ 
  //  Reset the radius of the G4CylindricalSurface
  //  Require radius to be non-negative
  if ( r >= 0.0 )
    radius = r;
  //  use old value (do not change radius) if out of the range, 
  //  but Print message
  else 
  {
    G4cerr << "Error in G4CylindricalSurface::SetRadius"
           << "--asked for negative radius\n"
           << "\tDefault radius of " << radius << " is used.\n";
  }
}


/*
G4double G4CylindricalSurface::gropeAlongHelix( const Helix* hx ) const
{  //  Grope for a solution of a Helix intersecting a G4CylindricalSurface.
   //  This function returns the turning angle (in radians) where the
   //  intersection occurs with only positive values allowed, or -1.0 if
   //  no intersection is found.
//  The idea is to start at the beginning of the Helix, then take steps
//  of some fraction of a turn.  If at the end of a Step, the current position
//  along the Helix and the previous position are on opposite sides of the
//  G4CylindricalSurface, then the solution must lie somewhere in between.
        G4int one_over_f = 8;  // one over fraction of a turn to go in each Step
        G4double turn_angle = 0.0;
        G4double dist_along = 0.0;
        G4double d_new;
        G4double fk = 1.0 / G4double( one_over_f );
        G4double scal = Scale();
        G4double d_old = HowNear( hx->position( dist_along ) );
        G4double rh = hx->GetRadius();  // radius of Helix
        G4Vector3D prp = hx->getPerp(); // perpendicular vector
        G4double prpmag = prp.Magnitude();
        G4double rhp = rh / prpmag;
        G4int max_iter = one_over_f * HELIX_MAX_TURNS;
//  Take up to a user-settable number of turns along the Helix,
//  groping for an intersection point.
        for ( G4int k = 1; k < max_iter; k++ ) {
                turn_angle = 2.0 * pi * k / one_over_f;
                dist_along = turn_angle * std::fabs( rhp );
                d_new = HowNear( hx->position( dist_along ) );
                if ( ( d_old < 0.0 && d_new > 0.0 ) ||
                     ( d_old > 0.0 && d_new < 0.0 ) ) {
                        d_old = d_new;
//  Old and new points are on opposite sides of the G4CylindricalSurface, therefore
//  a solution lies in between, use a binary search to pin the point down
//  to the surface precision, but don't do more than 50 iterations.
                        G4int itr = 0;
                        while ( std::fabs( d_new / scal ) > SURFACE_PRECISION ) {
                                itr++;
                                if ( itr > 50 )
                                        return turn_angle;
                                turn_angle -= fk * pi;
                                dist_along = turn_angle * std::fabs( rhp );
                                d_new = HowNear( hx->position( dist_along ) );
                                if ( ( d_old < 0.0 && d_new > 0.0 ) ||
                                     ( d_old > 0.0 && d_new < 0.0 ) )
                                        fk *= -0.5;
                                else
                                        fk *=  0.5;
                                d_old = d_new;
                        }  //  end of while loop
                        return turn_angle;  // this is the best solution
                }  //  end of if condition
        }  //  end of for loop
//  Get here only if no solution is found, so return -1.0 to indicate that.
        return -1.0;
}
*/