//----------------------------------------------------------------------------
//  EgalTech 2017-2022
//----------------------------------------------------------------------------
// File : IntersLineSurfStd.cpp      Data : 08.01.22        Versione : 2.4a3
// Contenuto : Implementazione delle funzioni di intersezione
//             componente lineare e superficie standard.
//
//
// Modifiche : 22.01.15 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

#include "stdafx.h"
#include "IntersLineSurfStd.h"
#include "IntersLineBox.h"
#include "/EgtDev/Include/EGkIntersLinePlane.h"
#include "/EgtDev/Include/EGkFrame3d.h"
#include "/EgtDev/Include/ENkPolynomialRoots.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include <algorithm>

using namespace std ;

//----------------------------------------------------------------------------
int
LineDisc( const Point3d& ptPLine, const Vector3d& vtVLine,
          const Point3d& ptPDisc, const Vector3d& vtVDisc, double dRad, 
          double& dU1, double& dU2) 
{
  // Il raggio del disco deve essere non nullo.
   if ( dRad < EPS_SMALL)
      return D_ERROR_INT ;
  // Si richiede che i vettori siano normalizzati
   if ( ! vtVDisc.IsNormalized()  ||  ! vtVLine.IsNormalized())
      return D_ERROR_INT ; 

  // Intersezione con il piano del disco
   Point3d ptInt ;
   Plane3d plDisc ; plDisc.Set( ptPDisc, vtVDisc) ;
   int nLpInt = IntersLinePlane( ptPLine, vtVLine, 1, plDisc, ptInt, false) ;
   if ( nLpInt == ILPT_NO)
      return D_NO_INTERS ;

  // Un punto di intersezione con il piano
   if ( nLpInt == ILPT_YES) {
      double dDistCdP = Dist( ptInt, ptPDisc) ;
     // Se distanza uguale al raggio, intersezione sul bordo 
      if ( abs( dDistCdP - dRad) < EPS_SMALL) {
         dU1 = ( ptInt - ptPLine) * vtVLine ;
         dU2 = dU1 ;
         return D_BOUNDARY_INT_LINE_NOT_IN_PLANE ;
      }
     // Se distanza superiore al raggio, nessuna intersezione
      if ( dDistCdP > dRad)
         return D_NO_INTERS ;
     // Distanza inferiore al raggio, intersezione interna
      dU1 = ( ptInt - ptPLine) * vtVLine ;
      dU2 = dU1 ;
      return D_INNER_INT_LINE_NOT_IN_PLANE ;
   }

  // La linea giace nel piano
   {
     // Proiezione del centro del disco sulla linea
      Point3d ptP = ptPLine + (( ptPDisc - ptPLine) * vtVLine) * vtVLine ;
     // Distanza di questo punto di proiezione dal centro del disco
      double dDist = Dist( ptPDisc, ptP) ;
     // Se distanza uguale al raggio, intersezione tangente
      if ( abs( dDist - dRad) < EPS_SMALL) {
         dU1 = ( ptP - ptPLine) * vtVLine ;
         dU2 = dU1 ;
         return D_ONE_INT_LINE_ON_PLANE ;
      }
     // Se distanza superiore al raggio, nessuna intersezione
      if ( dDist > dRad)
         return D_NO_INTERS ;
     // Distanza inferiore al raggio, due intersezioni secanti
      double dDist2 = sqrt( dRad * dRad - dDist * dDist) ;
      Point3d ptI1 = ptP - dDist2 * vtVLine ;
      dU1 = ( ptI1 - ptPLine) * vtVLine ;
      Point3d ptI2 = ptP + dDist2 * vtVLine ;
      dU2 = ( ptI2 - ptPLine) * vtVLine ;
      return D_INFINITE_INT_LINE_ON_PLANE ;
   }
}

//----------------------------------------------------------------------------
int
RayDisc( const Point3d& ptPLine, const Vector3d& vtVLine,
         const Point3d& ptPDisc,  const Vector3d& vtVDisc, double dRad,
         double& dU1, double& dU2)
{
   int nIntType = LineDisc( ptPLine, vtVLine, ptPDisc, vtVDisc, dRad, dU1, dU2) ;
   if ( nIntType == D_ERROR_INT)
      return nIntType ;

  // La retta corrispondente al raggio è nel 
  // piano e attraversa la circonferenza.
   if ( nIntType == D_INFINITE_INT_LINE_ON_PLANE) {
     // Poiché dU2 > dU1, se dU2 < 0 allora dU1 < 0.
     // La semiretta non interferisce col disco.
      if ( dU2 < - EPS_SMALL)
         nIntType = D_NO_INTERS ;
     // Se solo dU1 < 0, c'è intersezione e la 
     // semi-retta parte dall'interno del disco
      else if ( dU1 < - EPS_SMALL) {
         dU1 = 0 ;
         nIntType = D_INFINITE_INT_LINE_ON_PLANE ; 
      }
   } 
  // In tutti gli altri casi c'è un solo punto 
  // di contatto fra retta e disco.
   else { 
     // Se tale punto avviene a un parametro negativo,
     // la semi-retta non interseca il disco.
      if ( dU1 < - EPS_SMALL)
         nIntType = D_NO_INTERS ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int 
SegmentDisc( const Point3d& ptPLine, const Vector3d& vtVLine, double dLen,
             const Point3d& ptPDisc, const Vector3d& vtVDisc, double dRad, 
             double& dU1, double& dU2)
{
   int nIntType = LineDisc( ptPLine, vtVLine, ptPDisc, vtVDisc, dRad, dU1, dU2) ;
   if ( nIntType == D_ERROR_INT)
      return nIntType ;

  // La retta corrispondente al segmento è nel 
  // piano e attraversa la circonferenza.
   if ( nIntType == D_INFINITE_INT_LINE_ON_PLANE) {
     // Interferenza 
      if ( dU1 < dLen + EPS_SMALL  &&  dU2 > - EPS_SMALL) {
         if ( dU1 > dLen - EPS_SMALL)
            nIntType = D_ONE_INT_LINE_ON_PLANE ;
         else if ( dU2 < EPS_SMALL) {
             dU1 = 0 ;   
             nIntType = D_ONE_INT_LINE_ON_PLANE ;
         }
         else {
            Clamp( dU1, 0., dLen) ;
            Clamp( dU2, 0., dLen) ;
         }
      }
     // Il segmento non interagisce
      else
         nIntType = D_NO_INTERS ;
   }
  // Negli altri casi c'è una sola intersezione 
   else if ( dU1 < - EPS_SMALL  ||  dU1 > dLen + EPS_SMALL) 
      nIntType = D_NO_INTERS ;
   return nIntType ;
}

//----------------------------------------------------------------------------
int 
LineSphere( const Point3d& ptPLine, const Vector3d& vtVLine,
            const Point3d& ptSphC, double dRad, 
            double& dU1, double& dU2)
{
  // Il raggio della sfera deve essere non nullo.
   if ( dRad < EPS_SMALL) 
      return S_ERROR_INT ;   
  // Si richiede che il vettore sia normalizzato
   if ( ! vtVLine.IsNormalized()) 
      return S_ERROR_INT ;
        
  // Vettore congiungente il punto della retta con il centro della sfera
   Vector3d vtR = ( ptSphC - ptPLine) - ( ptSphC - ptPLine) * vtVLine * vtVLine ;
  // Quadrato del raggio della sfera
   double dSqRad = dRad * dRad ;
  // Quadrato della distanza punto retta
   double dSqDist = vtR * vtR ;
  // Studio delle soluzioni
   int nIntType ;
  // Retta esterna
   if ( dSqDist > dSqRad + 2 * dRad * EPS_SMALL) 
      nIntType = S_NO_INTERS ;
  // Retta tangente
   else if ( dSqDist > dSqRad - 2 * dRad * EPS_SMALL) {
      dU1 = ( ptSphC - ptPLine) * vtVLine ;
      nIntType = S_ONE_INT_TAN ;
   }
  // Retta secante
   else {
      double dSemiDistSol = sqrt( dSqRad - dSqDist) ;
      dU1 = ( ptSphC - ptPLine) * vtVLine - dSemiDistSol ;
      dU2 = ( ptSphC - ptPLine) * vtVLine + dSemiDistSol ;
      nIntType = S_TWO_INT ; 
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
RaySphere( const Point3d& ptPLine, const Vector3d& vtVLine, 
           const Point3d& ptSphC, double dRad,
           double& dU1, double& dU2)
{
   int nIntType = LineSphere( ptPLine, vtVLine, ptSphC, dRad, dU1, dU2) ;
   if ( nIntType == S_ERROR_INT)
      return nIntType ; 

  // Caso in cui retta associata al raggio è tangente 
   if ( nIntType == S_ONE_INT_TAN) {
     // Se la soluzione è negativa, il raggio 
     // non interferisce con la sfera.
      if ( dU1 < - EPS_SMALL)      
         nIntType = S_NO_INTERS ;   
   }
  // Caso in cui la retta associata al raggio è secante
   else if ( nIntType == S_TWO_INT) {
     // Se dU1 è negativo, vi è solo una soluzione,
     // e in questo caso il raggio non è tangente.
      if ( dU1 < - EPS_SMALL) {
         dU1 = dU2 ;
         nIntType = S_ONE_INT_SEC ;
      }
     // Se anche dU2 è negativo il raggio,
     // non interferisce.
      if ( dU2 < - EPS_SMALL)
         nIntType = S_NO_INTERS ;      
   } 
   return nIntType ;
}

//----------------------------------------------------------------------------
int
SegmentSphere( const Point3d& ptPLine, const Vector3d& vtVLine, double dLen, 
               const Point3d& ptSphC, double dRad, 
               double& dU1, double& dU2)
{
   int nIntType = LineSphere( ptPLine, vtVLine, ptSphC, dRad, dU1, dU2) ;
   if ( nIntType == S_ERROR_INT)
      return nIntType ;

  // La retta associata al segmento è tangente
   if ( nIntType == S_ONE_INT_TAN) {
     // Se il punto è fuori dal segmento,
     // non vi sono soluzioni
      if ( dU1 < - EPS_SMALL  ||  dU1 > dLen + EPS_SMALL)
         nIntType = S_NO_INTERS ; 
   }
  // La retta associata al segmento è secante 
   else if ( nIntType == S_TWO_INT) {
     // Segmento esterno
      if ( dU1 > dLen + EPS_SMALL) 
         nIntType = S_NO_INTERS ;    
      else if ( dU1 > - EPS_SMALL) {
        // Segmento secante con una intersezione
        // In questo caso è valida dU1
         if ( dU2 > dLen + EPS_SMALL)    
            nIntType = S_ONE_INT_SEC ;
      }
      else {
        // Il segmento non interferisce
         if ( dU2 > dLen + EPS_SMALL) {
            nIntType = S_NO_INTERS ; 
         }
        // Segmento secante con una intersezione
        // In questo caso è valida dU2, quindi 
        // assegnamo dU2 a dU1
         else if ( dU2 > - EPS_SMALL) {
            dU1 = dU2 ;
            nIntType = S_ONE_INT_SEC ;
         }
        // Il segmento non interferisce
         else 
            nIntType = S_NO_INTERS ; 
      }
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int 
LineSemiSphere( const Point3d& ptPLine, const Vector3d& vtVLine,
                const Point3d& ptSphC, const Vector3d& vtSSphOrient, double dRad, 
                double& dU1, double& dU2)
{
   int nIntType = LineSphere( ptPLine, vtVLine, ptSphC, dRad, dU1, dU2) ;
   if ( nIntType == S_ERROR_INT  ||  ! vtSSphOrient.IsNormalized())
      return S_ERROR_INT ;
  // Un punto di tangenza
   if ( nIntType == S_ONE_INT_TAN) {
      Point3d ptInt = ptPLine + dU1 * vtVLine ;
      Vector3d vtInt = ptInt - ptSphC ; 
      vtInt.Normalize() ;
     // Se il punto di contatto è nella calotta superiore : scarto la soluzione
      if ( vtInt * vtSSphOrient > EPS_ZERO) 
         nIntType = S_NO_INTERS ;
   }
  // Due punti di secanza
   else if ( nIntType == S_TWO_INT) {
      Point3d ptInt1 = ptPLine + dU1 * vtVLine ;
      Point3d ptInt2 = ptPLine + dU1 * vtVLine ;
      Vector3d vtInt1 = ptInt1 - ptSphC ; 
      Vector3d vtInt2 = ptInt2 - ptSphC ; 
      vtInt1.Normalize() ;
      vtInt2.Normalize() ;
      int nInt = 2 ;
     // Scarto le soluzioni sulla calotta superiore
      if ( vtInt1 * vtSSphOrient > EPS_ZERO) {
         dU1 = dU2 ;
         -- nInt ;
      } 
      if ( vtInt2 * vtSSphOrient > EPS_ZERO)
         -- nInt ;
     // In base al numero di soluzioni valide  
     // aggiorno la tipologia di intersezione. 
      if ( nInt == 1)
         nIntType = S_ONE_INT_SEC ;
      else if ( nInt == 0)
         nIntType = S_NO_INTERS ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
RaySemiSphere( const Point3d& ptPLine,  const Vector3d& vtVLine,
               const Point3d& ptSphC, const Vector3d& vtSSphOrient, double dRad, 
               double& dU1, double& dU2) 
{
   int nIntType = RaySphere( ptPLine, vtVLine, ptSphC, dRad, dU1, dU2) ;
   if ( nIntType == S_ERROR_INT  ||  ( ! vtSSphOrient.IsNormalized()))
      return S_ERROR_INT ;
  // Un Punto in cui la semi-retta interferisce con la semi-sfera
   if ( nIntType == S_ONE_INT_SEC  ||  nIntType == S_ONE_INT_TAN) {
      Point3d ptInt = ptPLine + dU1 * vtVLine ;
      Vector3d vtInt = ptInt - ptSphC ; 
      vtInt.Normalize() ;
     // Se il punto di contatto è nella calotta superiore
     // Scarto la soluzione
      if ( vtInt * vtSSphOrient > EPS_ZERO) 
         nIntType = S_NO_INTERS ;
   }
  // Due punti in cui la semi-retta è secante
   else if ( nIntType == S_TWO_INT) {
      Point3d ptInt1 = ptPLine + dU1 * vtVLine ;
      Point3d ptInt2 = ptPLine + dU1 * vtVLine ;
      Vector3d vtInt1 = ptInt1 - ptSphC ; 
      Vector3d vtInt2 = ptInt2 - ptSphC ; 
      vtInt1.Normalize() ;
      vtInt2.Normalize() ;
      int nInt = 2 ;
     // Scarto le soluzioni sulla calotta superiore
      if ( vtInt1 * vtSSphOrient > EPS_ZERO) {
         dU1 = dU2 ;
         -- nInt ;
      } 
      if ( vtInt2 * vtSSphOrient > EPS_ZERO)
         -- nInt ;
     // In base al numero di soluzioni valide  
     // aggiorno la tipologia di intersezione. 
      if ( nInt == 1)
         nIntType = S_ONE_INT_SEC ;
      else if ( nInt == 0)
         nIntType = S_NO_INTERS ;
   }
   return nIntType ; 
}

//----------------------------------------------------------------------------
int 
SegmentSemiSphere( const Point3d& ptPLine, const Vector3d& vtVLine, double dLen,
                   const Point3d& ptSphC, const Vector3d& vtSSphOrient, double dRad, 
                   double& dU1, double& dU2)
{
   int nIntType = SegmentSphere( ptPLine, vtVLine, dLen, ptSphC, dRad, dU1, dU2) ;
   if ( nIntType == S_ERROR_INT  ||  ( ! vtSSphOrient.IsNormalized()))
      return S_ERROR_INT ;
  // Un Punto in cui il segmento interferisce con la semi-sfera
   if ( nIntType == S_ONE_INT_SEC  ||  nIntType == S_ONE_INT_TAN) {
      Point3d ptInt = ptPLine + dU1 * vtVLine ;
      Vector3d vtInt = ptInt - ptSphC ; 
      vtInt.IsNormalized() ;
     // Se il punto di contatto è nella calotta superiore
     // Scarto la soluzione
      if ( vtInt * vtSSphOrient > EPS_ZERO) 
         nIntType = S_NO_INTERS ;
   }
  // Due punti in cui il segmento è secante
   else if ( nIntType == S_TWO_INT) {
      Point3d ptInt1 = ptPLine + dU1 * vtVLine ;
      Point3d ptInt2 = ptPLine + dU1 * vtVLine ;
      Vector3d vtInt1 = ptInt1 - ptSphC ; 
      Vector3d vtInt2 = ptInt2 - ptSphC ; 
      vtInt1.Normalize() ;
      vtInt2.Normalize() ;
      int nInt = 2 ;
     // Scarto le soluzioni sulla calotta superiore
      if ( vtInt1 * vtSSphOrient > EPS_ZERO) {
         dU1 = dU2 ;
         -- nInt ;
      } 
      if ( vtInt2 * vtSSphOrient > EPS_ZERO)
         -- nInt ;
     // In base al numero di soluzioni valide  
     // aggiorno la tipologia di intersezione. 
      if ( nInt == 1)
         nIntType = S_ONE_INT_SEC ;
      else if ( nInt == 0)
         nIntType = S_NO_INTERS ;
   }
   return nIntType ; 
}

//----------------------------------------------------------------------------
int 
LinCompSemiSphere( const Point3d& ptPLine, const Vector3d& vtVLine, double dLen, int nLinType, 
                   const Point3d& ptSphC, const Vector3d& vtSSphOrient, double dRad,              
                   double& dU1, double& dU2) 
{    
   int nIntType ;
  // Retta
   if ( nLinType == Line) 
      nIntType = LineSphere( ptPLine, vtVLine, ptSphC, dRad, dU1, dU2) ; 
  // Semi-retta
   else if ( nLinType == Ray) 
      nIntType = RaySphere( ptPLine, vtVLine, ptSphC, dRad, dU1, dU2) ;
  // Segmento
   else if ( nLinType == Segment) 
      nIntType = SegmentSphere( ptPLine, vtVLine, dLen, ptSphC, dRad, dU1, dU2) ;
  // Errore
   else 
      nIntType = S_ERROR_INT ;
  // Ancora errore
   if ( nIntType == S_ERROR_INT  ||  ( ! vtSSphOrient.IsNormalized()))
      return S_ERROR_INT ;
  // Un Punto in cui il segmento interferisce con la semi-sfera
   if ( nIntType == S_ONE_INT_SEC  ||  nIntType == S_ONE_INT_TAN) {
      Point3d ptInt = ptPLine + dU1 * vtVLine ;
      Vector3d vtInt = ptInt - ptSphC ; 
      vtInt.IsNormalized() ;
     // Se il punto di contatto è nella calotta superiore
     // Scarto la soluzione
      if ( vtInt * vtSSphOrient > EPS_ZERO) 
         nIntType = S_NO_INTERS ;
   }
  // Due punti in cui il segmento è secante
   else if ( nIntType == S_TWO_INT) {
      Point3d ptInt1 = ptPLine + dU1 * vtVLine ;
      Point3d ptInt2 = ptPLine + dU1 * vtVLine ;
      Vector3d vtInt1 = ptInt1 - ptSphC ; 
      Vector3d vtInt2 = ptInt2 - ptSphC ; 
      vtInt1.Normalize() ;
      vtInt2.Normalize() ;
      int nInt = 2 ;
     // Scarto le soluzioni sulla calotta superiore
      if ( vtInt1 * vtSSphOrient > EPS_ZERO) {
         dU1 = dU2 ;
         -- nInt ;
      } 
      if ( vtInt2 * vtSSphOrient > EPS_ZERO)
         -- nInt ;
     // In base al numero di soluzioni valide  
     // aggiorno la tipologia di intersezione. 
      if ( nInt == 1)
         nIntType = S_ONE_INT_SEC ;
      else if ( nInt == 0)
         nIntType = S_NO_INTERS ;
   }
   return nIntType ; 
}

//----------------------------------------------------------------------------
int
IntersLineInfiniteCylinder( const Point3d& ptPLine, const Vector3d& vtVLine,
                            const Point3d& ptPCyl, const Vector3d& vtVCyl, double dCylRad,
                            double& dU1, double& dU2)
{
  // Se il raggio non è significativamente maggiore di zero, vi è un errore.
   if ( dCylRad < EPS_SMALL)
      return CC_ERROR_INT ;
  // Si richiede che i vettori siano normalizzati
   if ( ! vtVCyl.IsNormalized()  ||  ! vtVLine.IsNormalized())
      return CC_ERROR_INT ;

  // La retta e l'asse del cilindro sono paralleli
   if ( AreSameOrOppositeVectorExact( vtVLine, vtVCyl)) {
      double dRadialDist = ( ( ptPLine - ptPCyl) - ( ( ptPLine - ptPCyl) * vtVCyl) * vtVCyl).Len() ;
     // La retta appartiene al cilindro 
      if ( abs( dRadialDist - dCylRad) < EPS_SMALL) {
         dU1 = - INFINITO ;
         dU2 = INFINITO ;
         return CC_INF_INT ;
      } 
     // non vi è intersezione
      else
         return CC_NO_INTERS ;
   }

  // Setto i coefficienti dell'equazione
   DBLVECTOR vdCoef( 3) ;
   vdCoef[0] = ( ptPLine - ptPCyl).SqLen() - ( ( ptPLine - ptPCyl) * vtVCyl) * ( ( ptPLine - ptPCyl) * vtVCyl) 
               - dCylRad * dCylRad ;
   vdCoef[1] = 2 * ( ( ptPLine - ptPCyl) * vtVLine - ( ( ptPLine - ptPCyl) * vtVCyl) * ( vtVCyl * vtVLine) ) ;
   vdCoef[2] = 1 - ( vtVCyl * vtVLine) * ( vtVCyl * vtVLine)  ;
  // Risolvo l'equazione
   DBLVECTOR vdRoots ;
   int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;
  // Studio delle soluzioni
   int nIntType = CC_NO_INTERS ;
  // C'è un punto di tangenza
   if ( nRoot == 1) {
      dU1 = vdRoots[0] ;
      nIntType = CC_ONE_INT_TAN ;
   }
  // L'equazione ha due soluzioni, bisogna valutare se sono ben distinte.
   else if ( nRoot == 2) {
      dU1 = vdRoots[0] ;
      dU2 = vdRoots[1] ; 
      if ( dU1 > dU2)
         swap( dU1, dU2) ;
     // Soluzioni sufficientemente lontane: ci sono due intersezioni
      if ( dU2 - dU1 > EPS_SMALL)                   
         nIntType = CC_TWO_INT ;    
     // Soluzioni non sufficientemente lontane: c'è un punto di tangenza.
      else 
         nIntType = CC_ONE_INT_TAN ;     
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
IntersRayInfiniteCylinder( const Point3d& ptPLine, const Vector3d& vtVLine,
                           const Point3d& ptPCyl, const Vector3d& vtVCyl, double dCylRad,
                           double& dU1, double& dU2)
{ 
   int nIntType = IntersLineInfiniteCylinder( ptPLine, vtVLine, ptPCyl, vtVCyl, dCylRad, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Se la retta associata appartiene al
  // cilindro, solo il parametro di inizio 
  // deve essere modificato
   if ( nIntType == CC_INF_INT) 
      dU1 = 0 ;
  // La retta associata è tangente: se la tangenza
  // ha parametro negativo, la semi-retta non interferisce
   else if ( nIntType == CC_ONE_INT_TAN) {
      if ( dU1 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ; 
   }
  // La retta associata è secante
   else if ( nIntType == CC_TWO_INT) {
     // Se dU1 è negativo, solo dU2 è accettabile
     // assegno dU2 a dU1
      if ( dU1 < - EPS_SMALL) {
         dU1 = dU2 ;
         nIntType = CC_ONE_INT_SEC ;
      }
     // Se anche dU2 è negativo la 
     // la semi-retta non interferisce
      if ( dU2 < - EPS_SMALL) 
         nIntType = CC_NO_INTERS ; 
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
IntersSegmentInfiniteCylinder( const Point3d& ptPLine, const Vector3d& vtVLine, double dLen,
                               const Point3d& ptPCyl, const Vector3d& vtVCyl, double dCylRad,
                               double& dU1, double& dU2)
{
   int nIntType = IntersLineInfiniteCylinder( ptPLine, vtVLine, ptPCyl, vtVCyl, dCylRad, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // La retta associata al segmento 
  // appartiene al cilindro.
   if ( nIntType == CC_INF_INT) {
      dU1 = 0 ;  
      dU2 = dLen ;
   }
  // La retta associata al segmento 
  // è tangente al cilindro.
   else if ( nIntType == CC_ONE_INT_TAN) {
     // Sel il parametro è fuori intervallo il
     // segmento non interferisce con il cilindro.
      if ( dU1 < - EPS_SMALL  ||  dU1 > dLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
  // La retta associata al segmento è 
  // secante al cilindro
   else if ( nIntType == CC_TWO_INT) {
     // Segmento non interagente
      if ( dU1 > dLen + EPS_SMALL) 
         nIntType = CC_NO_INTERS ; 
      else if ( dU1 > - EPS_SMALL) {
        // Segmento secante in dU1
         if ( dU2 > dLen + EPS_SMALL)    
            nIntType = CC_ONE_INT_SEC ;
      }
      else {
        // Segmento non interagente
         if ( dU2 > dLen + EPS_SMALL) {
            nIntType = CC_NO_INTERS ; 
         }
        // Segmento secante in dU2
         else if ( dU2 > - EPS_SMALL) {
            dU1 = dU2 ;
            nIntType = CC_ONE_INT_SEC ;
         }
        // Segmento non interagente
         else 
            nIntType = CC_NO_INTERS ;
      }
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
IntersLineSemiFiniteCylinder( const Point3d& ptPLine, const Vector3d& vtVLine,
                              const Point3d& ptPCyl, const Vector3d& vtVCyl, double dCylRad,
                              double& dU1, double& dU2)
{
   int nIntType = IntersLineInfiniteCylinder( ptPLine, vtVLine, ptPCyl, vtVCyl, dCylRad, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Caso parte della retta appartenente al cilindro
   if ( nIntType == CC_INF_INT) {
     // Retta e cilindro equiversi 
      if ( vtVCyl * vtVLine > 0)    
         dU1 = ( ptPLine - ptPCyl) * vtVCyl ;
     // Retta e cilindro controversi
      else {
         dU1 = - INFINITO ;
         dU2 = ( ptPLine - ptPCyl) * vtVCyl ;
      }      
   }   
  // Caso retta tangente
   else if ( nIntType == CC_ONE_INT_TAN) {
     // Se il punto di tangenza è fuori dalla parte 
     // ammissibile del cilindro, rigetto la soluzione
      Point3d ptInt = ptPLine + dU1 * vtVLine ;
      if ( ( ptInt - ptPCyl) * vtVCyl < 0)
         nIntType = CC_NO_INTERS ;
   }
  // Caso di due intersezioni
   else if ( nIntType == CC_TWO_INT) {
     // Rigetto le soluzioni fuori dalla  
     // regione ammissibile del cilindro.
      Point3d ptInt1 = ptPLine + dU1 * vtVLine ;
      Point3d ptInt2 = ptPLine + dU2 * vtVLine ;
      int nSolNum = 2 ;
      if ( ( ptInt1 - ptPCyl) * vtVCyl < 0) {
         -- nSolNum ;
         dU1 = dU2 ;
      }
      if ( ( ptInt2 - ptPCyl) * vtVCyl < 0) 
         -- nSolNum ; 
     // Dal numero di soluzioni rimaste 
     // aggiorno la tipologia di interferenza.
      if ( nSolNum == 1)
         nIntType = CC_ONE_INT_SEC ;
      else if ( nSolNum == 0) 
         nIntType = CC_NO_INTERS ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
RaySemiFiniteCylinder( const Point3d& ptPLine, const Vector3d& vtVLine, 
                       const Point3d& ptPCyl, const Vector3d& vtVCyl, double dCylRad, 
                       double& dU1, double& dU2)
{
   int nIntType = IntersLineSemiFiniteCylinder( ptPLine, vtVLine, ptPCyl, vtVCyl, dCylRad, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della semi-retta appartiene al cilindro  
   if ( nIntType == CC_INF_INT) {    
     // La semi-retta non interferisce
      if ( dU2 < - EPS_SMALL) 
         nIntType = CC_NO_INTERS ; 
     // La semi-retta incomincia nel cilindro
      else if ( dU1 < - EPS_SMALL)
            dU1 = 0 ;   
   } 
  // La retta associata è tangente  
  // o secante in un punto.
   else if ( nIntType == CC_ONE_INT_SEC  ||  
             nIntType ==  CC_ONE_INT_TAN) {
     // Se il parametro è negativo, 
     // scarto la soluzione.
      if ( dU1 < - EPS_SMALL) 
         nIntType = CC_NO_INTERS ;
   }
  // Retta associata secante
   else if ( nIntType == CC_TWO_INT) {
      int nSolNum = 2 ;
      if ( dU1 < - EPS_SMALL) {
         dU1 = dU2 ;
         -- nSolNum ;
      }
      if ( dU2 < - EPS_SMALL)
         -- nSolNum ;
     // Dal numero di soluzioni rimaste 
     // aggiorno la tipologia di interferenza.
      if ( nSolNum == 1)
         nIntType = CC_ONE_INT_SEC ;
      else if ( nSolNum == 0) 
         nIntType = CC_NO_INTERS ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
SegmentSemiFiniteCylinder( const Point3d& ptPLine, const Vector3d& vtVLine, double dLen, 
                           const Point3d& ptPCyl, const Vector3d& vtVCyl, double dCylRad,
                           double& dU1, double& dU2)
{
   int nIntType = IntersLineSemiFiniteCylinder( ptPLine, vtVLine, ptPCyl, vtVCyl, dCylRad, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta associata 
  // appartiene al cilindro.
   if ( nIntType == CC_INF_INT) {
     // Il segmento non interferisce 
      if ( dU1 > dLen + EPS_SMALL) 
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
        // Il segmento termina nel cilindro
         if ( dU2 > dLen + EPS_SMALL) 
            dU2 = dLen ;
      }
      else {
        // Il segmento contenuto totalmente
         if ( dU2 > dLen + EPS_SMALL) {
            dU1 = 0 ;
            dU2 = dLen ;
         }
        // Il segmento comincia nel cilindro
         else if ( dU2 > - EPS_SMALL) 
            dU1 = 0 ; 
        // Il segmeno non interferisce 
         else 
            nIntType = CC_NO_INTERS ;
      }
   }
  // La retta associata è tangente  
  // o secante in un punto.
   else if ( nIntType == CC_ONE_INT_SEC  ||  
             nIntType == CC_ONE_INT_TAN) {
     // Il sgmento non interferisce
      if ( dU1 < 0  ||  dU1 > dLen)
         nIntType = CC_NO_INTERS ;
   }
  // La retta associata è secante
   else if ( nIntType == CC_TWO_INT) {
     // Il segmento non interferisce
      if ( dU1 > dLen + EPS_SMALL) 
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
        // Segmento secante in dU1
         if ( dU2 > dLen + EPS_SMALL)    
            nIntType = CC_ONE_INT_SEC ;
      }
      else {
        // Il segmento non interagisce
         if ( dU2 > dLen + EPS_SMALL) {
            nIntType = CC_NO_INTERS ;
         }
        // Il segmento è secante in dU2, 
        // assegno dU2 a dU1.
         else if ( dU2 > - EPS_SMALL) {
            dU1 = dU2 ;
            nIntType = CC_ONE_INT_SEC ;
         }
        // Il segmento non interferisce
         else 
            nIntType = CC_NO_INTERS ;
      }
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
IntersLineCylinder( const Point3d& ptPLine, const Vector3d& vtVLine,
                    const Point3d& ptPCyl, const Vector3d& vtVCyl, double dCylRad, double dCylHeigth,
                    double& dU1, double& dU2) 
{
   int nIntType = IntersLineInfiniteCylinder( ptPLine, vtVLine, ptPCyl, vtVCyl, dCylRad, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta appartenente al cilindro
   if ( nIntType == CC_INF_INT) {
      dU1 = ( ptPCyl - ptPLine) * vtVLine ;
     // Retta e cilindro equiversi
      if ( vtVCyl * vtVLine > 0)      
         dU2 = dU1 + dCylHeigth ;           
     // Retta e cilindro controversi
      else        
         dU2 = dU1 - dCylHeigth ;   
   }   
  // Caso di un punto di tangenza
   else if ( nIntType == CC_ONE_INT_TAN) {
     // Se il punto di tangenza è fuori dalla parte 
     // ammissibile del cilindro, rigetto la soluzione
      Point3d ptInt = ptPLine + dU1 * vtVLine ;
      if ( ( ptInt - ptPCyl) * vtVCyl < 0  ||
           ( ptInt - ptPCyl) * vtVCyl > dCylHeigth)
         nIntType = CC_NO_INTERS ;
   }
  // Caso di due intersezioni
   else if ( nIntType == CC_TWO_INT) {
     // Elimino le soluzioni fuori dal cilindro finito
      Point3d ptInt1 = ptPLine + dU1 * vtVLine ;
      Point3d ptInt2 = ptPLine + dU2 * vtVLine ;
      int nSolNum = 2 ;
      if ( ( ptInt1 - ptPCyl) * vtVCyl < 0  ||
           ( ptInt1 - ptPCyl) * vtVCyl > dCylHeigth) {
         -- nSolNum ;
         dU1 = dU2 ;
      }
      if ( ( ptInt2 - ptPCyl) * vtVCyl < 0  ||
           ( ptInt2 - ptPCyl) * vtVCyl > dCylHeigth) 
         -- nSolNum ;
     // Aggiorno il tipo di intersezione
      if ( nSolNum == 1)
         nIntType = CC_ONE_INT_SEC ;
      else if ( nSolNum == 0) 
         nIntType = CC_NO_INTERS ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
RayCylinder( const Point3d& ptPLine, const Vector3d& vtVLine, 
             const Point3d& ptPCyl, const Vector3d& vtVCyl, double dCylRad, double dCylHeigth, 
             double& dU1, double& dU2)
{
   int nIntType = IntersLineCylinder( ptPLine, vtVLine, ptPCyl, vtVCyl, dCylRad, dCylHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta associata 
  // appartiene al cilindro
   if ( nIntType == CC_INF_INT) {
     // Semi-retta non interferisce
      if ( dU2 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ;
     // La semi-retta comincia nel cilindro
      else if ( dU1 < - EPS_SMALL)
         dU1 = 0 ;
   }
  // Retta associata tangente
   else if ( nIntType == CC_ONE_INT_TAN) {
      if ( dU1 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
   else  if ( nIntType ==  CC_TWO_INT) {
     // Scarto soluzioni con parametro negativo
      int nSolNum = 2 ;
      if ( dU1 < - EPS_SMALL) {
         dU1 = dU2 ;
         -- nSolNum ;
      }
      if ( dU2 < - EPS_SMALL) 
         -- nSolNum ;
     // Dal numero di soluzioni rimaste 
     // aggiorno la tipologia di interferenza.
      if ( nSolNum == 1)
         nIntType = CC_ONE_INT_SEC ;
      else if ( nSolNum == 0) 
         nIntType = CC_NO_INTERS ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
IntersSegmentCylinder( const Point3d& ptPLine, const Vector3d& vtVLine, double dLen,
                       const Point3d& ptPCyl,  const Vector3d& vtVCyl, double dCylRad, double dCylHeigth, 
                       double& dU1, double& dU2)
{
   int nIntType = IntersLineCylinder( ptPLine, vtVLine, ptPCyl, vtVCyl, dCylRad, dCylHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta associata appartiene al cilindro
   if ( nIntType == CC_INF_INT) {
     // Segmento non interferisce
      if ( dU1 > dLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
         if ( dU2 > dLen + EPS_SMALL)
            dU2 = dLen ;        
      }
      else {
        // Segmento non interferisce
         if ( dU2 < - EPS_SMALL)
            nIntType = CC_NO_INTERS ;
        // Il segmento comincia nel cilindro
         else {
            dU1 = 0 ;
           // Il segmento termina nel cilindro
            if ( dU2 > dLen + EPS_SMALL)
            dU2 = dLen ;  
         }
      }
   }
  // La retta associata è tangente al cilindro
   else if ( nIntType == CC_ONE_INT_TAN) {
     // Segmento non interferisce
      if ( dU1 < - EPS_SMALL  ||  dU1 > dLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
  // Retta associata è secante con due punti di intersezione
   else if ( nIntType == CC_TWO_INT) {
     // Il segmento non interferisce
      if ( dU1 > dLen + EPS_SMALL) 
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
        // Segmento secante in dU1
         if ( dU2 > dLen + EPS_SMALL)    
            nIntType = CC_ONE_INT_SEC ;
      }
      else {
        // Il segmento non interagisce
         if ( dU2 > dLen + EPS_SMALL) {
            nIntType = CC_NO_INTERS ;
         }
        // Il segmento è secante in dU2, 
        // assegno dU2 a dU1.
         else if ( dU2 > - EPS_SMALL) {
            dU1 = dU2 ;
            nIntType = CC_ONE_INT_SEC ;
         }
        // Il segmento non interferisce
         else 
            nIntType = CC_NO_INTERS ;
      }
   }
  // Retta associata è secante con un punto di intersezione
   else if ( nIntType == CC_ONE_INT_SEC) {
      if ( dU1 < -EPS_SMALL  ||  dU1 > dLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
LineInfiniteCone( const Point3d& ptPLine, const Vector3d& vtDLine, 
                  const Point3d& ptVCone, const Vector3d& vtDCone, double dConeRad, double dConeHeigth, 
                  double& dU1, double& dU2)
{
  // Raggio e altezza del cono devono essere maggiori di zero
   if ( dConeRad < EPS_SMALL  ||  dConeHeigth < EPS_SMALL)
      return CC_ERROR_INT ;
  // Si richiede che i vettori siano normalizzati
   if ( ! vtDCone.IsNormalized()  ||  ! vtDLine.IsNormalized())
      return CC_ERROR_INT ;

  // Tangente dell'angolo semi-apertura del cono
   double dTanTheta = dConeRad / dConeHeigth ;

  // Lunghezze delle componenti parallela e ortogonale del 
  // vettore della retta rispetto all'asse del cilindro.
   double dDLong = vtDLine * vtDCone ;
   double dDOrt = (vtDLine - dDLong * vtDCone).Len() ;
  // Flag parallelismo fra la retta e una retta generatrice del cono.
   bool bSameAngle = abs( dDLong) > EPS_ZERO  &&  abs( dDOrt / abs( dDLong) - dTanTheta) < EPS_ZERO ;
   if ( bSameAngle) {
     // Valuto se parte della retta appartiene al cono.
      Vector3d vtOO =  ptVCone - ptPLine ;
      Vector3d vtOOPerp = vtOO - ( vtOO * vtDLine) * vtDLine ;
     // Caso in cui la retta appartiene al cono
      if ( vtOOPerp.SqLen() < EPS_SMALL * EPS_SMALL) {
        // Calcolo il parametro a cui si trova il 
        // vertice del cono lungo la retta.
         if ( vtDLine * vtDCone > 0) {
            dU1 = vtOO * vtDLine ;
            dU2 = INFINITO ;
         }
         else {
            dU1 = - INFINITO ;
            dU2 = vtOO * vtDLine ;
         }
         return CC_INF_INT ;
      }
     // Vertice del cono non sta sulla retta
      else {
         double dC0 = ( ptPLine - ptVCone) * ( ptPLine - ptVCone) -
                      ( 1 + dTanTheta * dTanTheta) * ( ( ptPLine - ptVCone) * vtDCone) * ( ( ptPLine - ptVCone) * vtDCone) ;
         double dC1 = 2 * ( vtDLine * ( ptPLine - ptVCone) - 
                      ( 1 + dTanTheta * dTanTheta) * ( vtDCone * vtDLine) * ( ptPLine - ptVCone) * vtDCone) ;
         if ( abs( dC1) > EPS_ZERO) {
            dU1 = - dC0 / dC1 ;
            Point3d ptInt = ptPLine + dU1 * vtDLine ;
           // Se la soluzione è sulla falda negativa deve essere eliminata.
            if ( ( ptInt - ptVCone) * vtDCone < 0) 
               return CC_NO_INTERS ;
            return CC_ONE_INT_SEC ;
         }
         else if ( abs( dC0) > EPS_SMALL)
            return CC_NO_INTERS ;         
      }
   }

  // Setto i coefficienti dell'equazione
   DBLVECTOR vdCoef( 3) ;
   double dDiffDiff = ( ptPLine - ptVCone) * ( ptPLine - ptVCone) ;
   double dDiffDC = ( ptPLine - ptVCone) * vtDCone ;
   double dDiffDL = ( ptPLine - ptVCone) * vtDLine ;
   double dDCDL = vtDCone * vtDLine ;
   double dInvSqCosTheta = 1 + dTanTheta * dTanTheta ;
   vdCoef[0] = dDiffDiff - dInvSqCosTheta * dDiffDC * dDiffDC ; 
   vdCoef[1] = 2 * ( dDiffDL - dInvSqCosTheta * dDCDL * dDiffDC) ;
   vdCoef[2] = 1 - dInvSqCosTheta * dDCDL * dDCDL ;  

  // Risolvo l'equazione
   DBLVECTOR vdRoots ;
   int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;

  // Studio le soluzioni
   int nIntType = CC_NO_INTERS ;
   if ( nRoot == 1) {
      dU1 = vdRoots[0] ;
      Point3d ptInt = ptPLine + dU1 * vtDLine ;
     // Soluzione sul cono negativo
      if ( ( ptInt - ptVCone) * vtDCone < 0)
         nIntType = CC_NO_INTERS ; 
     // Soluzione sul vertice
      else if ( SqDist( ptInt, ptVCone) < EPS_SMALL * EPS_SMALL)
         nIntType = CC_ON_VERT ;
     // Retta parallela e secante in un punto
      else if ( bSameAngle) 
         nIntType = CC_ONE_INT_SEC ;
     // Retta tangente al cono
      else 
         nIntType = CC_ONE_INT_TAN ;
   }
   else if ( nRoot == 2) {
      dU1 = vdRoots[0] ;
      dU2 = vdRoots[1] ;
      if ( dU1 > dU2)
         swap( dU1, dU2) ;      
     // Se le soluzioni non sono sufficientemente lontane,
     // la soluzione è di tangenza
      if ( dU2 - dU1 < EPS_SMALL * EPS_SMALL) {
         Point3d ptInt1 = ptPLine + dU1 * vtDLine ;
        // Soluzione sul cono negativo
         if ( ( ptInt1 - ptVCone) * vtDCone < 0)
            nIntType = CC_NO_INTERS ; 
        // Soluzione sul vertice
         else if ( SqDist( ptInt1, ptVCone) < EPS_SMALL * EPS_SMALL)
            nIntType = CC_ON_VERT ;
        // Retta parallela e secante in un punto
         else if ( bSameAngle) 
            nIntType = CC_ONE_INT_SEC ;
        // Retta tangente al cono
         else 
            nIntType = CC_ONE_INT_TAN ;
      }
     // Soluzioni ben distinte
      else {
         Point3d ptInt1 = ptPLine + dU1 * vtDLine ;
         Point3d ptInt2 = ptPLine + dU2 * vtDLine ;
        // Elimino le soluzioni sul cono negativo
         int nNumSol = 2 ;
         if ( ( ptInt2 - ptVCone) * vtDCone < 0) 
            -- nNumSol ;
         if ( ( ptInt1 - ptVCone) * vtDCone < 0) {
            -- nNumSol ;          
            dU1 = dU2 ;
         }
        // Dal numero di soluzioni rimaste 
        // aggiorno la tipologia di interferenza.
         if ( nNumSol == 0)
            nIntType = CC_NO_INTERS ;
         else if ( nNumSol == 1)
            nIntType = CC_ONE_INT_SEC ;
         else if ( nNumSol == 2)
            nIntType = CC_TWO_INT ;
      }   
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int
RayInfiniteCone( const Point3d& ptPLine, const Vector3d& vtDLine,
                 const Point3d& ptVCone, const Vector3d& vtDCone, double dConeRad, double dConeHeigth, 
                 double& dU1, double& dU2)
{
   int nIntType = LineInfiniteCone( ptPLine, vtDLine, ptVCone, vtDCone, dConeRad, dConeHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta associata  
  // appartiene al cono.
   if ( nIntType == CC_INF_INT) {
      if ( vtDLine * vtDCone > 0) {
        // Semi-retta comincia nel cono
         if ( dU1 < 0)
            dU1 = 0 ;
      }
      else {
        // Semiretta non interferisce
        // con il cono.
         if ( dU2 < 0)
            nIntType = CC_NO_INTERS ; 
        // Semi-retta comincia nel cono
         else 
            dU1 = 0 ;
      }
   }
  // Retta associata interferisce in un punto
   else if ( nIntType == CC_ONE_INT_TAN  ||
             nIntType == CC_ONE_INT_SEC  ||
             nIntType == CC_ON_VERT) {
     // Semi-retta non interferisce
      if ( dU1 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ; 
   }
  // Retta associata interferisce in due punti.
   else if ( nIntType == CC_TWO_INT) {
     // Semi-retta non interferisce
      if ( dU2 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ; 
     // Semi-retta secante in un punto
      else if ( dU1 < - EPS_SMALL)
         nIntType = CC_ONE_INT_SEC ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int 
SegmentInfiniteCone( const Point3d& ptPLine, const Vector3d& vtDLine, double dSgLen,
                     const Point3d& ptVCone, const Vector3d& vtDCone, double dConeRad, double dConeHeigth,
                     double& dU1, double& dU2)
{
   int nIntType = LineInfiniteCone( ptPLine, vtDLine, ptVCone, vtDCone, dConeRad, dConeHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta associata 
  // appartiene al cono 
   if ( nIntType == CC_INF_INT) {
     // Il segmento non interferisce 
      if ( dU1 > dSgLen + EPS_SMALL) 
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
        // Il segmento termina nel cono
         if ( dU2 > dSgLen + EPS_SMALL) 
            dU2 = dSgLen ;
      }
      else {
        // Il segmento contenuto totalmente
         if ( dU2 > dSgLen + EPS_SMALL) {
            dU1 = 0 ;
            dU2 = dSgLen ;
         }
        // Il segmento comincia nel cono
         else if ( dU2 > - EPS_SMALL) 
            dU1 = 0 ; 
        // Il segmeno non interferisce 
         else 
            nIntType = CC_NO_INTERS ;
      }
   }
  // Retta associata interferisce in un punto
   else if ( nIntType == CC_ONE_INT_TAN  ||
             nIntType == CC_ONE_INT_SEC  ||
             nIntType == CC_ON_VERT) {
     // Segmento non interagente
      if ( dU1 < - EPS_SMALL  ||  dU1 > dSgLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
  // Retta associata interferisce in due punti.
   else if ( nIntType == CC_TWO_INT) {
     // Il segmento non interferisce
      if ( dU1 > dSgLen + EPS_SMALL) 
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
        // Segmento secante in dU1
         if ( dU2 > dSgLen + EPS_SMALL)    
            nIntType = CC_ONE_INT_SEC ;
      }
      else {
        // Il segmento non interagisce
         if ( dU2 > dSgLen + EPS_SMALL) {
            nIntType = CC_NO_INTERS ;
         }
        // Il segmento è secante in dU2, 
        // assegno dU2 a dU1.
         else if ( dU2 > - EPS_SMALL) {
            dU1 = dU2 ;
            nIntType = CC_ONE_INT_SEC ;
         }
        // Il segmento non interferisce
         else 
            nIntType = CC_NO_INTERS ;
      }     
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int 
LineCone( const Point3d& ptPLine, const Vector3d& vtDLine,
          const Point3d& ptVCone, const Vector3d& vtDCone, double dConeRad, double dConeHeigth, 
          double& dU1, double& dU2)
{
   int nIntType = LineInfiniteCone( ptPLine, vtDLine, ptVCone, vtDCone, dConeRad, dConeHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;
  // Tangente dell'angolo semi-apertura del cono  
   double dTanTheta = dConeRad / dConeHeigth ; 
  // Parte della retta associata giace sul cono
   if ( nIntType == CC_INF_INT) {
      double dApothem = dConeHeigth * sqrt( 1 + dTanTheta * dTanTheta) ;
     // Retta e asse del cono equiversi
      if ( vtDLine * vtDCone > 0) 
         dU2 = dU1 + dApothem ;        
     // Retta e asse del cono controversi
      else         
         dU1 = dU2 - dApothem ;     
   }
  // Retta e il cono infinito interagiscono in un punto
   else if ( nIntType == CC_ONE_INT_TAN  ||
             nIntType == CC_ONE_INT_SEC  ||
             nIntType == CC_ON_VERT) {
     // Se la soluzione è fuori dal cono la rigetto
      Point3d ptInt = ptPLine + dU1 * vtDLine ;    
      if ( ( ptInt - ptVCone) * vtDCone > dConeHeigth)
         nIntType = CC_NO_INTERS ; ;
    }
  // Retta e il cono infinito interagiscono in due punti
    else if ( nIntType == CC_TWO_INT) {
      // Rigetto le soluzioni fuori dal cono
       Point3d ptInt1 = ptPLine + dU1 * vtDLine ;
       Point3d ptInt2 = ptPLine + dU2 * vtDLine ; 
       int nNumSol = 2 ;
       if ( ( ptInt2 - ptVCone) * vtDCone > dConeHeigth) 
          -- nNumSol ;
       if ( ( ptInt1 - ptVCone) * vtDCone > dConeHeigth) {
          -- nNumSol ;         
       dU1 = dU2 ;
       }
     // Dal numero di soluzioni rimaste 
     // aggiorno la tipologia di interferenza.
      if ( nNumSol == 0)
         nIntType = CC_NO_INTERS ;
      else if ( nNumSol == 1)
         nIntType = CC_ONE_INT_SEC ;
      else if ( nNumSol == 2)
         nIntType = CC_TWO_INT ;
    }
    return nIntType ;
}

//----------------------------------------------------------------------------
int 
RayCone( const Point3d& ptPLine, const Vector3d& vtDLine,
         const Point3d& ptVCone, const Vector3d& vtDCone, double dConeRad, double dConeHeigth,
         double& dU1, double& dU2)
{
   int nIntType = LineCone( ptPLine, vtDLine, ptVCone, vtDCone, dConeRad, dConeHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;
  // Parte della retta associata 
  // giace sul cono
   if ( nIntType == CC_INF_INT) { 
      double dApothem = sqrt( dConeRad * dConeRad + dConeHeigth * dConeHeigth) ; 
     // Semi-retta esterna
      if ( dU2 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ;
     // Una soluzione tangente
      else if ( dU2 < EPS_SMALL) {
          dU1 = 0 ;
          nIntType = CC_ONE_INT_TAN ;
      }   
     // Semi-retta comincia nel cono
      else if ( dU2 < dApothem)
         dU1 = 0 ;
   } 
  // Retta associata interagisce in un punto
   else if ( nIntType == CC_ONE_INT_TAN  ||
             nIntType == CC_ONE_INT_SEC  ||
             nIntType == CC_ON_VERT) {
     // Semi-retta esterna
      if ( dU1 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
  // Retta associata interagisce in due punti
   else if ( nIntType == CC_TWO_INT) {
     // Se dU1 è negativo, solo dU2 è accettabile
     // assegno dU2 a dU1
      if ( dU1 < - EPS_SMALL) {
         dU1 = dU2 ;
         nIntType = CC_ONE_INT_SEC ;
      }
     // Se anche dU2 è negativo la 
     // la semi-retta non interferisce
      if ( dU2 < - EPS_SMALL) 
         nIntType = CC_NO_INTERS ; 
   }
   return nIntType ;
}
                                                          
//----------------------------------------------------------------------------
int
SegmentCone( const Point3d& ptPLine, const Vector3d& vtDLine, double dSgLen,
             const Point3d& ptVCone, const Vector3d& vtDCone, double dConeRad, double dConeHeigth, 
             double& dU1, double& dU2) 
{
  // Verifico se il segmento interseca un insieme di box impilati che includono il cono
   Frame3d frCone ; frCone.Set( ptVCone, vtDCone) ;
   Point3d ptLineL = ptPLine ; ptLineL.ToLoc( frCone) ;
   Vector3d vtLineL = vtDLine ; vtLineL.ToLoc( frCone) ;
   const int N_SLABS = 4 ;
   double dH = dConeHeigth / N_SLABS ;
   double dR = dConeRad / N_SLABS ;
   int nMiss = 0 ;
   for ( int i = 1 ; i <= N_SLABS ; ++ i) {
      double dPar1, dPar2 ;
      ptLineL.z -= ( i - 1) * dH ; 
      double dCurrR = i * dR ;
      if ( ! IntersLineBox( ptLineL, vtLineL, Point3d( -dCurrR, -dCurrR, 0), Point3d( dCurrR, dCurrR, dH), dPar1, dPar2)  ||
           dPar2 < -EPS_SMALL  ||  dPar1 > dSgLen + EPS_SMALL)
         ++ nMiss ;
   }
   if ( nMiss == N_SLABS)
      return CC_NO_INTERS ;

  // Calcolo l'intersezione con la linea illimitata
   int nIntType = LineCone( ptPLine, vtDLine, ptVCone, vtDCone, dConeRad, dConeHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

   if ( nIntType == CC_INF_INT) {    
      // Segmento non interferisce
      if ( dU1 > dSgLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
         if ( dU2 > dSgLen + EPS_SMALL)
            dU2 = dSgLen ;        
      }
      else {
        // Segmento non interferisce
         if ( dU2 < - EPS_SMALL)
            nIntType = CC_NO_INTERS ;
        // Il segmento comincia nel cono
         else {
            dU1 = 0 ;
           // Il segmento termina nel cono
            if ( dU2 > dSgLen + EPS_SMALL)
            dU2 = dSgLen ;  
         }
      }
   }
  // Retta associata interagisce in un punto
   else if ( nIntType == CC_ONE_INT_TAN  ||
             nIntType == CC_ONE_INT_SEC  ||
             nIntType == CC_ON_VERT) {
     // Segmento non interferisce
      if ( dU1 < - EPS_SMALL  ||  dU1 > dSgLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
  // Retta associata interagisce in due punti
   else if ( nIntType == CC_TWO_INT) {
      // Il segmento non interferisce
      if ( dU1 > dSgLen + EPS_SMALL) 
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
        // Segmento secante in dU1
         if ( dU2 > dSgLen + EPS_SMALL)    
            nIntType = CC_ONE_INT_SEC ;
      }
      else {
        // Il segmento non interagisce
         if ( dU2 > dSgLen + EPS_SMALL) {
            nIntType = CC_NO_INTERS ;
         }
        // Il segmento è secante in dU2, 
        // assegno dU2 a dU1.
         else if ( dU2 > - EPS_SMALL) {
            dU1 = dU2 ;
            nIntType = CC_ONE_INT_SEC ;
         }
        // Il segmento non interferisce
         else 
            nIntType = CC_NO_INTERS ;
      }     
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int 
LineConeFrustum( const Point3d& ptPLine, const Vector3d& vtDLine,
                 const Point3d& ptMinBase, const Vector3d& vtDCone, double dMinRad, double dMaxRad, double dHeigth, 
                 double& dU1, double& dU2)
{
  // Controlli sull'ammissibilità del tronco di cono
   if ( ! ( dMinRad > EPS_SMALL  &&  dMaxRad > EPS_SMALL  &&  dHeigth > EPS_SMALL))
      return CC_ERROR_INT ;

   double dDeltaH = dMinRad * dHeigth / ( dMaxRad - dMinRad) ; 
   Point3d ptVCone = ptMinBase - dDeltaH * vtDCone ;

   int nIntType = LineInfiniteCone( ptPLine, vtDLine, ptVCone, vtDCone, dMinRad, dDeltaH, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta associata giace sul cono
   if ( nIntType == CC_INF_INT) {
    
      double dMinApo = sqrt( dDeltaH * dDeltaH + dMinRad * dMinRad) ;
      double dMaxApo = sqrt( ( dDeltaH + dHeigth) * ( dDeltaH + dHeigth) + dMaxRad * dMaxRad) ;
     // Retta e vettore cono equiversi
      if ( vtDLine * vtDCone > 0) {
         double dSignedLenOO = ( ptPLine - ptVCone) * vtDLine ;   
         dU1 = dMinApo - dSignedLenOO ;
         dU2 = dMaxApo - dSignedLenOO ;       
      }
     // Retta e vettore cono controversi
      else {
         double dSignedLenOO = ( ptPLine - ptVCone) * vtDLine ;   
         dU1 = - dMaxApo - dSignedLenOO ;
         dU2 = - dMinApo - dSignedLenOO ; 
      }
   }
  // Retta associata interagisce in un punto
   else if ( nIntType == CC_ONE_INT_TAN  ||
             nIntType == CC_ONE_INT_SEC  ||
             nIntType == CC_ON_VERT) {

      Point3d ptInt = ptPLine + dU1 * vtDLine ;

      if ( ( ptInt - ptVCone) * vtDCone < dDeltaH - EPS_SMALL ||
           ( ptInt - ptVCone) * vtDCone > dDeltaH + dHeigth + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
  // Retta associata interagisce in due punti
   else if ( nIntType == CC_TWO_INT) {
     // Elimino le soluzioni da fuori dal tronco
      Point3d ptInt1 = ptPLine + dU1 * vtDLine ;
      Point3d ptInt2 = ptPLine + dU2 * vtDLine ;
      int nNumSol = 2 ;
      if ( ( ptInt2 - ptVCone) * vtDCone < dDeltaH  ||
           ( ptInt2 - ptVCone) * vtDCone > dDeltaH + dHeigth)
         -- nNumSol ;
      if ( ( ptInt1 - ptVCone) * vtDCone < dDeltaH  ||
           ( ptInt1 - ptVCone) * vtDCone > dDeltaH + dHeigth) {
         -- nNumSol ;     
         dU1 = dU2 ;
      }
     // Dal numero di soluzioni rimaste 
     // aggiorno la tipologia di interferenza.
      if ( nNumSol == 0)
         nIntType = CC_NO_INTERS ;
      else if ( nNumSol == 1)
         nIntType = CC_ONE_INT_SEC ;
      else if ( nNumSol == 2)
         nIntType = CC_TWO_INT ;
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int 
RayConeFrustum( const Point3d& ptPLine, const Vector3d& vtDLine, 
                const Point3d& ptMinBase, const Vector3d& vtDCone, double dMinRad, double dMaxRad, double dHeigth,
                double& dU1, double& dU2)
{
   int nIntType = LineConeFrustum( ptPLine, vtDLine, ptMinBase, vtDCone, dMinRad, dMaxRad, dHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta associata giace sul tronco di cono
   if ( nIntType == CC_INF_INT) {
      double dDeltaH = dMinRad * dHeigth / ( dMaxRad - dMinRad) ; 
      double dMinApo = sqrt( dDeltaH * dDeltaH + dMinRad * dMinRad) ;
      double dMaxApo = sqrt( ( dDeltaH + dHeigth) * ( dDeltaH + dHeigth) + dMaxRad * dMaxRad) ;
     // Semi-retta non interagisce 
      if ( dU2 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ;
     // Semi-retta parte sul bordo del tronco 
      else if ( dU2 < EPS_SMALL) {
         dU1 = 0 ;
         nIntType = CC_ONE_INT_TAN ;
      }
     // Semi-retta comincia nel tronco
      else if ( dU2 < dMaxApo - dMinApo)
         dU1 = 0 ;
   }
  // Retta associata interagisce in un punto
   else if ( nIntType == CC_ONE_INT_TAN  ||
             nIntType == CC_ONE_INT_SEC  ||
             nIntType == CC_ON_VERT) {
     // Semi-retta non interagisce
      if ( dU1 < - EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
  // Retta associata interagisce in un punto
   else if ( nIntType == CC_TWO_INT) {
     // Se dU1 è negativo, solo dU2 è accettabile
     // assegno dU2 a dU1
      if ( dU1 < - EPS_SMALL) {
         dU1 = dU2 ;
         nIntType = CC_ONE_INT_SEC ;
      }
     // Se anche dU2 è negativo la 
     // la semi-retta non interferisce
      if ( dU2 < - EPS_SMALL) 
         nIntType = CC_NO_INTERS ;      
   }
    return nIntType ;
}

//----------------------------------------------------------------------------
int
SegmentConeFrustum( const Point3d& ptPLine, const Vector3d& vtDLine, double dSgLen, 
                    const Point3d& ptMinBase, const Vector3d& vtDCone, double dMinRad, double dMaxRad, double dHeigth, 
                    double& dU1, double& dU2)
{
  // Verifico se il segmento interseca un insieme di box impilati che includono il tronco cono
   Frame3d frCone ; frCone.Set( ptMinBase, vtDCone) ;
   Point3d ptLineL = ptPLine ; ptLineL.ToLoc( frCone) ;
   Vector3d vtLineL = vtDLine ; vtLineL.ToLoc( frCone) ;
   const int N_SLABS = 4 ;
   double dH = dHeigth / N_SLABS ;
   double dR = ( dMaxRad - dMinRad) / N_SLABS ;
   int nMiss = 0 ;
   for ( int i = 1 ; i <= N_SLABS ; ++ i) {
      double dPar1, dPar2 ;
      ptLineL.z -= ( i - 1) * dH ; 
      double dCurrR = dMinRad + i * dR ;
      if ( ! IntersLineBox( ptLineL, vtLineL, Point3d( -dCurrR, -dCurrR, 0), Point3d( dCurrR, dCurrR, dH), dPar1, dPar2)  ||
           dPar2 < -EPS_SMALL  ||  dPar1 > dSgLen + EPS_SMALL)
         ++ nMiss ;
   }
   if ( nMiss == N_SLABS)
      return CC_NO_INTERS ;

  // Calcolo l'intersezione con la linea illimitata
   int nIntType = LineConeFrustum( ptPLine, vtDLine, ptMinBase, vtDCone, dMinRad, dMaxRad, dHeigth, dU1, dU2) ;
   if ( nIntType == CC_ERROR_INT)
      return nIntType ;

  // Parte della retta associata giace sul tronco
   if ( nIntType == CC_INF_INT) {
      // Segmento non interferisce
      if ( dU1 > dSgLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
         if ( dU2 > dSgLen + EPS_SMALL)
            dU2 = dSgLen ;        
      }
      else {
        // Segmento non interferisce
         if ( dU2 < - EPS_SMALL)
            nIntType = CC_NO_INTERS ;
        // Il segmento comincia nel cilindro
         else {
            dU1 = 0 ;
           // Il segmento termina nel cilindro
            if ( dU2 > dSgLen + EPS_SMALL)
            dU2 = dSgLen ;  
         }
      }
   }
  // Retta associata interagisce in un punto
   else if (  nIntType == CC_ONE_INT_TAN  ||
              nIntType == CC_ONE_INT_SEC  ||
              nIntType == CC_ON_VERT) {
      // Segmento non interferisce
      if ( dU1 < - EPS_SMALL  ||  dU1 > dSgLen + EPS_SMALL)
         nIntType = CC_NO_INTERS ;
   }
  // Retta associata interagisce in due punti
   else if ( nIntType == CC_TWO_INT) {
     // Il segmento non interferisce
      if ( dU1 > dSgLen + EPS_SMALL) 
         nIntType = CC_NO_INTERS ;
      else if ( dU1 > - EPS_SMALL) {
        // Segmento secante in dU1
         if ( dU2 > dSgLen + EPS_SMALL)    
            nIntType = CC_ONE_INT_SEC ;
      }
      else {
        // Il segmento non interagisce
         if ( dU2 > dSgLen + EPS_SMALL) {
            nIntType = CC_NO_INTERS ;
         }
        // Il segmento è secante in dU2, 
        // assegno dU2 a dU1.
         else if ( dU2 > - EPS_SMALL) {
            dU1 = dU2 ;
            nIntType = CC_ONE_INT_SEC ;
         }
        // Il segmento non interferisce
         else 
            nIntType = CC_NO_INTERS ;
      }     
   }
   return nIntType ;
}

//----------------------------------------------------------------------------
int 
LineTorus( const Point3d& ptLine, const Vector3d& vtLine,
           const Point3d& ptOTorus, const Vector3d& vtAxTorus, double dMinRad, double dMaxRad,
           BOOLVECTOR& vbType, DBLVECTOR& vdPar)
{
  // I Raggi devono essere positiv 
   if ( dMinRad < EPS_SMALL  ||  dMaxRad < EPS_SMALL)
      return T_ERROR ;
  // Si richiede che i vettori siano normalizzati
   if ( ! vtAxTorus.IsNormalized()  ||  ! vtLine.IsNormalized())
      return T_ERROR ;

  // Sistema di riferimento del toro
   Frame3d frTorusFrame ;
   frTorusFrame.Set( ptOTorus, vtAxTorus) ;
  // Porto la retta nel sistema di riferimento del toro
   Point3d ptLn = ptLine ;
   ptLn.ToLoc( frTorusFrame) ;
   Vector3d vtLn = vtLine ;
   vtLn.ToLoc( frTorusFrame) ;
  
  // Setto i coefficienti dell'equazione
   DBLVECTOR vdCoef( 5) ;
   double dSqDistPO = ptLn.x * ptLn.x + ptLn.y * ptLn.y + ptLn.z * ptLn.z ;
   double dPtVt = ptLn.x * vtLn.x + ptLn.y * vtLn.y + ptLn.z * vtLn.z ;
   double dSqMaxR = dMaxRad * dMaxRad ;
   double dDeltaSqR = dSqMaxR - dMinRad * dMinRad ;
   double dB = 2 * dDeltaSqR - 4 * dSqMaxR ; 
   vdCoef[0] = dSqDistPO * dSqDistPO + dB * ( dSqDistPO - ptLn.z * ptLn.z) +
               2 * dDeltaSqR * ptLn.z * ptLn.z + dDeltaSqR * dDeltaSqR ;
   vdCoef[1] = 4 * dPtVt * dSqDistPO + 2 * dB * ( dPtVt - ptLn.z * vtLn.z) + 4 * dDeltaSqR * ptLn.z * vtLn.z ; 
   vdCoef[2] = 4 * dPtVt * dPtVt + dB * ( 1 - vtLn.z * vtLn.z) + 2 * ( dSqDistPO + dDeltaSqR * vtLn.z * vtLn.z) ;
   vdCoef[3] = 4 * dPtVt ;
   vdCoef[4] = 1 ;

  // Risolvo l'equazione e ridimensiono il vettore dei parametri
   DBLVECTOR vdRoots ;
   PolynomialRoots( 4, vdCoef, vdRoots) ;
   vdPar = vdRoots ;

  // Riordino le soluzioni
   for ( int ni = 0 ; ni < int( vdPar.size()) - 1 ; ++ ni) {
      for ( int nj = ni ; nj < int( vdPar.size()) ; ++ nj) {
         if ( vdPar[ni] > vdPar[nj]) {
            swap( vdPar[ni], vdPar[nj]) ;
         }
      }
   }

  // Studio le soluzioni
   int nIntType = T_ERROR ;
   if ( vdPar.empty()) 
      nIntType = T_NO_INT ;
   else if ( vdPar.size() == 1) {
      nIntType = T_ONE_TAN ;
      vbType.push_back( false) ;
   }
   else if ( vdPar.size() == 2) {
     // Soluzioni sufficientemente distinte
      if ( vdPar[1] - vdPar[0] > EPS_SMALL) {
         nIntType = T_TWO_SEC ;
         vbType.push_back( true) ;
         vbType.push_back( true) ;
      }        
     // Soluzioni coincidenti
      else {
         vdPar.resize( 1) ;
         nIntType = T_ONE_TAN ;
         vbType.push_back( false) ;
      }
   }
   else if ( vdPar.size() == 3) {
      nIntType = T_THREE_ONE_TAN_TWO_SEC ;
   }
   else if ( vdPar.size() == 4) {
     // Prime due soluzioni distinte
      if ( vdPar[1] - vdPar[0] > EPS_SMALL) {
        // Seconda e terza soluzione distinta
         if ( vdPar[2] - vdPar[1] > EPS_SMALL) {
           // Quarta e terza soluzione distinte
           // Qauttro secanti
            if ( vdPar[3] - vdPar[2] > EPS_SMALL) {
               nIntType = T_FOUR_SEC ;
               vbType.push_back( true) ;
               vbType.push_back( true) ;
               vbType.push_back( true) ;
               vbType.push_back( true) ;
            }
           // Quarta e terza soluzione coincidenti
           // Secante secante tangente
            else {
               nIntType = T_THREE_ONE_TAN_TWO_SEC ;
               vdPar.resize( 3) ;
               vbType.push_back( true) ;
               vbType.push_back( true) ;
               vbType.push_back( false) ;
            }
         }
        // Seconda e terza soluzione coincidenti
         else {
           // Secante tangente secante
            if ( vdPar[3] - vdPar[2] > EPS_SMALL) {
               nIntType = T_THREE_ONE_TAN_TWO_SEC ;
               vdPar[2] = vdPar[3] ;
               vdPar.resize( 3) ;
               vbType.push_back( true) ;
               vbType.push_back( false) ;
               vbType.push_back( true) ;
            }
         }
      }
     // Prime due soluzioni conicidenti
      else {
        // Terza e quarta soluzione distinte
        // Tangente secante secante
         if ( vdPar[3] - vdPar[2] > EPS_SMALL) {
            nIntType = T_THREE_ONE_TAN_TWO_SEC ;
            vdPar[1] = vdPar[2] ;
            vdPar[2] = vdPar[3] ;
            vdPar.resize( 3) ;
            vbType.push_back( false) ;
            vbType.push_back( true) ;
            vbType.push_back( true) ;
         }
        // Terza e quarta soluzione distinte
         else {
           // Secona e terza soluzione distinte 
           // Tanbente tangente
            if ( vdPar[2] - vdPar[1] > EPS_SMALL) {
               nIntType = T_TWO_TAN ;
               vdPar[1] = vdPar[2] ;
               vdPar.resize( 2) ;
               vbType.push_back( false) ;
               vbType.push_back( false) ;
            }
           // Seconda e terza soluzione coincidenti
           // Una soluzione tangente
            else {
               nIntType = T_ONE_TAN ;
               vdPar[0] = 0.5 * ( vdPar[0] + vdPar[3]) ;
               vdPar.resize( 1) ;
               vbType.push_back( false) ;
            }
         }
      }
   }
   return nIntType ; 
}

//----------------------------------------------------------------------------
int
RayTorus( const Point3d& ptLine, const Vector3d& vtLine,
          const Point3d& ptOTorus, const Vector3d& vtAxTorus, double dMinRad, double dMaxRad,
          BOOLVECTOR& vbType, DBLVECTOR& vdPar)
{
   int nIntType = LineTorus( ptLine, vtLine, ptOTorus, vtAxTorus, dMinRad, dMaxRad, vbType, vdPar) ;
   if ( nIntType == T_ERROR  ||  nIntType == T_NO_INT)
      return nIntType ;

   int nSize = int( vdPar.size()) ;
  // Ciclo sui parametri
   for ( int n = 0 ; n < nSize ; ++ n) {
     // Se il parametro è negativo lo elimino
      if ( vdPar[n] < - EPS_SMALL) {
         for ( int m = n ; m < nSize - 1 ; ++ m) {
            vdPar[m] = vdPar[m+1] ;
            vbType[m] = vbType[m+1] ;
         }
        // Ridimensiono
         -- nSize ;
         vdPar.resize( nSize) ;
         vbType.resize( nSize) ;
        // Aggiorno l'indice
         -- n ;
      }
   }

   int nSecNum = 0 ;
  // Ciclo sui parametri
   for ( int n= 0 ; n < nSize ; ++ n) {
     // Se secante ne aumento il numero
      if ( vbType[n])
         ++ nSecNum ;
   }
  
  // Numero di puti di tangenza
   int nNumTan = nSize - nSecNum ;
  // Nessun contatto
   if ( nSize == 0)
      nIntType = T_NO_INT ;
   else if ( nSize == 1) {
     // Semi-retta secante in un punto
      if ( vbType[0])
         nIntType = T_ONE_SEC ;
      else 
        // Semi-retta tangente in un punto
         nIntType = T_ONE_TAN ;
   }
   else if ( nSize == 2) {     
      // Semi-retta secante in due punti
      if ( nNumTan == 0) 
         nIntType = T_TWO_SEC ;
     // Semi-retta secante in un punto e tangente in un altro
      else if ( nNumTan == 1) 
         nIntType = T_TWO_TAN_SEC    ;
     // Semi-retta tangente in due punti
      else 
         nIntType = T_TWO_TAN ;
   }
   else if ( nSize == 3) {
     // Semi-retta secante in tre punti
      if ( nNumTan == 0)
         nIntType = T_THREE_SEC ;
     // Semi-retta secante in due punti e tangente in un punto
      else if ( nNumTan == 1) 
         nIntType = T_THREE_ONE_TAN_TWO_SEC ;
     // Errore
      else 
         nIntType = T_ERROR ;
   }
   else if ( nSize == 4) {
     // Semi-retta secante in quattro punti
      if ( nNumTan == 0)
         nIntType = T_FOUR_SEC ;
     // Errore
      else 
         nIntType = T_ERROR ;
   }
  // Errore
   else 
      nIntType = T_ERROR ;

   return nIntType ;
}

//----------------------------------------------------------------------------
int 
SegmentTorus( const Point3d& ptLine, const Vector3d& vtLine, double dSgLen,
              const Point3d& ptOTorus, const Vector3d& vtAxTorus, double dMinRad, double dMaxRad,
              BOOLVECTOR& vbType, DBLVECTOR& vdPar)
{
  // Verifico se il segmento interseca il box del toro
   Frame3d frTorus ; frTorus.Set( ptOTorus, vtAxTorus) ;
   Point3d ptLineL = ptLine ; ptLineL.ToLoc( frTorus) ;
   Vector3d vtLineL = vtLine ; vtLineL.ToLoc( frTorus) ;
   double dTotRad = dMaxRad + dMinRad ;
   double dU1, dU2 ;
   if ( ! IntersLineBox( ptLineL, vtLineL, Point3d( -dTotRad, -dTotRad, -dMinRad), Point3d( dTotRad, dTotRad, dMinRad), dU1, dU2)  ||
        dU2 < -EPS_SMALL  ||  dU1 > dSgLen + EPS_SMALL)
      return T_NO_INT ;

  // Calcolo l'intersezione con la linea illimitata
   int nIntType = LineTorus( ptLine, vtLine, ptOTorus, vtAxTorus, dMinRad, dMaxRad, vbType, vdPar) ;
   if ( nIntType == T_ERROR  ||  nIntType == T_NO_INT)
      return nIntType ;

  // Elimino le intersezioni fuori dal segmento
   int nSize = int( vdPar.size()) ;
   for ( int n = 0 ; n < nSize ; ++ n) {
      if ( vdPar[n] < -EPS_SMALL  ||  vdPar[n] > dSgLen + EPS_SMALL) {
        // tolgo l'intersezione dai vettori delle intersezioni
         vdPar.erase( vdPar.begin() + n) ;
         vbType.erase( vbType.begin() + n) ;
        // aggiorno la dimensione dei vettori
         -- nSize ;
        // aggiorno l'indice
         -- n ;
      }
   }

  // Determino il numero delle intersezioni secanti
   int nSecNum = 0 ;
   for ( int n = 0 ; n < nSize ; ++ n) {
      if ( vbType[n])
         ++ nSecNum ;
   }
  
  // Numero di punti di tangenza
   int nNumTan = nSize - nSecNum ;
  // Nessun contatto
   if ( nSize == 0)
      nIntType = T_NO_INT ;
   else if ( nSize == 1) {
     // Segmento secante in un punto
      if ( vbType[0])
         nIntType = T_ONE_SEC ;
      else 
        // Segmento tangente in un punto
         nIntType = T_ONE_TAN ;
   }
   else if ( nSize == 2) {     
      // Segmento secante in due punti
      if ( nNumTan == 0) 
         nIntType = T_TWO_SEC ;
     // Segmento secante in un punto e tangente in un altro
      else if ( nNumTan == 1) 
         nIntType = T_TWO_TAN_SEC    ;
     // Segmento tangente in due punti
      else 
         nIntType = T_TWO_TAN ;
   }
   else if ( nSize == 3) {
     // Segmento secante in tre punti
      if ( nNumTan == 0)
         nIntType = T_THREE_SEC ;
     // Segmento secante in due punti e tangente in un punto
      else if ( nNumTan == 1) 
         nIntType = T_THREE_ONE_TAN_TWO_SEC ;
     // Errore
      else 
         nIntType = T_ERROR ;
   }
   else if ( nSize == 4) {
     // Segmento secante in quattro punti
      if ( nNumTan == 0)
         nIntType = T_FOUR_SEC ;
     // Errore
      else 
         nIntType = T_ERROR ;
   }
  // Errore
   else 
      nIntType = T_ERROR ;

   return nIntType ;
}

//----------------------------------------------------------------------------
// Valuta la posizione reciproca fra un componente lineare e la parte inferiore 
// ed esterna rispetto al centro di una superficie torica.
int 
LinCompTorusExtInt( const Point3d& ptLine, const Vector3d& vtLine, double dSgLen, int nLinType, 
                    const Point3d& ptOTorus, const Vector3d& vtAxTorus, double dMinRad, double dMaxRad,                                            
                    BOOLVECTOR& vbType, DBLVECTOR& vdPar) 
{
   int nIntType ;
  // Retta
   if ( nLinType == Line) 
      nIntType = LineTorus( ptLine, vtLine, ptOTorus, vtAxTorus, dMinRad, dMaxRad, vbType, vdPar) ;
  // Semi-retta
   else if ( nLinType == Ray) 
      nIntType = RayTorus( ptLine, vtLine, ptOTorus, vtAxTorus, dMinRad, dMaxRad, vbType, vdPar) ;
  // Segmento
   else if ( nLinType == Segment) 
      nIntType = SegmentTorus( ptLine, vtLine, dSgLen, ptOTorus, vtAxTorus, dMinRad, dMaxRad, vbType, vdPar) ;
  // Errore
   else 
      nIntType = T_ERROR ;
  // Ritorno errore
   if ( nIntType == T_ERROR)
      return nIntType ;

   int nSize = int( vdPar.size()) ;
  // Ciclo sui punti di contatto
   for ( int n = 0 ; n < nSize ; ++ n) {
      Vector3d vtInt = ptLine + vdPar[n] * vtLine - ptOTorus ;
      Vector3d vtIntPlane = vtInt - vtInt * vtAxTorus * vtAxTorus ;
      vtInt.Normalize() ;
     // Se il punto di contatto è fuori dalla regione ammissibile lo elimino
      if ( vtIntPlane.Len() < dMaxRad - EPS_SMALL  ||
           vtInt * vtAxTorus > EPS_ZERO) {
         for ( int m = n ; m < nSize - 1 ; ++ m) {
            vdPar[m] = vdPar[m+1] ;
            vbType[m] = vbType[m+1] ;
         }
        // Ridimensiono
         -- nSize ;
         vdPar.resize( nSize) ;
         vbType.resize( nSize) ;
        // Aggiorno l'indice
         -- n ;
      }
   }

   int nSecNum = 0 ;
  // Ciclo sui parametri
   for ( int n = 0 ; n < nSize ; ++ n) {
     // Se secante ne aumento il numero
      if ( vbType[n])
         ++ nSecNum ;
   }
  
  // Numero di puti di tangenza
   int nNumTan = nSize - nSecNum ;
  // Nessun contatto
   if ( nSize == 0)
      nIntType = T_NO_INT ;
   else if ( nSize == 1) {
     // Segmento secante in un punto
      if ( vbType[0])
         nIntType = T_ONE_SEC ;
      else 
        // Segmento tangente in un punto
         nIntType = T_ONE_TAN ;
   }
   else if ( nSize == 2) {     
      // Segmento secante in due punti
      if ( nNumTan == 0) 
         nIntType = T_TWO_SEC ;
     // Segmento secante in un punto e tangente in un altro
      else if ( nNumTan == 1) 
         nIntType = T_TWO_TAN_SEC    ;
     // Segmento tangente in due punti
      else 
         nIntType = T_TWO_TAN ;
   }
   else if ( nSize == 3) {
     // Segmento secante in tre punti
      if ( nNumTan == 0)
         nIntType = T_THREE_SEC ;
     // Segmento secante in due punti e tangente in un punto
      else if ( nNumTan == 1) 
         nIntType = T_THREE_ONE_TAN_TWO_SEC ;
     // Errore
      else 
         nIntType = T_ERROR ;
   }
   else if ( nSize == 4) {
     // Segmento secante in quattro punti
      if ( nNumTan == 0)
         nIntType = T_FOUR_SEC ;
     // Errore
      else 
         nIntType = T_ERROR ;
   }
  // Errore
   else 
      nIntType = T_ERROR ;

   return nIntType ;
}

//----------------------------------------------------------------------------
// Valuta la posizione reciproca fra un componente lineare e la parte superiore 
// ed interna di una superficie torica.
int 
LinCompTorusInnUpInt( const Point3d& ptLine, const Vector3d& vtLine, double dSgLen, int nLinType,  
                      const Point3d& ptOTorus, const Vector3d& vtAxTorus, double dMinRad, double dMaxRad,
                      BOOLVECTOR& vbType, DBLVECTOR& vdPar) 
{
   int nIntType ;
  // Retta
   if ( nLinType == Line)
      nIntType = LineTorus( ptLine, vtLine, ptOTorus, vtAxTorus, dMinRad, dMaxRad, vbType, vdPar) ;
  // Semi-retta
   else if ( nLinType == Ray)
      nIntType = RayTorus( ptLine, vtLine, ptOTorus, vtAxTorus, dMinRad, dMaxRad, vbType, vdPar) ;
  // Segmento
   else if ( nLinType == Segment)
      nIntType = SegmentTorus( ptLine, vtLine, dSgLen, ptOTorus, vtAxTorus, dMinRad, dMaxRad, vbType, vdPar) ;
  // Errore
   else
      nIntType = S_ERROR_INT ;
  // Ritorno errore
   if ( nIntType == S_ERROR_INT)
      return nIntType ;

   int nSize = int( vdPar.size()) ;
  // Ciclo sui punti di contatto
   for ( int n = 0 ; n < nSize ; ++ n) {
      Vector3d vtInt = ptLine + vdPar[n] * vtLine - ptOTorus ;
      Vector3d vtIntPlane = vtInt - vtInt * vtAxTorus * vtAxTorus ;
      vtInt.Normalize() ;
     // Se il punto di contatto è fuori dalla regione ammissibile lo elimino
      if ( vtIntPlane.Len() > dMaxRad + EPS_SMALL  ||
           vtInt * vtAxTorus < - EPS_ZERO) {
         for ( int m = n ; m < nSize - 1 ; ++ m) {
            vdPar[m] = vdPar[m+1] ;
            vbType[m] = vbType[m+1] ;
         }
        // Ridimensiono
         -- nSize ;
         vdPar.resize( nSize) ;
         vbType.resize( nSize) ;
        // Aggiorno l'indice
         -- n ;
      }
   }

   int nSecNum = 0 ;
  // Ciclo sui parametri
   for ( int n = 0 ; n < nSize ; ++ n) {
     // Se secante ne aumento il numero
      if ( vbType[n])
         ++ nSecNum ;
   }
  
  // Numero di puti di tangenza
   int nNumTan = nSize - nSecNum ;
  // Nessun contatto
   if ( nSize == 0)
      nIntType = T_NO_INT ;
   else if ( nSize == 1) {
     // Segmento secante in un punto
      if ( vbType[0])
         nIntType = T_ONE_SEC ;
      else
        // Segmento tangente in un punto
         nIntType = T_ONE_TAN ;
   }
   else if ( nSize == 2) {
      // Segmento secante in due punti
      if ( nNumTan == 0) 
         nIntType = T_TWO_SEC ;
     // Segmento secante in un punto e tangente in un altro
      else if ( nNumTan == 1)
         nIntType = T_TWO_TAN_SEC ;
     // Segmento tangente in due punti
      else
         nIntType = T_TWO_TAN ;
   }
   else if ( nSize == 3) {
     // Segmento secante in tre punti
      if ( nNumTan == 0)
         nIntType = T_THREE_SEC ;
     // Segmento secante in due punti e tangente in un punto
      else if ( nNumTan == 1) 
         nIntType = T_THREE_ONE_TAN_TWO_SEC ;
     // Errore
      else
         nIntType = T_ERROR ;
   }
   else if ( nSize == 4) {
     // Segmento secante in quattro punti
      if ( nNumTan == 0)
         nIntType = T_FOUR_SEC ;
     // Errore
      else
         nIntType = T_ERROR ;
   }
  // Errore
   else
      nIntType = T_ERROR ;

   return nIntType ;
}