break;
             //
             case TopAbs_EDGE: {
-              Standard_Boolean bHasSameBounds;
-              Standard_Integer aNbComPrt2;
-              //
-              aNbComPrt2=aCPart.Ranges2().Length();
-              if (aNbComPrt2>1){
+              if (aNbCPrts > 1) {
                 break;
               }
-              //// <-LXBR   
+              //
+              Standard_Boolean bHasSameBounds;
               bHasSameBounds=aPB1->HasSameBounds(aPB2);
               if (!bHasSameBounds) {
                 break;
 
             }
             //
             const gp_Pnt& aPnew = BRep_Tool::Pnt(aVnew);
-            if (!myContext->IsValidPointForFace(aPnew, aF, aTolE)) {
+            if (!myContext->IsValidPointForFace(aPnew, aF, aTolE+aTolF)) {
               continue;
             }
             //
 
     case GeomAbs_Cylinder: quad1.SetValue(HS1->Surface().Cylinder()); break;
     case GeomAbs_Cone:     quad1.SetValue(HS1->Surface().Cone()); break;
     case GeomAbs_Sphere:   quad1.SetValue(HS1->Surface().Sphere()); break;
+    case GeomAbs_Torus:    quad1.SetValue(HS1->Surface().Torus()); break;
     default: Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
   }
 }
 
     case GeomAbs_Cylinder: quad1.SetValue(myHS1->Surface().Cylinder()); break;
     case GeomAbs_Cone:     quad1.SetValue(myHS1->Surface().Cone()); break;
     case GeomAbs_Sphere:   quad1.SetValue(myHS1->Surface().Sphere()); break;
+    case GeomAbs_Torus:    quad1.SetValue(myHS1->Surface().Torus()); break;
     default: Standard_ConstructionError::Raise("GeomInt_LineConstructor::Parameters");
   }
   switch (myHS2->Surface().GetType())
     case GeomAbs_Cylinder: quad2.SetValue(myHS2->Surface().Cylinder()); break;
     case GeomAbs_Cone:     quad2.SetValue(myHS2->Surface().Cone()); break;
     case GeomAbs_Sphere:   quad2.SetValue(myHS2->Surface().Sphere()); break;
+    case GeomAbs_Torus:    quad2.SetValue(myHS2->Surface().Torus()); break;
     default: Standard_ConstructionError::Raise("GeomInt_LineConstructor::Parameters");
   }
   quad1.Parameters(Ptref,U1,V1);
 
 uses Pln         from gp,
      Cylinder    from gp,
      Cone        from gp,
-     Sphere      from gp,
+     Sphere      from gp, 
+     Torus       from gp,
      Pnt         from gp,
      Lin         from gp,
      Circ        from gp,
      
      
 raises  NotDone      from StdFail,
-       DomainError  from Standard,
-        OutOfRange   from Standard
+    DomainError  from Standard,
+    OutOfRange   from Standard
 is
 
     Create
-    
-       ---Purpose: Empty constructor.
-       
-       returns QuadQuadGeo from IntAna;
+    ---Purpose: Empty constructor.
+    returns QuadQuadGeo from IntAna;
+
 
-       
     Create(P1,P2        : Pln    from gp; 
            TolAng, Tol  : Real   from Standard)
-    
-       ---Purpose: Creates the intersection between two planes.
-       --          TolAng is the angular tolerance used to determine
-       --          if the planes are parallel.
-       --          Tol is the tolerance used to determine if the planes
-       --          are identical (only when they are parallel).
-    
-       returns QuadQuadGeo from IntAna;
+    ---Purpose: Creates the intersection between two planes.
+    --          TolAng is the angular tolerance used to determine
+    --          if the planes are parallel.
+    --          Tol is the tolerance used to determine if the planes
+    --          are identical (only when they are parallel).
+    returns QuadQuadGeo from IntAna;
 
 
     Perform(me          : in out; 
             P1,P2       : Pln    from gp; 
             TolAng, Tol : Real   from Standard)
-    
-       ---Purpose: Intersects two planes.
-       --          TolAng is the angular tolerance used to determine
-       --          if the planes are parallel.
-       --          Tol is the tolerance used to determine if the planes
-       --          are identical (only when they are parallel).
-    
-       is static;
+    ---Purpose: Intersects two planes.
+    --          TolAng is the angular tolerance used to determine
+    --          if the planes are parallel.
+    --          Tol is the tolerance used to determine if the planes
+    --          are identical (only when they are parallel).
+    is static;
 
 
     Create(P : Pln      from gp; 
     --          center is less than Tol, the result will be the circle. 
     --          H is the height of the cylinder <Cyl>. It is  used to check
     --          whether the plane and cylinder are parallel.
-
     returns QuadQuadGeo from IntAna;
-       
+
 
     Perform(me: in out; 
             P : Pln      from gp; 
             C : Cylinder from gp; 
             Tolang,Tol: Real from Standard; 
             H :Real from Standard = 0)
-    
     ---Purpose: Intersects a plane and a cylinder.
     --          TolAng is the angular tolerance used to determine
     --          if the axis of the cylinder is parallel to the plane.
     --          center is less than Tol, the result will be the circle.
     --          H is the height of the cylinder <Cyl>. It is  used to check
     --          whether the plane and cylinder are parallel.
-
     is static;
-       
+
 
     Create(P : Pln    from gp; 
            S : Sphere from gp)
+    ---Purpose: Creates the intersection between a plane and a sphere.
+    returns QuadQuadGeo from IntAna;
 
-       ---Purpose: Creates the intersection between a plane and a sphere.
 
-       returns QuadQuadGeo from IntAna;
-       
-       
     Perform(me: in out; 
             P : Pln    from gp; 
             S : Sphere from gp)
+    ---Purpose: Intersects a plane and a sphere. 
+    is static;
 
-       ---Purpose: Intersects a plane and a sphere.
 
-       is static;
-       
-       
     Create(P : Pln  from gp; 
            C : Cone from gp; 
            Tolang,Tol: Real from Standard)
-    
-       ---Purpose: Creates the intersection between a plane and a cone.
-       --          TolAng is the angular tolerance used to determine
-       --          if the axis of the cone is parallel or perpendicular
-       --          to the plane, and if the generating line of the cone
-       --          is parallel to the plane.
-       --          Tol is the tolerance used to determine if the apex
-       --          of the cone is in the plane.
+    ---Purpose: Creates the intersection between a plane and a cone.
+    --          TolAng is the angular tolerance used to determine
+    --          if the axis of the cone is parallel or perpendicular
+    --          to the plane, and if the generating line of the cone
+    --          is parallel to the plane.
+    --          Tol is the tolerance used to determine if the apex
+    --          of the cone is in the plane.
+    returns QuadQuadGeo from IntAna;
 
-       returns QuadQuadGeo from IntAna;
 
-               
     Perform(me: in out; 
             P : Pln  from gp; 
             C : Cone from gp; 
             Tolang,Tol: Real from Standard)
-    
-       ---Purpose: Intersects a plane and a cone.
-       --          TolAng is the angular tolerance used to determine
-       --          if the axis of the cone is parallel or perpendicular
-       --          to the plane, and if the generating line of the cone
-       --          is parallel to the plane.
-       --          Tol is the tolerance used to determine if the apex
-       --          of the cone is in the plane.
-
-       is static;
+    ---Purpose: Intersects a plane and a cone.
+    --          TolAng is the angular tolerance used to determine
+    --          if the axis of the cone is parallel or perpendicular
+    --          to the plane, and if the generating line of the cone
+    --          is parallel to the plane.
+    --          Tol is the tolerance used to determine if the apex
+    --          of the cone is in the plane.
+    is static;
 
 
     Create(Cyl1,Cyl2: Cylinder from gp; 
            Tol      : Real     from Standard)
+    ---Purpose: Creates the intersection between two cylinders.
+    returns QuadQuadGeo from IntAna;
 
-       ---Purpose: Creates the intersection between two cylinders.
-    
-       returns QuadQuadGeo from IntAna;
-       
 
     Perform(me       : in out; 
             Cyl1,Cyl2: Cylinder from gp; 
             Tol      : Real     from Standard)
+    ---Purpose: Intersects two cylinders
+    is static;
 
-       ---Purpose: Intersects two cylinders
-    
-       is static;
-       
 
     Create(Cyl: Cylinder from gp;  
            Sph: Sphere   from gp; 
            Tol: Real     from Standard)
-          
-       ---Purpose: Creates the intersection between a Cylinder and a Sphere.
-    
-        returns QuadQuadGeo from IntAna;
-       
-       
+    ---Purpose: Creates the intersection between a Cylinder and a Sphere.
+    returns QuadQuadGeo from IntAna;
+
+
     Perform(me : in out; 
             Cyl: Cylinder from gp;  
             Sph: Sphere   from gp; 
             Tol: Real     from Standard)
+    ---Purpose: Intersects a cylinder and a sphere.
+    is static;
 
-       ---Purpose: Intersects a cylinder and a sphere.
 
-       is static;
-       
-       
     Create(Cyl: Cylinder from gp;  
            Con: Cone     from gp; 
            Tol: Real     from Standard)
+    ---Purpose: Creates the intersection between a Cylinder and a Cone
+    returns QuadQuadGeo from IntAna;
 
-        ---Purpose: Creates the intersection between a Cylinder and a Cone
-            
-        returns QuadQuadGeo from IntAna;
 
-       
     Perform(me : in out; 
             Cyl: Cylinder from gp;  
             Con: Cone     from gp; 
             Tol: Real     from Standard)
+    ---Purpose: Intersects a cylinder and a cone.
+    is static;
 
-       ---Purpose: Intersects a cylinder and a cone.
-
-       is static;
 
-       
     Create(Sph1: Sphere from gp;  
            Sph2: Sphere from gp; 
            Tol : Real   from Standard)
+    ---Purpose: Creates the intersection between two Spheres.    
+    returns QuadQuadGeo from IntAna;
 
-        ---Purpose: Creates the intersection between two Spheres.    
 
-        returns QuadQuadGeo from IntAna;
-       
-       
     Perform(me  : in out; 
             Sph1: Sphere from gp;  
             Sph2: Sphere from gp; 
             Tol : Real   from Standard)
+    ---Purpose: Intersects a two spheres.
+    is static;
 
-       ---Purpose: Intersects a two spheres.
 
-       is static;
-       
-       
     Create(Sph: Sphere from gp;  
            Con: Cone   from gp; 
            Tol: Real   from Standard)
-    
-        ---Purpose: Creates the intersection beween a Sphere and a Cone.
-    
-        returns QuadQuadGeo from IntAna;
+    ---Purpose: Creates the intersection beween a Sphere and a Cone.
+    returns QuadQuadGeo from IntAna;
+
 
-       
     Perform(me : in out; 
             Sph: Sphere from gp;  
             Con: Cone   from gp; 
             Tol: Real   from Standard)
-
-       ---Purpose: Intersects a sphere and a cone.
-
-       is static;
+    ---Purpose: Intersects a sphere and a cone.
+    is static;
 
 
-    Create(Con1: Cone from gp;  Con2: Cone from gp; Tol:Real from Standard)
-    
-        ---Purpose: Creates the intersection beween two cones.
-    
-        returns QuadQuadGeo from IntAna;
+    Create(Con1: Cone from gp; 
+           Con2: Cone from gp;  
+           Tol : Real from Standard)
+    ---Purpose: Creates the intersection beween two cones.
+    returns QuadQuadGeo from IntAna;
 
 
     Perform(me  : in out; 
             Con1: Cone from gp;  
             Con2: Cone from gp; 
             Tol :Real from Standard)
+    ---Purpose: Intersects two cones.
+    is static;
+
+
+    Create(Pln : Pln   from gp;
+           Tor : Torus from gp;
+           Tol : Real  from Standard)
+    ---Purpose: Creates the intersection beween plane and torus.
+    returns QuadQuadGeo from IntAna;
 
-       ---Purpose: Intersects two cones.
 
-       is static;
+    Perform(me  :in out; 
+            Pln : Pln   from gp;
+            Tor : Torus from gp;
+            Tol : Real  from Standard)
+    ---Purpose: Intersects plane and torus.
+    is static;
 
 
-    IsDone(me)
+    Create(Cyl : Cylinder from gp;
+           Tor : Torus    from gp;
+           Tol : Real     from Standard)
+    ---Purpose: Creates the intersection beween cylinder and torus.
+    returns QuadQuadGeo from IntAna;
+
+
+    Perform(me  : in out; 
+            Cyl : Cylinder from gp;
+            Tor : Torus    from gp;
+            Tol : Real     from Standard)
+    ---Purpose: Intersects cylinder and torus.
+    is static;
+
+
+    Create(Con : Cone  from gp;
+           Tor : Torus from gp;
+           Tol : Real  from Standard)
+    ---Purpose: Creates the intersection beween cone and torus.
+    returns QuadQuadGeo from IntAna;
+
+
+    Perform(me  : in out; 
+            Con : Cone  from gp;
+            Tor : Torus from gp;
+            Tol : Real  from Standard)
+    ---Purpose: Intersects cone and torus.
+    is static;
+
+
+    Create(Sph : Sphere from gp;
+           Tor : Torus  from gp;
+           Tol : Real   from Standard)
+    ---Purpose: Creates the intersection beween sphere and torus.
+    returns QuadQuadGeo from IntAna;
+
+
+    Perform(me  : in out; 
+            Sph : Sphere from gp;
+            Tor : Torus  from gp;
+            Tol : Real   from Standard)
+    ---Purpose: Intersects sphere and torus.
+    is static;
 
-       ---Purpose: Returns Standard_True if the computation was successful.
-       --          
-       ---C++: inline
 
-       returns Boolean from Standard
-       is static;
+    Create(Tor1 : Torus from gp;
+           Tor2 : Torus from gp;
+           Tol  : Real  from Standard)
+    ---Purpose: Creates the intersection beween two toruses.
+    returns QuadQuadGeo from IntAna;
+
+
+    Perform(me   : in out; 
+            Tor1 : Torus from gp;
+            Tor2 : Torus from gp;
+            Tol  : Real  from Standard)
+    ---Purpose: Intersects two toruses.
+    is static;
+
+
+    IsDone(me)
+    ---Purpose: Returns Standard_True if the computation was successful.
+    --          
+    ---C++: inline
+    returns Boolean from Standard
+    is static;
 
 
     TypeInter(me)
-    
-       ---Purpose: Returns the type of intersection.
-       --          
-       ---C++: inline
-    
-       returns ResultType from IntAna
-       
-       raises NotDone from StdFail
-       --- The exception NotDone is raised if IsDone return Standard_False.
-       
-       is static;
-       
-       
-     NbSolutions(me) 
-     
-       ---Purpose: Returns the number of interesections.
-       --          The possible intersections are :
-       --           - 1 point
-       --           - 1 or 2 line(s)
-       --           - 1 Point and 1 Line
-       --           - 1 circle
-       --           - 1 ellipse
-       --           - 1 parabola
-       --           - 1 or 2 hyperbola(s).
-       --          
-       ---C++: inline
-     
-        returns Integer from Standard
-       
-       raises NotDone from StdFail
-       --- The exception NotDone is raised if IsDone returns Standard_False.
-       
-       is static;
-        
-       
-     Point(me; Num: Integer from Standard)
-     
-       ---Purpose: Returns the point solution of range Num.
-     
-       returns Pnt from gp
-       
-       raises DomainError from Standard,
-              OutOfRange  from Standard,
-              NotDone     from StdFail
-       --- The exception NotDone is raised if IsDone return Standard_False.
-       --  The exception DomainError is raised if TypeInter does not return
-       --  IntAna_Point or TypeInter does not return IntAna_PointAndCircle.
-       --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
-
-       is static;
-       
-       
-     Line(me; Num: Integer from Standard)
-     
-       ---Purpose: Returns the line solution of range Num.
-     
-       returns Lin from gp
-       
-       raises DomainError from Standard,
-              OutOfRange  from Standard,
-              NotDone     from StdFail
-       --- The exception NotDone is raised if IsDone return Standard_False.
-       --  The exception DomainError is raised if TypeInter does not return
-       --  IntAna_Line.
-       --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
-       
-       is static;
-       
-       
-     Circle(me; Num: Integer from Standard)
-     
-       ---Purpose: Returns the circle solution of range Num.
-     
-       returns Circ from gp
-       
-       raises DomainError from Standard,
-              OutOfRange  from Standard,
-              NotDone     from StdFail
-       --- The exception NotDone is raised if IsDone return Standard_False.
-       --  The exception DomainError is raised if TypeInter does not return
-       --  IntAna_Circle or TypeInter does not return IntAna_PointAndCircle.
-       --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
-
-       is static;
-       
-       
-     Ellipse(me; Num: Integer from Standard)
-
-       ---Purpose: Returns the ellipse solution of range Num.
-     
-       returns Elips from gp
-       
-       raises DomainError from Standard,
-              OutOfRange  from Standard,
-              NotDone     from StdFail
-       --- The exception NotDone is raised if IsDone return Standard_False.
-       --  The exception DomainError is raised if TypeInter does not return
-       --  IntAna_Ellipse.
-       --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
-
-       is static;
-          
-
-     Parabola(me; Num: Integer from Standard)
-
-       ---Purpose: Returns the parabola solution of range Num.
-     
-       returns Parab from gp
-       
-       raises DomainError from Standard,
-              OutOfRange  from Standard,
-              NotDone     from StdFail
-       --- The exception NotDone is raised if IsDone return Standard_False.
-       --  The exception DomainError is raised if TypeInter does not return
-       --  IntAna_Parabola.
-       --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
+    ---Purpose: Returns the type of intersection.
+    --          
+    ---C++: inline
+    returns ResultType from IntAna
+    raises NotDone from StdFail
+    -- The exception NotDone is raised if IsDone return Standard_False.
+    is static;
 
-       is static;
 
+    NbSolutions(me) 
+    ---Purpose: Returns the number of interesections.
+    --          The possible intersections are :
+    --           - 1 point
+    --           - 1 or 2 line(s)
+    --           - 1 Point and 1 Line
+    --           - 1 circle
+    --           - 1 ellipse
+    --           - 1 parabola
+    --           - 1 or 2 hyperbola(s).
+    --          
+    ---C++: inline
+    returns Integer from Standard
+    raises NotDone from StdFail
+    -- The exception NotDone is raised if IsDone returns Standard_False.
+    is static;
+ 
+
+    Point(me; Num: Integer from Standard)
+    ---Purpose: Returns the point solution of range Num.
+    returns Pnt from gp
+    raises DomainError from Standard,
+           OutOfRange  from Standard,
+           NotDone     from StdFail
+    --  The exception NotDone is raised if IsDone return Standard_False.
+    --  The exception DomainError is raised if TypeInter does not return
+    --  IntAna_Point or TypeInter does not return IntAna_PointAndCircle.
+    --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
+    is static;
 
-    Hyperbola(me; Num: Integer from Standard)
 
-       ---Purpose: Returns the hyperbola solution of range Num.
-     
-       returns Hypr from gp
-       
-       raises DomainError from Standard,
-              OutOfRange  from Standard,
-              NotDone     from StdFail
-       --- The exception NotDone is raised if IsDone return Standard_False.
-       --  The exception DomainError is raised if TypeInter does not return
-       --  IntAna_Hyperbola.
-       --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
-
-       is static; 
-        
-    HasCommonGen(me)  returns  Boolean  from  Standard;         
+    Line(me; Num: Integer from Standard)
+    ---Purpose: Returns the line solution of range Num.
+    returns Lin from gp
+    raises DomainError from Standard,
+           OutOfRange  from Standard,
+           NotDone     from StdFail
+    --  The exception NotDone is raised if IsDone return Standard_False.
+    --  The exception DomainError is raised if TypeInter does not return
+    --  IntAna_Line.
+    --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
+    is static;
+
+
+    Circle(me; Num: Integer from Standard)
+    ---Purpose: Returns the circle solution of range Num.
+    returns Circ from gp
+    raises DomainError from Standard,
+           OutOfRange  from Standard,
+           NotDone     from StdFail
+    --  The exception NotDone is raised if IsDone return Standard_False.
+    --  The exception DomainError is raised if TypeInter does not return
+    --  IntAna_Circle or TypeInter does not return IntAna_PointAndCircle.
+    --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
+    is static;
+
+
+    Ellipse(me; Num: Integer from Standard)
+    ---Purpose: Returns the ellipse solution of range Num.
+    returns Elips from gp
+    raises DomainError from Standard,
+           OutOfRange  from Standard,
+           NotDone     from StdFail
+    --  The exception NotDone is raised if IsDone return Standard_False.
+    --  The exception DomainError is raised if TypeInter does not return
+    --  IntAna_Ellipse.
+    --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
+    is static;
+
+
+    Parabola(me; Num: Integer from Standard)
+    ---Purpose: Returns the parabola solution of range Num.
+    returns Parab from gp
+    raises DomainError from Standard,
+           OutOfRange  from Standard,
+           NotDone     from StdFail
+    --  The exception NotDone is raised if IsDone return Standard_False.
+    --  The exception DomainError is raised if TypeInter does not return
+    --  IntAna_Parabola.
+    --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
+    is static;
+
+
+    Hyperbola(me; Num: Integer from Standard)
+    ---Purpose: Returns the hyperbola solution of range Num.
+    returns Hypr from gp
+    raises DomainError from Standard,
+           OutOfRange  from Standard,
+           NotDone     from StdFail
+    --  The exception NotDone is raised if IsDone return Standard_False.
+    --  The exception DomainError is raised if TypeInter does not return
+    --  IntAna_Hyperbola.
+    --  The exception OutOfRange is raised if Num < 1 or Num > NbSolutions.
+    is static; 
+ 
+    HasCommonGen(me)  returns  Boolean  from  Standard; 
     PChar(me)  returns  Pnt  from  gp; 
     ---C++:  return  const&
-          
+
     InitTolerances(me:out)  
-        
-       ---Purpose: Initialize the values of inner tolerances. 
-       
-       is protected; 
-        
+    ---Purpose: Initialize the values of inner tolerances. 
+    is protected; 
+ 
 fields
 
     done       :   Boolean      from Standard  is protected;
     
     pt1        :   Pnt          from gp        is protected;
     pt2        :   Pnt          from gp        is protected;
+    pt3        :   Pnt          from gp        is protected;
+    pt4        :   Pnt          from gp        is protected;
 
     dir1       :   Dir          from gp        is protected;
     dir2       :   Dir          from gp        is protected;
+    dir3       :   Dir          from gp        is protected;
+    dir4       :   Dir          from gp        is protected;
 
     param1     :   Real         from Standard  is protected;
     param2     :   Real         from Standard  is protected;
+    param3     :   Real         from Standard  is protected;
+    param4     :   Real         from Standard  is protected;
     param1bis  :   Real         from Standard  is protected;
     param2bis  :   Real         from Standard  is protected;
     -- 
     --  
     myCommonGen  :  Boolean    from  Standard  is  protected; 
     myPChar      :  Pnt        from  gp        is protected;
-               
+
 end QuadQuadGeo;
 
   AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
 
   void Distance(Standard_Real& dist,
-               Standard_Real& Param1,
-               Standard_Real& Param2);
+                Standard_Real& Param1,
+                Standard_Real& Param2);
   
   gp_Pnt PtIntersect()              { 
     return ptintersect;
 
  protected:
   Standard_Real Det33(const Standard_Real a11,
-                     const Standard_Real a12,
-                     const Standard_Real a13,
-                     const Standard_Real a21,
-                     const Standard_Real a22,
-                     const Standard_Real a23,
-                     const Standard_Real a31,
-                     const Standard_Real a32,
-                     const Standard_Real a33) {
+                      const Standard_Real a12,
+                      const Standard_Real a13,
+                      const Standard_Real a21,
+                      const Standard_Real a22,
+                      const Standard_Real a23,
+                      const Standard_Real a31,
+                      const Standard_Real a32,
+                      const Standard_Real a33) {
     Standard_Real theReturn =  
       a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;   
     return theReturn ;
   }
   else {   
     thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
-                                                          Axe2.Location())));
+                                                           Axe2.Location())));
   }
   //--- check if Axis are Coplanar
   Standard_Real D33;
   if(thedistance<myEPSILON_DISTANCE) {
     D33=Det33(V1.X(),V1.Y(),V1.Z()
-             ,V2.X(),V2.Y(),V2.Z()
-             ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
+              ,V2.X(),V2.Y(),V2.Z()
+              ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
     if(Abs(D33)<=myEPSILON_DISTANCE) { 
       thecoplanar=Standard_True;
     }
       A=(smy * V2.X() - smx * V2.Y())/Det1;
     }
     else if((Det2!=0.0) 
-           && ((Abs(Det2) >= Abs(Det1))
-               &&(Abs(Det2) >= Abs(Det3)))) {
+            && ((Abs(Det2) >= Abs(Det1))
+                &&(Abs(Det2) >= Abs(Det3)))) {
       A=(smz * V2.Y() - smy * V2.Z())/Det2;
     }
     else {
       A=(smz * V2.X() - smx * V2.Z())/Det3;
     }
     ptintersect.SetCoord( P1.X() + A * V1.X()
-                        ,P1.Y() + A * V1.Y()
-                        ,P1.Z() + A * V1.Z());
+                         ,P1.Y() + A * V1.Y()
+                         ,P1.Z() + A * V1.Z());
   }
   else { 
     ptintersect.SetCoord(0,0,0);  //-- Pour eviter des FPE
   
   gp_Dir N  = U1.Crossed(U2);
   Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
-                         U1.Y(),U2.Y(),N.Y(),
-                         U1.Z(),U2.Z(),N.Z());
+                          U1.Y(),U2.Y(),N.Y(),
+                          U1.Z(),U2.Z(),N.Z());
   if(D) { 
     dist = Det33(U1.X(),U2.X(),O1O2.X(),
-                U1.Y(),U2.Y(),O1O2.Y(),
-                U1.Z(),U2.Z(),O1O2.Z()) / D;
+                 U1.Y(),U2.Y(),O1O2.Y(),
+                 U1.Z(),U2.Z(),O1O2.Z()) / D;
     Param1 = Det33(O1O2.X(),U2.X(),N.X(),
-                  O1O2.Y(),U2.Y(),N.Y(),
-                  O1O2.Z(),U2.Z(),N.Z()) / (-D);
+                   O1O2.Y(),U2.Y(),N.Y(),
+                   O1O2.Z(),U2.Z(),N.Z()) / (-D);
     //------------------------------------------------------------
     //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
     //-- soit : Segment perpendiculaire : O1+P1 D1
     //--                                  O2-P2 D2
     Param2 = Det33(U1.X(),O1O2.X(),N.X(),
-                  U1.Y(),O1O2.Y(),N.Y(),
-                  U1.Z(),O1O2.Z(),N.Z()) / (D);
+                   U1.Y(),O1O2.Y(),N.Y(),
+                   U1.Z(),O1O2.Z(),N.Z()) / (D);
   }
 }
 //=======================================================================
       typeres(IntAna_Empty),
       pt1(0,0,0),
       pt2(0,0,0),
+      pt3(0,0,0),
+      pt4(0,0,0),
       param1(0),
       param2(0),
+      param3(0),
+      param4(0),
       param1bis(0),
       param2bis(0),
       myCommonGen(Standard_False),
 //purpose  : Pln  Pln 
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1, 
-                                        const gp_Pln& P2,
-                                        const Standard_Real TolAng,
-                                        const Standard_Real Tol)
+                                         const gp_Pln& P2,
+                                         const Standard_Real TolAng,
+                                         const Standard_Real Tol)
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1, 
-                                   const gp_Pln& P2,
-                                   const Standard_Real TolAng,
-                                   const Standard_Real Tol)
+                                    const gp_Pln& P2,
+                                    const Standard_Real TolAng,
+                                    const Standard_Real Tol)
 {
   done=Standard_False;
   //
       typeres(IntAna_Empty),
       pt1(0,0,0),
       pt2(0,0,0),
+      pt3(0,0,0),
+      pt4(0,0,0),
       param1(0),
       param2(0),
+      param3(0),
+      param4(0),
       param1bis(0),
       param2bis(0),
       myCommonGen(Standard_False),
       Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
       Standard_Real dangle = Abs( dang );
       if( dangle > Tolang )
-       {
-         Standard_Real sinda = Abs( Sin( dangle ) );
-         Standard_Real dif = Abs( sinda - Tol );
-         if( dif < Tol )
-           {
-             tolang = sinda * 2.;
-             newparams = Standard_True;
-           }
-       }
+        {
+          Standard_Real sinda = Abs( Sin( dangle ) );
+          Standard_Real dif = Abs( sinda - Tol );
+          if( dif < Tol )
+            {
+              tolang = sinda * 2.;
+              newparams = Standard_True;
+            }
+        }
     }
 
   nbint = 0;
 
     if (Abs(Abs(dist)-radius) < Tol)
       {
-       nbint = 1;
-       pt1.SetXYZ(omega);
-
-       if( newparams )
-         {
-           gp_XYZ omegaXYZ(X,Y,Z);
-           gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
-           Standard_Real Xt,Yt,Zt,distt;
-           omegaXYZtrnsl.Coord(Xt,Yt,Zt);
-           distt = A*Xt + B*Yt + C*Zt + D;
-           gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
-           gp_Pnt ppt1;
-           ppt1.SetXYZ( omega1 );
-           gp_Vec vv1(pt1,ppt1);
-           gp_Dir dd1( vv1 );
-           dir1 = dd1;
-         }
-       else
-         dir1 = axec.Direction();
+        nbint = 1;
+        pt1.SetXYZ(omega);
+
+        if( newparams )
+          {
+            gp_XYZ omegaXYZ(X,Y,Z);
+            gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
+            Standard_Real Xt,Yt,Zt,distt;
+            omegaXYZtrnsl.Coord(Xt,Yt,Zt);
+            distt = A*Xt + B*Yt + C*Zt + D;
+            gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
+            gp_Pnt ppt1;
+            ppt1.SetXYZ( omega1 );
+            gp_Vec vv1(pt1,ppt1);
+            gp_Dir dd1( vv1 );
+            dir1 = dd1;
+          }
+        else
+          dir1 = axec.Direction();
     }
     else if (Abs(dist) < radius)
       {
-       nbint = 2;
-       h = Sqrt(radius*radius - dist*dist);
-       axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
-
-       pt1.SetXYZ(omega - h*axey);
-       pt2.SetXYZ(omega + h*axey);
-
-       if( newparams )
-         { 
-           gp_XYZ omegaXYZ(X,Y,Z);
-           gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
-           Standard_Real Xt,Yt,Zt,distt,ht;
-           omegaXYZtrnsl.Coord(Xt,Yt,Zt);
-           distt = A*Xt + B*Yt + C*Zt + D;
-           //      ht = Sqrt(radius*radius - distt*distt);
-           Standard_Real anSqrtArg = radius*radius - distt*distt;
-           ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
-
-           gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
-           gp_Pnt ppt1,ppt2;
-           ppt1.SetXYZ( omega1 - ht*axey);
-           ppt2.SetXYZ( omega1 + ht*axey);
-           gp_Vec vv1(pt1,ppt1);
-           gp_Vec vv2(pt2,ppt2);
-           gp_Dir dd1( vv1 );
-           gp_Dir dd2( vv2 );
-           dir1 = dd1;
-           dir2 = dd2;
-         }
-       else
-         {
-           dir1 = axec.Direction();
-           dir2 = axec.Direction();
-         }
+        nbint = 2;
+        h = Sqrt(radius*radius - dist*dist);
+        axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
+
+        pt1.SetXYZ(omega - h*axey);
+        pt2.SetXYZ(omega + h*axey);
+
+        if( newparams )
+          { 
+            gp_XYZ omegaXYZ(X,Y,Z);
+            gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
+            Standard_Real Xt,Yt,Zt,distt,ht;
+            omegaXYZtrnsl.Coord(Xt,Yt,Zt);
+            distt = A*Xt + B*Yt + C*Zt + D;
+            //             ht = Sqrt(radius*radius - distt*distt);
+            Standard_Real anSqrtArg = radius*radius - distt*distt;
+            ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
+
+            gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
+            gp_Pnt ppt1,ppt2;
+            ppt1.SetXYZ( omega1 - ht*axey);
+            ppt2.SetXYZ( omega1 + ht*axey);
+            gp_Vec vv1(pt1,ppt1);
+            gp_Vec vv2(pt2,ppt2);
+            gp_Dir dd1( vv1 );
+            gp_Dir dd2( vv2 );
+            dir1 = dd1;
+            dir2 = dd2;
+          }
+        else
+          {
+            dir1 = axec.Direction();
+            dir2 = axec.Direction();
+          }
     }
     //  else nbint = 0
 
 //purpose  : Pln Cone
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
-                                        const gp_Cone& Co,
-                                        const Standard_Real Tolang,
-                                        const Standard_Real Tol) 
+                                         const gp_Cone& Co,
+                                         const Standard_Real Tolang,
+                                         const Standard_Real Tol) 
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
-                                  const gp_Cone& Co,
-                                  const Standard_Real Tolang,
-                                  const Standard_Real Tol) 
+                                   const gp_Cone& Co,
+                                   const Standard_Real Tolang,
+                                   const Standard_Real Tol) 
 {
 
   done = Standard_False;
       distance = apex.Distance(center);
 
       if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
-       axex.Reverse();
-       axey.Reverse();
+        axex.Reverse();
+        axey.Reverse();
       }
 
       if (Abs(costa) < Tolang) {  // plan parallele a une generatrice
-       typeres = IntAna_Parabola;
-       nbint = 1;
-       deltacenter = distance/2./cosa;
-       axex.Normalize();
-       pt1.SetXYZ(center.XYZ()-deltacenter*axex);
-       dir1 = normp;
-       dir2.SetXYZ(axex);
-       param1 = deltacenter*sina*sina;
+        typeres = IntAna_Parabola;
+        nbint = 1;
+        deltacenter = distance/2./cosa;
+        axex.Normalize();
+        pt1.SetXYZ(center.XYZ()-deltacenter*axex);
+        dir1 = normp;
+        dir2.SetXYZ(axex);
+        param1 = deltacenter*sina*sina;
       }
       else if (sint  < Tolang) {            // plan perpendiculaire a l axe
-       typeres = IntAna_Circle;
-       nbint = 1;
-       pt1 = center;
-       dir1 = Co.Position().Direction();
-       dir2 = Co.Position().XDirection();
-       param1 = apex.Distance(center)*Abs(Tan(angl));
+        typeres = IntAna_Circle;
+        nbint = 1;
+        pt1 = center;
+        dir1 = Co.Position().Direction();
+        dir2 = Co.Position().XDirection();
+        param1 = apex.Distance(center)*Abs(Tan(angl));
       }
       else if (cost < sina ) {
-       typeres = IntAna_Hyperbola;
-       nbint = 2;
-       axex.Normalize();
-
-       deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
-       pt1.SetXYZ(center.XYZ() - deltacenter*axex);
-       pt2 = pt1;
-       dir1 = normp;
-       dir2.SetXYZ(axex);
-       param1    = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
-       param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
+        typeres = IntAna_Hyperbola;
+        nbint = 2;
+        axex.Normalize();
+
+        deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
+        pt1.SetXYZ(center.XYZ() - deltacenter*axex);
+        pt2 = pt1;
+        dir1 = normp;
+        dir2.SetXYZ(axex);
+        param1    = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
+        param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
 
       }
       else {                    // on a alors cost > sina
-       typeres = IntAna_Ellipse;
-       nbint = 1;
-       Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
-       deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
-       axex.Normalize();
-       pt1.SetXYZ(center.XYZ() + deltacenter*axex);
-       dir1 = normp;
-       dir2.SetXYZ(axex);
-       param1    = radius;
-       param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
+        typeres = IntAna_Ellipse;
+        nbint = 1;
+        Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
+        deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
+        axex.Normalize();
+        pt1.SetXYZ(center.XYZ() + deltacenter*axex);
+        dir1 = normp;
+        dir2.SetXYZ(axex);
+        param1    = radius;
+        param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
       }
     }
   }
 //purpose  : Pln Sphere
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
-                                        const gp_Sphere& S)
+                                         const gp_Sphere& S)
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
-                                  ,const gp_Sphere& S) 
+                                   ,const gp_Sphere& S) 
 {
   
   done = Standard_False;
 //purpose  : Cylinder - Cylinder
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
-                                        const gp_Cylinder& Cyl2,
-                                        const Standard_Real Tol) 
+                                         const gp_Cylinder& Cyl2,
+                                         const Standard_Real Tol) 
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
-                    const gp_Cylinder& Cyl2,
-                    const Standard_Real Tol) 
+                     const gp_Cylinder& Cyl2,
+                     const Standard_Real Tol) 
 {
   done=Standard_True;
   //---------------------------- Parallel axes -------------------------
   if(A1A2.Parallel()) {
     if(DistA1A2<=Tol) {
       if(RmR<=Tol) {
-       typeres=IntAna_Same;
+        typeres=IntAna_Same;
       }
       else {
-       typeres=IntAna_Empty;
+        typeres=IntAna_Empty;
       }
     }
     else {  //-- DistA1A2 > Tol
       Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
       
       P2.SetCoord( P2t.X() - ProjP2OnDirCyl1*DirCyl.X()
-                 ,P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y()
-                 ,P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
+                  ,P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y()
+                  ,P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
       //-- 
       Standard_Real R1pR2=R1+R2;
       if(DistA1A2>(R1pR2+Tol)) { 
-       typeres=IntAna_Empty;
-       nbint=0;
+        typeres=IntAna_Empty;
+        nbint=0;
       }
       else if(DistA1A2>(R1pR2)) {
-       //-- 1 Tangent line -------------------------------------OK
-       typeres=IntAna_Line;
-
-       nbint=1;
-       dir1=DirCyl;
-       Standard_Real R1_R1pR2=R1/R1pR2;
-       pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X())
-                    ,P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y())
-                    ,P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
-       
+        //-- 1 Tangent line -------------------------------------OK
+        typeres=IntAna_Line;
+
+        nbint=1;
+        dir1=DirCyl;
+        Standard_Real R1_R1pR2=R1/R1pR2;
+        pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X())
+                     ,P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y())
+                     ,P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
+        
       }
       else if(DistA1A2>RmR) {
-       //-- 2 lines ---------------------------------------------OK
-       typeres=IntAna_Line;
-       nbint=2;
-       dir1=DirCyl;
-       gp_Vec P1P2(P1,P2);
-       gp_Dir DirA1A2=gp_Dir(P1P2);
-       gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
-       dir2=dir1;
-       Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);       
-
-//     Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
-       Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
-       Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
-       
-       if((Beta+Beta)<Tol) { 
-         nbint=1;
-         pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
-                      ,P1.Y() + Alpha*DirA1A2.Y()
-                      ,P1.Z() + Alpha*DirA1A2.Z());
-       }
-       else { 
-         pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X()
-                      ,P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y()
-                      ,P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z() );
-         
-         pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X()
-                      ,P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y()
-                      ,P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
-       }
+        //-- 2 lines ---------------------------------------------OK
+        typeres=IntAna_Line;
+        nbint=2;
+        dir1=DirCyl;
+        gp_Vec P1P2(P1,P2);
+        gp_Dir DirA1A2=gp_Dir(P1P2);
+        gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
+        dir2=dir1;
+        Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);       
+
+//         Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
+        Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
+        Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
+        
+        if((Beta+Beta)<Tol) { 
+          nbint=1;
+          pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
+                       ,P1.Y() + Alpha*DirA1A2.Y()
+                       ,P1.Z() + Alpha*DirA1A2.Z());
+        }
+        else { 
+          pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X()
+                       ,P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y()
+                       ,P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z() );
+          
+          pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X()
+                       ,P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y()
+                       ,P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
+        }
       }
       else if(DistA1A2>(RmR-Tol)) {
-       //-- 1 Tangent ------------------------------------------OK
-       typeres=IntAna_Line;
-       nbint=1;
-       dir1=DirCyl;
-       Standard_Real R1_RmR=R1/RmR;
+        //-- 1 Tangent ------------------------------------------OK
+        typeres=IntAna_Line;
+        nbint=1;
+        dir1=DirCyl;
+        Standard_Real R1_RmR=R1/RmR;
 
-       if(R1 < R2) R1_RmR = -R1_RmR;
+        if(R1 < R2) R1_RmR = -R1_RmR;
 
-       pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X())
-                    ,P1.Y() + R1_RmR * (P2.Y()-P1.Y())
-                    ,P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
+        pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X())
+                     ,P1.Y() + R1_RmR * (P2.Y()-P1.Y())
+                     ,P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
       }
       else {
-       nbint=0;
-       typeres=IntAna_Empty;
+        nbint=0;
+        typeres=IntAna_Empty;
       }
     }
   }
       A=Abs(Sin(0.5*A));
       
       if(A==0.0 || B==0.0) {
-       typeres=IntAna_Same;
-       return;
+        typeres=IntAna_Same;
+        return;
       }
       
       
       gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
       dir1 = gp_Dir(dircyl1.Added(dircyl2));
       dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
-       
+        
       param2   = Cyl1.Radius() / A;
       param1   = Cyl1.Radius() / B;
       param2bis= param1bis = Cyl1.Radius();
       if(param1 < param1bis) {
-       A=param1; param1=param1bis; param1bis=A;
+        A=param1; param1=param1bis; param1bis=A;
       }
       if(param2 < param2bis) {
-       A=param2; param2=param2bis; param2bis=A;
+        A=param2; param2=param2bis; param2bis=A;
       }
     }
     else {
       if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol) { 
-       typeres = IntAna_Point;
-       Standard_Real d,p1,p2;
-
-       gp_Dir D1 = Cyl1.Axis().Direction();
-       gp_Dir D2 = Cyl2.Axis().Direction();
-       A1A2.Distance(d,p1,p2);
-       gp_Pnt P = Cyl1.Axis().Location();
-       gp_Pnt P1(P.X() - p1*D1.X(),
-                 P.Y() - p1*D1.Y(),
-                 P.Z() - p1*D1.Z());
-       P = Cyl2.Axis().Location();
-       gp_Pnt P2(P.X() - p2*D2.X(),
-                 P.Y() - p2*D2.Y(),
-                 P.Z() - p2*D2.Z());
-       gp_Vec P1P2(P1,P2);
-       D1=gp_Dir(P1P2);
-       p1=Cyl1.Radius();
-       pt1.SetCoord(P1.X() + p1*D1.X(),
-                    P1.Y() + p1*D1.Y(),
-                    P1.Z() + p1*D1.Z());
-       nbint = 1;
+        typeres = IntAna_Point;
+        Standard_Real d,p1,p2;
+
+        gp_Dir D1 = Cyl1.Axis().Direction();
+        gp_Dir D2 = Cyl2.Axis().Direction();
+        A1A2.Distance(d,p1,p2);
+        gp_Pnt P = Cyl1.Axis().Location();
+        gp_Pnt P1(P.X() - p1*D1.X(),
+                  P.Y() - p1*D1.Y(),
+                  P.Z() - p1*D1.Z());
+        P = Cyl2.Axis().Location();
+        gp_Pnt P2(P.X() - p2*D2.X(),
+                  P.Y() - p2*D2.Y(),
+                  P.Z() - p2*D2.Z());
+        gp_Vec P1P2(P1,P2);
+        D1=gp_Dir(P1P2);
+        p1=Cyl1.Radius();
+        pt1.SetCoord(P1.X() + p1*D1.X(),
+                     P1.Y() + p1*D1.Y(),
+                     P1.Z() + p1*D1.Z());
+        nbint = 1;
       }
       else {
-       typeres=IntAna_NoGeometricSolution;
+        typeres=IntAna_NoGeometricSolution;
       }
     }
   }
 //purpose  : Cylinder - Cone
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
-                                        const gp_Cone& Con,
-                                        const Standard_Real Tol) 
+                                         const gp_Cone& Con,
+                                         const Standard_Real Tol) 
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
-                                  const gp_Cone& Con,
-                                  const Standard_Real ) 
+                                   const gp_Cone& Con,
+                                   const Standard_Real ) 
 {
   done=Standard_True;
   AxeOperator A1A2(Cyl.Axis(),Con.Axis());
     Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
     gp_Dir dir=Cyl.Position().Direction();
     pt1.SetCoord( Pt.X() + dist*dir.X()
-                ,Pt.Y() + dist*dir.Y()
-                ,Pt.Z() + dist*dir.Z());
+                 ,Pt.Y() + dist*dir.Y()
+                 ,Pt.Z() + dist*dir.Z());
     pt2.SetCoord( Pt.X() - dist*dir.X()
-                ,Pt.Y() - dist*dir.Y()
-                ,Pt.Z() - dist*dir.Z());
+                 ,Pt.Y() - dist*dir.Y()
+                 ,Pt.Z() - dist*dir.Z());
     dir1=dir2=dir;
     param1=param2=Cyl.Radius();
     nbint=2;
 //purpose  : Cylinder - Sphere
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
-                                        const gp_Sphere& Sph,
-                                        const Standard_Real Tol) 
+                                         const gp_Sphere& Sph,
+                                         const Standard_Real Tol) 
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
-                                  ,const gp_Sphere& Sph
-                                  ,const Standard_Real)
+                                   ,const gp_Sphere& Sph
+                                   ,const Standard_Real)
 {
   done=Standard_True;
   gp_Pnt Pt=Sph.Location();
       dir1 = dir2 = dir;
       typeres=IntAna_Circle;
       pt1.SetCoord( Pt.X() + dist*dir.X()
-                  ,Pt.Y() + dist*dir.Y()
-                  ,Pt.Z() + dist*dir.Z());
+                   ,Pt.Y() + dist*dir.Y()
+                   ,Pt.Z() + dist*dir.Z());
       nbint=1;
       param1 = Cyl.Radius();
       if(dist>RealEpsilon()) {
-       pt2.SetCoord( Pt.X() - dist*dir.X()
-                    ,Pt.Y() - dist*dir.Y()
-                    ,Pt.Z() - dist*dir.Z());
-       param2=Cyl.Radius();
-       nbint=2;
+        pt2.SetCoord( Pt.X() - dist*dir.X()
+                     ,Pt.Y() - dist*dir.Y()
+                     ,Pt.Z() - dist*dir.Z());
+        param2=Cyl.Radius();
+        nbint=2;
       }
     }
   }
 //purpose  : Cone - Cone
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
-                                        const gp_Cone& Con2,
-                                        const Standard_Real Tol) 
+                                         const gp_Cone& Con2,
+                                         const Standard_Real Tol) 
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
-                                  const gp_Cone& Con2,
-                                  const Standard_Real Tol) 
+                                   const gp_Cone& Con2,
+                                   const Standard_Real Tol) 
 {
   done=Standard_True;
   //
     
     if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) { 
       if (fabs(d) < TOL_APEX_CONF) {
-       typeres = IntAna_Point;
-       nbint = 1;
-       pt1 = P;
-       return;
+        typeres = IntAna_Point;
+        nbint = 1;
+        pt1 = P;
+        return;
       }
       x=(d*tg2)/(tg1+tg2);
       pt1.SetCoord( P.X() + x*D.X()
-                  ,P.Y() + x*D.Y()
-                  ,P.Z() + x*D.Z());
+                   ,P.Y() + x*D.Y()
+                   ,P.Z() + x*D.Z());
       param1=Abs(x*tg1);
 
       x=(d*tg2)/(tg2-tg1);
       pt2.SetCoord( P.X() + x*D.X()
-                  ,P.Y() + x*D.Y()
-                  ,P.Z() + x*D.Z());
+                   ,P.Y() + x*D.Y()
+                   ,P.Z() + x*D.Z());
       param2=Abs(x*tg1);
       dir1 = dir2 = D;
       nbint=2;
     }
     else {
       if (fabs(d) < TOL_APEX_CONF) { 
-       typeres=IntAna_Same;
+        typeres=IntAna_Same;
       }
       else {
-       typeres=IntAna_Circle;
-       nbint=1;
-       x=d*0.5;
-       pt1.SetCoord( P.X() + x*D.X()
-                    ,P.Y() + x*D.Y()
-                    ,P.Z() + x*D.Z());
-       param1 = Abs(x * tg1);
-       dir1 = D;
+        typeres=IntAna_Circle;
+        nbint=1;
+        x=d*0.5;
+        pt1.SetCoord( P.X() + x*D.X()
+                     ,P.Y() + x*D.Y()
+                     ,P.Z() + x*D.Z());
+        param1 = Abs(x * tg1);
+        dir1 = D;
       }
     }
   } //-- fin A1A2.Same
     Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
 
     gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
-                     O1O2n.Y()-O1O2_DA1*DA1.Y(),
-                     O1O2n.Z()-O1O2_DA1*DA1.Z());
+                      O1O2n.Y()-O1O2_DA1*DA1.Y(),
+                      O1O2n.Z()-O1O2_DA1*DA1.Z());
     gp_Dir DB1=gp_Dir(O1_Proj_A2);
 
     Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
     Standard_Real X1 = X2+yO1O2;
     
     gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
-             Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()), 
-             Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
+              Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()), 
+              Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
 
     gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
-                0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
-                0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
+                 0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
+                 0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
     gp_Vec P1MO1O2(P1,MO1O2);
     
     gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
     IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
       if(INTER_QUAD_PLN.IsDone()) {
       switch(INTER_QUAD_PLN.TypeInter()) {
-      case IntAna_Ellipse:     {
-       typeres=IntAna_Ellipse;
-       gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
-       pt1 = E.Location();
-       dir1 = E.Position().Direction();
-       dir2 = E.Position().XDirection();
-       param1 = E.MajorRadius();
-       param1bis = E.MinorRadius();
-       nbint = 1;
-       break; 
+      case IntAna_Ellipse:         {
+        typeres=IntAna_Ellipse;
+        gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
+        pt1 = E.Location();
+        dir1 = E.Position().Direction();
+        dir2 = E.Position().XDirection();
+        param1 = E.MajorRadius();
+        param1bis = E.MinorRadius();
+        nbint = 1;
+        break; 
       }
       case IntAna_Circle: {
-       typeres=IntAna_Circle;
-       gp_Circ C=INTER_QUAD_PLN.Circle(1);
-       pt1 = C.Location();
-       dir1 = C.Position().XDirection();
-       dir2 = C.Position().YDirection();
-       param1 = C.Radius();
-       nbint = 1;
-       break;
+        typeres=IntAna_Circle;
+        gp_Circ C=INTER_QUAD_PLN.Circle(1);
+        pt1 = C.Location();
+        dir1 = C.Position().XDirection();
+        dir2 = C.Position().YDirection();
+        param1 = C.Radius();
+        nbint = 1;
+        break;
       }
       case IntAna_Hyperbola: {
-       typeres=IntAna_Hyperbola;
-       gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
-       pt1 = pt2 = H.Location();
-       dir1 = H.Position().Direction();
-       dir2 = H.Position().XDirection();
-       param1 = param2 = H.MajorRadius();
-       param1bis = param2bis = H.MinorRadius();
-       nbint = 2;
-       break;
+        typeres=IntAna_Hyperbola;
+        gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
+        pt1 = pt2 = H.Location();
+        dir1 = H.Position().Direction();
+        dir2 = H.Position().XDirection();
+        param1 = param2 = H.MajorRadius();
+        param1bis = param2bis = H.MinorRadius();
+        nbint = 2;
+        break;
       }
       case IntAna_Line: {
-       typeres=IntAna_Line;
-       gp_Lin H=INTER_QUAD_PLN.Line(1);
-       pt1 = pt2 = H.Location();
-       dir1 = dir2 = H.Position().Direction();
-       param1 = param2 = 0.0;
-       param1bis = param2bis = 0.0;
-       nbint = 2;
-       break;
+        typeres=IntAna_Line;
+        gp_Lin H=INTER_QUAD_PLN.Line(1);
+        pt1 = pt2 = H.Location();
+        dir1 = dir2 = H.Position().Direction();
+        param1 = param2 = 0.0;
+        param1bis = param2bis = 0.0;
+        nbint = 2;
+        break;
       }
       default:
-       typeres=IntAna_NoGeometricSolution; 
+        typeres=IntAna_NoGeometricSolution; 
       }
     }
   }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
     aQApex1=Con1.Apex();
     aD3Ax1=aAx1.Direction(); 
     aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
-                 aQApex1.Y()+aD1*aD3Ax1.Y(),
-                 aQApex1.Z()+aD1*aD3Ax1.Z());
+                  aQApex1.Y()+aD1*aD3Ax1.Y(),
+                  aQApex1.Z()+aD1*aD3Ax1.Z());
     //
     aDx=aD3Ax1.Dot(aAx2.Direction());
     if (aDx<0.) {
     aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
     //
     aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
-                 aQApex1.Y()+aD2*aD3Ax2.Y(),
-                 aQApex1.Z()+aD2*aD3Ax2.Z());
+                  aQApex1.Y()+aD2*aD3Ax2.Y(),
+                  aQApex1.Z()+aD2*aD3Ax2.Z());
     //
     gp_Pln aPln1(aQA1, aD3Ax1);
     gp_Pln aPln2(aQA2, aD3Ax2);
 
       if(INTER_QUAD_PLN.IsDone()) {
       switch(INTER_QUAD_PLN.TypeInter()) {
-      case IntAna_Ellipse:     {
-       typeres=IntAna_Ellipse;
-       gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
-       pt1 = E.Location();
-       dir1 = E.Position().Direction();
-       dir2 = E.Position().XDirection();
-       param1 = E.MajorRadius();
-       param1bis = E.MinorRadius();
-       nbint = 1;
-       break; 
+      case IntAna_Ellipse:         {
+        typeres=IntAna_Ellipse;
+        gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
+        pt1 = E.Location();
+        dir1 = E.Position().Direction();
+        dir2 = E.Position().XDirection();
+        param1 = E.MajorRadius();
+        param1bis = E.MinorRadius();
+        nbint = 1;
+        break; 
       }
       case IntAna_Circle: {
-       typeres=IntAna_Circle;
-       gp_Circ C=INTER_QUAD_PLN.Circle(1);
-       pt1 = C.Location();
-       dir1 = C.Position().XDirection();
-       dir2 = C.Position().YDirection();
-       param1 = C.Radius();
-       nbint = 1;
-       break;
+        typeres=IntAna_Circle;
+        gp_Circ C=INTER_QUAD_PLN.Circle(1);
+        pt1 = C.Location();
+        dir1 = C.Position().XDirection();
+        dir2 = C.Position().YDirection();
+        param1 = C.Radius();
+        nbint = 1;
+        break;
       }
       case IntAna_Parabola: {
-       typeres=IntAna_Parabola;
-       gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
-       pt1 = Prb.Location();
-       dir1 = Prb.Position().Direction();
-       dir2 = Prb.Position().XDirection();
-       param1 = Prb.Focal();
-       nbint = 1;
-       break;
+        typeres=IntAna_Parabola;
+        gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
+        pt1 = Prb.Location();
+        dir1 = Prb.Position().Direction();
+        dir2 = Prb.Position().XDirection();
+        param1 = Prb.Focal();
+        nbint = 1;
+        break;
       }
       case IntAna_Hyperbola: {
-       typeres=IntAna_Hyperbola;
-       gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
-       pt1 = pt2 = H.Location();
-       dir1 = H.Position().Direction();
-       dir2 = H.Position().XDirection();
-       param1 = param2 = H.MajorRadius();
-       param1bis = param2bis = H.MinorRadius();
-       nbint = 2;
-       break;
+        typeres=IntAna_Hyperbola;
+        gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
+        pt1 = pt2 = H.Location();
+        dir1 = H.Position().Direction();
+        dir2 = H.Position().XDirection();
+        param1 = param2 = H.MajorRadius();
+        param1bis = param2bis = H.MinorRadius();
+        nbint = 2;
+        break;
       }
       default:
-       typeres=IntAna_NoGeometricSolution; 
+        typeres=IntAna_NoGeometricSolution; 
       }
     }
   }
 //purpose  : Sphere - Cone
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
-                                        const gp_Cone& Con,
-                                        const Standard_Real Tol) 
+                                         const gp_Cone& Con,
+                                         const Standard_Real Tol) 
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
-                                  const gp_Cone& Con,
-                                  const Standard_Real)
+                                   const gp_Cone& Con,
+                                   const Standard_Real)
 {
   
   //
     //--                tga = y / (x+dApexSphCenter)
     Standard_Real tgatga = tga * tga;
     math_DirectPolynomialRoots Eq( 1.0+tgatga
-                                 ,2.0*tgatga*dApexSphCenter
-                                 ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
+                                  ,2.0*tgatga*dApexSphCenter
+                                  ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
     if(Eq.IsDone()) {
       Standard_Integer nbsol=Eq.NbSolutions();
       if(nbsol==0) {
-       typeres=IntAna_Empty;
+        typeres=IntAna_Empty;
       }
       else { 
-       typeres=IntAna_Circle;
-       if(nbsol>=1) {
-         Standard_Real x = Eq.Value(1);
-         Standard_Real dApexSphCenterpx = dApexSphCenter+x;
-         nbint=1;
-         pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
-                      ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
-                      ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
-         param1 = tga * dApexSphCenterpx;
-         param1 = Abs(param1);
-         dir1 = ConDir;
-         if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
-           typeres=IntAna_PointAndCircle;
-           param1=0.0;
-         }
-       }
-       if(nbsol>=2) {
-         Standard_Real x=Eq.Value(2);
-         Standard_Real dApexSphCenterpx = dApexSphCenter+x;
-         nbint=2;
-         pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
-                      ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
-                      ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
-         param2 = tga * dApexSphCenterpx;
-         param2 = Abs(param2);
-         dir2=ConDir;
-         if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
-           typeres=IntAna_PointAndCircle;
-           param2=0.0;
-         }
-       }
+        typeres=IntAna_Circle;
+        if(nbsol>=1) {
+          Standard_Real x = Eq.Value(1);
+          Standard_Real dApexSphCenterpx = dApexSphCenter+x;
+          nbint=1;
+          pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
+                       ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
+                       ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
+          param1 = tga * dApexSphCenterpx;
+          param1 = Abs(param1);
+          dir1 = ConDir;
+          if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
+            typeres=IntAna_PointAndCircle;
+            param1=0.0;
+          }
+        }
+        if(nbsol>=2) {
+          Standard_Real x=Eq.Value(2);
+          Standard_Real dApexSphCenterpx = dApexSphCenter+x;
+          nbint=2;
+          pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
+                       ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
+                       ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
+          param2 = tga * dApexSphCenterpx;
+          param2 = Abs(param2);
+          dir2=ConDir;
+          if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
+            typeres=IntAna_PointAndCircle;
+            param2=0.0;
+          }
+        }
       }
     }
     else {
 //purpose  : Sphere - Sphere
 //=======================================================================
   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
-                                        ,const gp_Sphere& Sph2
-                                        ,const Standard_Real Tol) 
+                                         ,const gp_Sphere& Sph2
+                                         ,const Standard_Real Tol) 
 : done(Standard_False),
   nbint(0),
   typeres(IntAna_Empty),
   pt1(0,0,0),
   pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
   param1(0),
   param2(0),
+  param3(0),
+  param4(0),
   param1bis(0),
   param2bis(0),
   myCommonGen(Standard_False),
 //purpose  : 
 //=======================================================================
   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
-                                  const gp_Sphere& Sph2,
-                                  const Standard_Real Tol)   
+                                   const gp_Sphere& Sph2,
+                                   const Standard_Real Tol)   
 {
   done=Standard_True;
   gp_Pnt O1=Sph1.Location();
       Standard_Real t2;
       if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
       else         t2=(-R1+(dO1O2-R2))*0.5;
-       
+        
       pt1.SetCoord( O1.X() + t2*Dir.X()
-                  ,O1.Y() + t2*Dir.Y()
-                  ,O1.Z() + t2*Dir.Z());
+                   ,O1.Y() + t2*Dir.Y()
+                   ,O1.Z() + t2*Dir.Z());
     }
     else  {
       //-----------------------------------------------------------------
       //--                                            
       //-----------------------------------------------------------------
       if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
-       typeres=IntAna_Empty;
+        typeres=IntAna_Empty;
       }
       else {
-       //---------------------------------------------------------------
-       //--     
-       //--
-       //---------------------------------------------------------------
-       Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);       
-       Standard_Real Beta = R1*R1-Alpha*Alpha;
-       Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
-       
-       if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) { 
-         typeres = IntAna_Point;
-         Alpha = (R1 + (dO1O2 - R2)) * 0.5;
-       }
-       else { 
-         typeres = IntAna_Circle;
-         dir1 = Dir;
-         param1 = Beta;
-       }         
-       pt1.SetCoord( O1.X() + Alpha*Dir.X()
-                    ,O1.Y() + Alpha*Dir.Y()
-                    ,O1.Z() + Alpha*Dir.Z());
-       
-       nbint=1;
+        //---------------------------------------------------------------
+        //--     
+        //--
+        //---------------------------------------------------------------
+        Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);       
+        Standard_Real Beta = R1*R1-Alpha*Alpha;
+        Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
+        
+        if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) { 
+          typeres = IntAna_Point;
+          Alpha = (R1 + (dO1O2 - R2)) * 0.5;
+        }
+        else { 
+          typeres = IntAna_Circle;
+          dir1 = Dir;
+          param1 = Beta;
+        }          
+        pt1.SetCoord( O1.X() + Alpha*Dir.X()
+                     ,O1.Y() + Alpha*Dir.Y()
+                     ,O1.Z() + Alpha*Dir.Z());
+        
+        nbint=1;
       }
     }
   }
 }
+
+//=======================================================================
+//function : IntAna_QuadQuadGeo
+//purpose  : Plane - Torus
+//=======================================================================
+IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
+                                       const gp_Torus& Tor,
+                                       const Standard_Real Tol) 
+: done(Standard_False),
+  nbint(0),
+  typeres(IntAna_Empty),
+  pt1(0,0,0),
+  pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
+  param1(0),
+  param2(0),
+  param3(0),
+  param4(0),
+  param1bis(0),
+  param2bis(0),
+  myCommonGen(Standard_False),
+  myPChar(0,0,0)
+{
+  InitTolerances();
+  Perform(Pln,Tor,Tol);
+}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
+                                 const gp_Torus& Tor,
+                                 const Standard_Real Tol)
+{
+  done = Standard_True;
+  //
+  Standard_Real aRMin, aRMaj;
+  //
+  aRMin = Tor.MinorRadius();
+  aRMaj = Tor.MajorRadius();
+  if (aRMin >= aRMaj) {
+    typeres = IntAna_NoGeometricSolution; 
+    return;
+  }
+  //
+  const gp_Ax1 aPlnAx = Pln.Axis();
+  const gp_Ax1 aTorAx = Tor.Axis();
+  //
+  Standard_Boolean bParallel, bNormal;
+  //
+  bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
+  bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
+  if (!bNormal && !bParallel) {
+    typeres = IntAna_NoGeometricSolution; 
+    return;
+  }
+  //
+  Standard_Real aDist;
+  //
+  gp_Pnt aTorLoc = aTorAx.Location();
+  if (bParallel) {
+    Standard_Real aDt, X, Y, Z, A, B, C, D;
+    //
+    Pln.Coefficients(A,B,C,D);
+    aTorLoc.Coord(X,Y,Z);
+    aDist = A*X + B*Y + C*Z + D;
+    //
+    if ((Abs(aDist) - aRMin) > Tol) {
+      typeres=IntAna_Empty;
+      return;
+    }
+    //
+    typeres = IntAna_Circle;
+    //
+    pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
+    aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
+    param1 = aRMaj + aDt;
+    dir1 = aTorAx.Direction();
+    nbint = 1;
+    if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
+      pt2 = pt1;
+      param2 = aRMaj - aDt;
+      dir2 = dir1;
+      nbint = 2;
+    }
+  }
+  //
+  else {
+    aDist = Pln.Distance(aTorLoc);
+    if (aDist > myEPSILON_DISTANCE) {
+      typeres = IntAna_NoGeometricSolution; 
+      return;
+    }
+    //
+    typeres = IntAna_Circle;
+    param2 = param1 = aRMin;
+    dir2 = dir1 = aPlnAx.Direction();
+    nbint = 2;
+    //
+    gp_Dir aDir = aTorAx.Direction()^dir1;
+    pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
+    pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
+  }
+}
+
+//=======================================================================
+//function : IntAna_QuadQuadGeo
+//purpose  : Cylinder - Torus
+//=======================================================================
+IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
+                                       const gp_Torus& Tor,
+                                       const Standard_Real Tol) 
+: done(Standard_False),
+  nbint(0),
+  typeres(IntAna_Empty),
+  pt1(0,0,0),
+  pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
+  param1(0),
+  param2(0),
+  param3(0),
+  param4(0),
+  param1bis(0),
+  param2bis(0),
+  myCommonGen(Standard_False),
+  myPChar(0,0,0)
+{
+  InitTolerances();
+  Perform(Cyl,Tor,Tol);
+}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
+                                 const gp_Torus& Tor,
+                                 const Standard_Real Tol) 
+{
+  done = Standard_True;
+  //
+  Standard_Real aRMin, aRMaj;
+  //
+  aRMin = Tor.MinorRadius();
+  aRMaj = Tor.MajorRadius();
+  if (aRMin >= aRMaj) {
+    typeres = IntAna_NoGeometricSolution; 
+    return;
+  }
+  //
+  const gp_Ax1 aCylAx = Cyl.Axis();
+  const gp_Ax1 aTorAx = Tor.Axis();
+  //
+  const gp_Lin aLin(aTorAx);
+  const gp_Pnt aLocCyl = Cyl.Location();
+  //
+  if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
+      (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
+    typeres = IntAna_NoGeometricSolution; 
+    return;
+  }
+  //
+  Standard_Real aRCyl;
+  //
+  aRCyl = Cyl.Radius();
+  if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
+    typeres = IntAna_Empty;
+    return;
+  }
+  //
+  typeres = IntAna_Circle;
+  //
+  Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
+  gp_XYZ aTorLoc = aTorAx.Location().XYZ();
+  //
+  dir1 = aTorAx.Direction();
+  pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
+  param1 = aRCyl;
+  nbint = 1;
+  if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
+                       (aRCyl < (aRMaj + aRMin))) {
+    dir2 = dir1;
+    pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
+    param2 = param1;
+    nbint = 2;
+  }
+}
+
+//=======================================================================
+//function : IntAna_QuadQuadGeo
+//purpose  : Cone - Torus
+//=======================================================================
+IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
+                                       const gp_Torus& Tor,
+                                       const Standard_Real Tol) 
+: done(Standard_False),
+  nbint(0),
+  typeres(IntAna_Empty),
+  pt1(0,0,0),
+  pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
+  param1(0),
+  param2(0),
+  param3(0),
+  param4(0),
+  param1bis(0),
+  param2bis(0),
+  myCommonGen(Standard_False),
+  myPChar(0,0,0)
+{
+  InitTolerances();
+  Perform(Con,Tor,Tol);
+}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
+                                 const gp_Torus& Tor,
+                                 const Standard_Real Tol) 
+{
+  done = Standard_True;
+  //
+  Standard_Real aRMin, aRMaj;
+  //
+  aRMin = Tor.MinorRadius();
+  aRMaj = Tor.MajorRadius();
+  if (aRMin >= aRMaj) {
+    typeres = IntAna_NoGeometricSolution; 
+    return;
+  }
+  //
+  const gp_Ax1 aConAx = Con.Axis();
+  const gp_Ax1 aTorAx = Tor.Axis();
+  //
+  const gp_Lin aLin(aTorAx);
+  const gp_Pnt aConApex = Con.Apex();
+  //
+  if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
+      (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
+    typeres = IntAna_NoGeometricSolution; 
+    return;
+  }
+  //
+  Standard_Real anAngle, aDist, aParam[4];
+  Standard_Integer i;
+  gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
+  gp_Dir aDir[4];
+  //
+  anAngle = Con.SemiAngle();
+  aTorLoc = aTorAx.Location();
+  //
+  aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
+  gp_Dir aDN (gp_Vec(aTorLoc, aPN));
+  gp_Ax1 anAxCLRot(aConApex, aDN);
+  gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
+  gp_Dir aDL = aConL.Position().Direction();
+  gp_Dir aXDir = Tor.XAxis().Direction();
+  //
+  typeres = IntAna_Empty;
+  //
+  for (i = 0; i < 2; ++i) {
+    if (i) {
+      aXDir.Reverse();
+    }
+    aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
+    //
+    aDist = aConL.Distance(aPCT);
+    if (aDist > aRMin+Tol) {
+      continue;
+    }
+    //
+    typeres = IntAna_Circle;
+    //
+    gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
+    aDist = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
+    //
+    gp_Pnt aP;
+    gp_XYZ aDVal = aDist*aDL.XYZ();
+    aP.SetXYZ(aPh + aDVal);
+    aParam[nbint] = aLin.Distance(aP);
+    aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
+    aDir[nbint] = aTorAx.Direction();
+    ++nbint;
+    if ((aDist < aRMin) && (aDVal.Modulus() > Tol)) {
+      aP.SetXYZ(aPh - aDVal);
+      aParam[nbint] = aLin.Distance(aP);
+      aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
+      aDir[nbint] = aDir[nbint-1];
+      ++nbint;
+    }
+  }
+  //
+  for (i = 0; i < nbint; ++i) {
+    switch (i) {
+    case 0:{
+      pt1 = aPt[i];
+      param1 = aParam[i];
+      dir1 = aDir[i];
+      break;
+    }
+    case 1:{
+      pt2 = aPt[i];
+      param2 = aParam[i];
+      dir2 = aDir[i];
+      break;
+    }
+    case 2:{
+      pt3 = aPt[i];
+      param3 = aParam[i];
+      dir3 = aDir[i];
+      break;
+    }
+    case 3:{
+      pt4 = aPt[i];
+      param4 = aParam[i];
+      dir4 = aDir[i];
+      break;
+    }
+    default:
+      break;
+    }
+  }
+}
+
+//=======================================================================
+//function : IntAna_QuadQuadGeo
+//purpose  : Sphere - Torus
+//=======================================================================
+IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
+                                       const gp_Torus& Tor,
+                                       const Standard_Real Tol) 
+: done(Standard_False),
+  nbint(0),
+  typeres(IntAna_Empty),
+  pt1(0,0,0),
+  pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
+  param1(0),
+  param2(0),
+  param3(0),
+  param4(0),
+  param1bis(0),
+  param2bis(0),
+  myCommonGen(Standard_False),
+  myPChar(0,0,0)
+{
+  InitTolerances();
+  Perform(Sph,Tor,Tol);
+}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
+                                 const gp_Torus& Tor,
+                                 const Standard_Real Tol) 
+{
+  done = Standard_True;
+  //
+  Standard_Real aRMin, aRMaj;
+  //
+  aRMin = Tor.MinorRadius();
+  aRMaj = Tor.MajorRadius();
+  if (aRMin >= aRMaj) {
+    typeres = IntAna_NoGeometricSolution; 
+    return;
+  }
+  //
+  const gp_Ax1 aTorAx = Tor.Axis();
+  const gp_Lin aLin(aTorAx);
+  const gp_Pnt aSphLoc = Sph.Location();
+  //
+  if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
+    typeres = IntAna_NoGeometricSolution;
+    return;
+  }
+  //
+  Standard_Real aRSph, aDist;
+  gp_Pnt aTorLoc;
+  //
+  gp_Dir aXDir = Tor.XAxis().Direction();
+  aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
+  aRSph = Sph.Radius();
+  //
+  gp_Vec aVec12(aTorLoc, aSphLoc);
+  aDist = aVec12.Magnitude();
+  if (((aDist - Tol) > (aRMin + aRSph)) || 
+      ((aDist + Tol) < Abs(aRMin - aRSph))) {
+    typeres = IntAna_Empty;
+    return;
+  }
+  //
+  typeres = IntAna_Circle;
+  //
+  Standard_Real anAlpha, aBeta;
+  //
+  anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
+  aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
+  //
+  gp_Dir aDir12(aVec12);
+  gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
+  gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
+  //
+  gp_Pnt aP;
+  gp_XYZ aDVal = aBeta*aDC.XYZ();
+  aP.SetXYZ(aPh + aDVal);
+  param1 = aLin.Distance(aP);
+  pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
+  dir1 = aTorAx.Direction();
+  nbint = 1;
+  if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) && 
+      (aDVal.Modulus() > Tol)) {
+    aP.SetXYZ(aPh - aDVal);
+    param2 = aLin.Distance(aP);
+    pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
+    dir2 = dir1;
+    nbint = 2;
+  }
+}
+
+//=======================================================================
+//function : IntAna_QuadQuadGeo
+//purpose  : Torus - Torus
+//=======================================================================
+IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
+                                       const gp_Torus& Tor2,
+                                       const Standard_Real Tol) 
+: done(Standard_False),
+  nbint(0),
+  typeres(IntAna_Empty),
+  pt1(0,0,0),
+  pt2(0,0,0),
+  pt3(0,0,0),
+  pt4(0,0,0),
+  param1(0),
+  param2(0),
+  param3(0),
+  param4(0),
+  param1bis(0),
+  param2bis(0),
+  myCommonGen(Standard_False),
+  myPChar(0,0,0)
+{
+  InitTolerances();
+  Perform(Tor1,Tor2,Tol);
+}
+//=======================================================================
+//function : Perform
+//purpose  : 
+//=======================================================================
+void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
+                                 const gp_Torus& Tor2,
+                                 const Standard_Real Tol) 
+{
+  done = Standard_True;
+  //
+  Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
+  //
+  aRMin1 = Tor1.MinorRadius();
+  aRMaj1 = Tor1.MajorRadius();
+  aRMin2 = Tor2.MinorRadius();
+  aRMaj2 = Tor2.MajorRadius();
+  if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
+    typeres = IntAna_NoGeometricSolution;
+    return;
+  }
+  //
+  const gp_Ax1 anAx1 = Tor1.Axis();
+  const gp_Ax1 anAx2 = Tor2.Axis();
+  //
+  gp_Lin aL1(anAx1);
+  if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
+      (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
+    typeres = IntAna_NoGeometricSolution;
+    return;
+  }
+  //
+  gp_Pnt aLoc1, aLoc2;
+  //
+  aLoc1 = anAx1.Location();
+  aLoc2 = anAx2.Location();
+  //
+  if (aLoc1.IsEqual(aLoc2, Tol) &&
+      (Abs(aRMin1 - aRMin2) <= Tol) && 
+      (Abs(aRMaj1 - aRMaj2) <= Tol)) {
+    typeres = IntAna_Same;
+    return;
+  }
+  //
+  Standard_Real aDist;
+  gp_Pnt aP1, aP2;
+  //
+  gp_Dir aXDir1 = Tor1.XAxis().Direction();
+  aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
+  aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
+  //
+  gp_Vec aV12(aP1, aP2);
+  aDist = aV12.Magnitude();
+  if (((aDist - Tol) > (aRMin1 + aRMin2)) || 
+      ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
+    typeres = IntAna_Empty;
+    return;
+  }
+  //
+  typeres = IntAna_Circle;
+  //
+  Standard_Real anAlpha, aBeta;
+  //
+  anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
+  aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
+  //
+  gp_Dir aDir12(aV12);
+  gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
+  gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
+  //
+  gp_Pnt aP;
+  gp_XYZ aDVal = aBeta*aDC.XYZ();
+  aP.SetXYZ(aPh + aDVal);
+  param1 = aL1.Distance(aP);
+  pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
+  dir1 = anAx1.Direction();
+  nbint = 1;
+  if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) && 
+      aDVal.Modulus() > Tol) {
+    aP.SetXYZ(aPh - aDVal);
+    param2 = aL1.Distance(aP);
+    pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
+    dir2 = dir1;
+    nbint = 2;
+  }
+}
+
 //=======================================================================
 //function : Point
 //purpose  : Returns a Point
   else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
     Standard_DomainError::Raise();
     }
-  if(n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));   }
-  else {     return(gp_Circ(DirToAx2(pt2,dir2),param2));   }
+  if      (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
+  else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
+  else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
+  else           { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
 }
 
 //=======================================================================
     Standard_OutOfRange::Raise();
   }
   return(gp_Parab(gp_Ax2( pt1
-                        ,dir1
-                        ,dir2)
-                 ,param1));
+                         ,dir1
+                         ,dir2)
+                  ,param1));
 }
 //=======================================================================
 //function : Hyperbola
     }
   if(n==1) {
     return(gp_Hypr(gp_Ax2( pt1
-                         ,dir1
-                         ,dir2)
-                  ,param1,param1bis));
+                          ,dir1
+                          ,dir2)
+                   ,param1,param1bis));
   }
   else {
     return(gp_Hypr(gp_Ax2( pt2
-                         ,dir1
-                         ,dir2.Reversed())
-                  ,param2,param2bis));
+                          ,dir1
+                          ,dir2.Reversed())
+                   ,param2,param2bis));
   }
 }
 //=======================================================================
       m=(aC[k]>0.);
       aNum=(m)? aC[k] : -aC[k];
       if (aNum>aR1 && aNum<aR2) {
-       if (m) {
-         aC[k]=1.;
-       }         
-       else {
-         aC[k]=-1.;
-       }
-       //
-       aC[(k+1)%3]=0.;
-       aC[(k+2)%3]=0.;
-       break;
+        if (m) {
+          aC[k]=1.;
+        }          
+        else {
+          aC[k]=-1.;
+        }
+        //
+        aC[(k+1)%3]=0.;
+        aC[(k+2)%3]=0.;
+        break;
       }
     }
     aDir.SetCoord(aC[0], aC[1], aC[2]);
 
 IntPatch_ImpImpIntersection_3.gxx
 IntPatch_ImpImpIntersection_4.gxx
 IntPatch_ImpImpIntersection_5.gxx
+IntPatch_ImpImpIntersection_6.gxx
 
   case GeomAbs_Sphere:
     pu1=M_PI+M_PI;
     break;
+  case GeomAbs_Torus:
+    pu1=pv1=M_PI+M_PI;
+    break;
   default:
     break;
   }
   case GeomAbs_Sphere:
     pu2=M_PI+M_PI;
     break;
+  case GeomAbs_Torus:
+    pu2=pv2=M_PI+M_PI;
+    break;
   default:
     break;
   }
 
 #include <IntPatch_ImpImpIntersection_3.gxx>
 #include <IntPatch_ImpImpIntersection_4.gxx>
 #include <IntPatch_ImpImpIntersection_5.gxx>
+#include <IntPatch_ImpImpIntersection_6.gxx>
 
                                Standard_Boolean&,
                                IntPatch_SequenceOfLine&,
                                IntPatch_SequenceOfPoint&);
+
+//torus
+static Standard_Boolean IntPTo(const IntSurf_Quadric&,
+                               const IntSurf_Quadric&,
+                               const Standard_Real,
+                               const Standard_Boolean,
+                               Standard_Boolean&,
+                               IntPatch_SequenceOfLine&);
+
+static Standard_Boolean IntCyTo(const IntSurf_Quadric&,
+                                const IntSurf_Quadric&,
+                                const Standard_Real,
+                                const Standard_Boolean,
+                                Standard_Boolean&,
+                                IntPatch_SequenceOfLine&);
+
+static Standard_Boolean IntCoTo(const IntSurf_Quadric&,
+                                const IntSurf_Quadric&,
+                                const Standard_Real,
+                                const Standard_Boolean,
+                                Standard_Boolean&,
+                                IntPatch_SequenceOfLine&);
+
+static Standard_Boolean IntSpTo(const IntSurf_Quadric&,
+                                const IntSurf_Quadric&,
+                                const Standard_Real,
+                                const Standard_Boolean,
+                                Standard_Boolean&,
+                                IntPatch_SequenceOfLine&);
+
+static Standard_Boolean IntToTo(const IntSurf_Quadric&,
+                                const IntSurf_Quadric&,
+                                const Standard_Real,
+                                Standard_Boolean&,
+                                Standard_Boolean&,
+                                IntPatch_SequenceOfLine&);
 
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
+static 
+  Standard_Integer SetQuad(const Handle(Adaptor3d_HSurface)& theS,
+                           GeomAbs_SurfaceType& theTS,
+                           IntSurf_Quadric& theQuad);
+
 //=======================================================================
 //function : IntPatch_ImpImpIntersection
 //purpose  : 
 //=======================================================================
 IntPatch_ImpImpIntersection::IntPatch_ImpImpIntersection
        (const Handle(Adaptor3d_HSurface)&  S1,
-       const Handle(Adaptor3d_TopolTool)& D1,
+        const Handle(Adaptor3d_TopolTool)& D1,
         const Handle(Adaptor3d_HSurface)&  S2,
-       const Handle(Adaptor3d_TopolTool)& D2,
-       const Standard_Real TolArc,
-       const Standard_Real TolTang)
+        const Handle(Adaptor3d_TopolTool)& D2,
+        const Standard_Real TolArc,
+        const Standard_Real TolTang)
 {
   Perform(S1,D1,S2,D2,TolArc,TolTang);
 }
 //purpose  : 
 //=======================================================================
 void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
-                                         const Handle(Adaptor3d_TopolTool)& D1,
-                                         const Handle(Adaptor3d_HSurface)&  S2,
-                                         const Handle(Adaptor3d_TopolTool)& D2,
-                                         const Standard_Real TolArc,
-                                         const Standard_Real TolTang) {
+                                          const Handle(Adaptor3d_TopolTool)& D1,
+                                          const Handle(Adaptor3d_HSurface)&  S2,
+                                          const Handle(Adaptor3d_TopolTool)& D2,
+                                          const Standard_Real TolArc,
+                                          const Standard_Real TolTang) {
   done = Standard_False;
   spnt.Clear();
   slin.Clear();
   IntPatch_SequenceOfPathPointOfTheSOnBounds pnt1,pnt2;
   //
   // On commence par intersecter les supports des surfaces
-  IntSurf_Quadric quad1;
-  IntSurf_Quadric quad2;
+  IntSurf_Quadric quad1, quad2;
   IntPatch_ArcFunction AFunc;
   const Standard_Real Tolang = 1.e-8;
-  GeomAbs_SurfaceType typs1 = S1->GetType();
-  GeomAbs_SurfaceType typs2 = S2->GetType();
+  GeomAbs_SurfaceType typs1, typs2;
+  Standard_Boolean bEmpty = Standard_False;
   //
-  switch (typs1) {
-
-  case GeomAbs_Plane :
-    {
-      quad1.SetValue(S1->Plane());
-
-      switch (typs2) {
-
-      case GeomAbs_Plane:
-       {
-         quad2.SetValue(S2->Plane());
-         if (!IntPP(quad1,quad2,Tolang,TolTang,SameSurf,slin)) {
-           return;
-         }
-       }
-       break;
-       
-       
-      case GeomAbs_Cylinder:
-       {
-         quad2.SetValue(S2->Cylinder());
-          Standard_Real VMin, VMax, H;
-          //
-          VMin = S1->FirstVParameter();
-          VMax = S1->LastVParameter();
-          H = (Precision::IsNegativeInfinite(VMin) || 
-               Precision::IsPositiveInfinite(VMax)) ? 0 : (VMax - VMin);
-         if (!IntPCy(quad1,quad2,Tolang,TolTang,Standard_False,empt,slin,H)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Sphere:
-       {
-         quad2.SetValue(S2->Sphere());
-         //modified by NIZNHY-PKV Tue Sep 20 09:03:06 2011f
-         if (!IntPSp(quad1,quad2,Tolang,TolTang,Standard_False,empt,slin,spnt)) {
-         //if (!IntPSp(quad1,quad2,TolTang,Standard_False,empt,slin,spnt)) {
-           //modified by NIZNHY-PKV Tue Sep 20 09:03:10 2011t
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-       
-      case GeomAbs_Cone:
-       {
-         quad2.SetValue(S2->Cone());
-         if (!IntPCo(quad1,quad2,Tolang,TolTang,Standard_False,
-                     empt,multpoint,slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-      default: 
-       {
-         Standard_ConstructionError::Raise();
-         break;
-       }
+  const Standard_Integer iT1 = SetQuad(S1, typs1, quad1);
+  const Standard_Integer iT2 = SetQuad(S2, typs2, quad2);
+  //
+  if (!iT1 || !iT2) {
+    Standard_ConstructionError::Raise();
+    return;
+  }
+  //
+  const Standard_Boolean bReverse = iT1 > iT2;
+  const Standard_Integer iTT = iT1*10 + iT2;
+  //
+  switch (iTT) {
+    case 11: { // Plane/Plane
+      if (!IntPP(quad1, quad2, Tolang, TolTang, SameSurf, slin)) {
+        return;
       }
+      break;
     }
-    break;
-
-  case GeomAbs_Cylinder:
-    {
-      quad1.SetValue(S1->Cylinder());
-      switch (typs2){
-
-      case GeomAbs_Plane:
-       {
-         quad2.SetValue(S2->Plane());
-          Standard_Real VMin, VMax, H;
-          //
-          VMin = S1->FirstVParameter();
-          VMax = S1->LastVParameter();
-          H = (Precision::IsNegativeInfinite(VMin) || 
-               Precision::IsPositiveInfinite(VMax)) ? 0 : (VMax - VMin);
-         if (!IntPCy(quad1,quad2,Tolang,TolTang,Standard_True,empt,slin,H)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Cylinder:
-       {
-         quad2.SetValue(S2->Cylinder());
-         if (!IntCyCy(quad1,quad2,TolTang,empt,SameSurf,multpoint,slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Sphere:
-       {
-         quad2.SetValue(S2->Sphere());
-         if (!IntCySp(quad1,quad2,TolTang,Standard_False,empt,multpoint,
-                      slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Cone:
-       {
-         quad2.SetValue(S2->Cone());
-         if (!IntCyCo(quad1,quad2,TolTang,Standard_False,empt,multpoint,
-                      slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-      default: 
-       {
-         Standard_ConstructionError::Raise();
-         break;
-       }
+    //
+    case 12:
+    case 21: { // Plane/Cylinder
+      Standard_Real VMin, VMax, H;
+      //
+      const Handle(Adaptor3d_HSurface)& aSCyl = bReverse ? S2 : S1;
+      VMin = aSCyl->FirstVParameter();
+      VMax = aSCyl->LastVParameter();
+      H = (Precision::IsNegativeInfinite(VMin) || 
+           Precision::IsPositiveInfinite(VMax)) ? 0 : (VMax - VMin);
+      //
+      if (!IntPCy(quad1, quad2, Tolang, TolTang, bReverse, empt, slin, H)) {
+        return;
       }
-
+      bEmpty = empt;
+      break;
     }
-    break;
-
-  case GeomAbs_Sphere:
-    {
-      quad1.SetValue(S1->Sphere());
-
-      switch (typs2){
-
-      case GeomAbs_Plane:
-       {
-         quad2.SetValue(S2->Plane());
-         //modified by NIZNHY-PKV Tue Sep 20 09:03:35 2011f
-         if (!IntPSp(quad1,quad2,Tolang,TolTang,Standard_True,empt,slin,spnt)) {
-         //if (!IntPSp(quad1,quad2,TolTang,Standard_True,empt,slin,spnt)) {
-           //modified by NIZNHY-PKV Tue Sep 20 09:03:38 2011t
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Cylinder:
-       {
-         quad2.SetValue(S2->Cylinder());
-         if (!IntCySp(quad1,quad2,TolTang,Standard_True,empt,multpoint,
-                      slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Sphere:
-       {
-         quad2.SetValue(S2->Sphere());
-         if (!IntSpSp(quad1,quad2,TolTang,empt,SameSurf,slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Cone:
-       {
-         quad2.SetValue(S2->Cone());
-         if (!IntCoSp(quad1,quad2,TolTang,Standard_True,empt,multpoint,
-                      slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-      default: 
-       {
-         Standard_ConstructionError::Raise();
-         break;
-       }
+    //
+    case 13:
+    case 31: { // Plane/Cone
+      if (!IntPCo(quad1, quad2, Tolang, TolTang, bReverse, empt, multpoint, slin, spnt)) {
+        return;
       }
-
+      bEmpty = empt;
+      break;
     }
-    break;
-
-  case GeomAbs_Cone:
-    {
-      quad1.SetValue(S1->Cone());
-
-      switch (typs2){
-
-      case GeomAbs_Plane:
-       {
-         quad2.SetValue(S2->Plane());
-         if (!IntPCo(quad1,quad2,Tolang,TolTang,Standard_True,
-                     empt,multpoint,slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Cylinder:
-       {
-         quad2.SetValue(S2->Cylinder());
-         if (!IntCyCo(quad1,quad2,TolTang,Standard_True,empt,multpoint,
-                      slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Sphere:
-       {
-         quad2.SetValue(S2->Sphere());
-         if (!IntCoSp(quad1,quad2,TolTang,Standard_False,empt,multpoint,
-                      slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      case GeomAbs_Cone:
-       {
-         quad2.SetValue(S2->Cone());
-         if (!IntCoCo(quad1,quad2,TolTang,empt,SameSurf,multpoint,
-                      slin,spnt)) {
-           return;
-         }
-         if (empt) {
-           done = Standard_True;
-           return;
-         }
-       }
-       break;
-
-      default:
-       {
-         Standard_ConstructionError::Raise();
-         break;
-       }
+    //
+    case 14:
+    case 41: { // Plane/Sphere
+      if (!IntPSp(quad1, quad2, Tolang, TolTang, bReverse, empt, slin, spnt)) {
+        return;
       }
-
+      bEmpty = empt;
+      break;
     }
-    break;
-    default:
-    {
+    //
+    case 15:
+    case 51: { // Plane/Torus
+      if (!IntPTo(quad1, quad2, TolTang, bReverse, empt, slin)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 22: { // Cylinder/Cylinder
+      if (!IntCyCy(quad1, quad2, TolTang, empt, SameSurf, multpoint, slin, spnt)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 23:
+    case 32: { // Cylinder/Cone
+      if (!IntCyCo(quad1, quad2, TolTang, bReverse, empt, multpoint, slin, spnt)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 24:
+    case 42: { // Cylinder/Sphere
+      if (!IntCySp(quad1, quad2, TolTang, bReverse, empt, multpoint, slin, spnt)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 25:
+    case 52: { // Cylinder/Torus
+      if (!IntCyTo(quad1, quad2, TolTang, bReverse, empt, slin)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 33: { // Cone/Cone
+      if (!IntCoCo(quad1, quad2, TolTang, empt, SameSurf, multpoint, slin, spnt)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 34:
+    case 43: { // Cone/Sphere
+      if (!IntCoSp(quad1, quad2, TolTang, bReverse, empt, multpoint, slin, spnt)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 35:
+    case 53: { // Cone/Torus
+      if (!IntCoTo(quad1, quad2, TolTang, bReverse, empt, slin)) {
+        return;
+      }
+      break;
+    }
+    //
+    case 44: { // Sphere/Sphere
+      if (!IntSpSp(quad1, quad2, TolTang, empt, SameSurf, slin, spnt)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 45:
+    case 54: { // Sphere/Torus
+      if (!IntSpTo(quad1, quad2, TolTang, bReverse, empt, slin)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    case 55: { // Torus/Torus
+      if (!IntToTo(quad1, quad2, TolTang, SameSurf, empt, slin)) {
+        return;
+      }
+      bEmpty = empt;
+      break;
+    }
+    //
+    default: {
       Standard_ConstructionError::Raise();
       break;
     }
-  } //switch (typs1) {
+  }
+  //
+  if (bEmpty) {
+    done = Standard_True;
+    return;
+  }
   //
   if (!SameSurf) {
     AFunc.SetQuadric(quad2);
       Ptreference = ElSLib::Value(0.,10.,S1->Cone());
     }
       break;
+    case GeomAbs_Torus:      {
+      Ptreference = ElSLib::Value(0.,0.,S1->Torus());
+    }
+      break;
     default: 
       break;
     }
     empt = Standard_False;
     // C est la qu il faut commencer a bosser...
     PutPointsOnLine(S1,S2,pnt1, slin, Standard_True, D1, quad1,quad2,
-                   multpoint,TolArc);
+                    multpoint,TolArc);
     
     PutPointsOnLine(S1,S2,pnt2, slin, Standard_False,D2, quad2,quad1,
-                   multpoint,TolArc);
+                    multpoint,TolArc);
 
     if (edg1.Length() != 0) {
       ProcessSegments(edg1,slin,quad1,quad2,Standard_True,TolArc);
       aState2=D2->Classify(aP2D, TolArc);
       //
       if(aState1!=TopAbs_OUT && aState2!=TopAbs_OUT) {
-       aSIP.Append(aIP);
+        aSIP.Append(aIP);
       }
     }
     //
       const Handle(IntPatch_GLine)& glin = *((Handle(IntPatch_GLine)*)&slin.Value(i));
       nbv = glin->NbVertex();
       if(glin->NbVertex() == 0) { 
-       gp_Circ Circ = glin->Circle();
-       P=ElCLib::Value(0.0,Circ);
-       quad1.Parameters(P,u1,v1);
-       quad2.Parameters(P,u2,v2);
-       point.SetValue(P,TolArc,Standard_False);
-       point.SetParameters(u1,v1,u2,v2);
-       point.SetParameter(0.0);
-       glin->AddVertex(point);
-
-       P=ElCLib::Value(0.0,Circ);
-       quad1.Parameters(P,u1,v1);
-       quad2.Parameters(P,u2,v2);
-       point.SetValue(P,TolArc,Standard_False);
-       point.SetParameters(u1,v1,u2,v2);
-       point.SetParameter(M_PI+M_PI);
-       glin->AddVertex(point);
+        gp_Circ Circ = glin->Circle();
+        P=ElCLib::Value(0.0,Circ);
+        quad1.Parameters(P,u1,v1);
+        quad2.Parameters(P,u2,v2);
+        point.SetValue(P,TolArc,Standard_False);
+        point.SetParameters(u1,v1,u2,v2);
+        point.SetParameter(0.0);
+        glin->AddVertex(point);
+
+        P=ElCLib::Value(0.0,Circ);
+        quad1.Parameters(P,u1,v1);
+        quad2.Parameters(P,u2,v2);
+        point.SetValue(P,TolArc,Standard_False);
+        point.SetParameters(u1,v1,u2,v2);
+        point.SetParameter(M_PI+M_PI);
+        glin->AddVertex(point);
       }
     }
     
       const Handle(IntPatch_GLine)& glin = *((Handle(IntPatch_GLine)*)&slin.Value(i));
       nbv = glin->NbVertex();
       if(glin->NbVertex() == 0) { 
-       gp_Elips Elips = glin->Ellipse();
-       P=ElCLib::Value(0.0,Elips);
-       quad1.Parameters(P,u1,v1);
-       quad2.Parameters(P,u2,v2);
-       point.SetValue(P,TolArc,Standard_False);
-       point.SetParameters(u1,v1,u2,v2);
-       point.SetParameter(0.0);
-       glin->AddVertex(point);
-
-       P=ElCLib::Value(0.0,Elips);
-       quad1.Parameters(P,u1,v1);
-       quad2.Parameters(P,u2,v2);
-       point.SetValue(P,TolArc,Standard_False);
-       point.SetParameters(u1,v1,u2,v2);
-       point.SetParameter(M_PI+M_PI);
-       glin->AddVertex(point);
+        gp_Elips Elips = glin->Ellipse();
+        P=ElCLib::Value(0.0,Elips);
+        quad1.Parameters(P,u1,v1);
+        quad2.Parameters(P,u2,v2);
+        point.SetValue(P,TolArc,Standard_False);
+        point.SetParameters(u1,v1,u2,v2);
+        point.SetParameter(0.0);
+        glin->AddVertex(point);
+
+        P=ElCLib::Value(0.0,Elips);
+        quad1.Parameters(P,u1,v1);
+        quad2.Parameters(P,u2,v2);
+        point.SetValue(P,TolArc,Standard_False);
+        point.SetParameters(u1,v1,u2,v2);
+        point.SetParameter(M_PI+M_PI);
+        glin->AddVertex(point);
       }
     }
   }
   done = Standard_True;
 }
 
+//=======================================================================
+//function : SetQuad
+//purpose  : 
+//=======================================================================
+Standard_Integer SetQuad(const Handle(Adaptor3d_HSurface)& theS,
+                         GeomAbs_SurfaceType& theTS,
+                         IntSurf_Quadric& theQuad)
+{
+  theTS = theS->GetType();
+  Standard_Integer iRet = 0;
+  switch (theTS) {
+  case GeomAbs_Plane:
+    theQuad.SetValue(theS->Plane());
+    iRet = 1;
+    break;
+  case GeomAbs_Cylinder:
+    theQuad.SetValue(theS->Cylinder());
+    iRet = 2;
+    break;
+  case GeomAbs_Cone:
+    theQuad.SetValue(theS->Cone());
+    iRet = 3;
+    break;
+  case GeomAbs_Sphere:
+    theQuad.SetValue(theS->Sphere());
+    iRet = 4;
+    break;
+  case GeomAbs_Torus:
+    theQuad.SetValue(theS->Torus());
+    iRet = 5;
+    break;
+  default:
+    break;
+  }
+  //
+  return iRet;
+}
+
 
 //modified by NIZNHY-PKV Thu Sep 15 11:09:12 2011
 static 
   void SeamPosition(const gp_Pnt& aPLoc, 
-                   const gp_Ax3& aPos,
-                   gp_Ax2& aSeamPos);
+                    const gp_Ax3& aPos,
+                    gp_Ax2& aSeamPos);
 static
   void AdjustToSeam (const gp_Cylinder& aQuad,
-                    gp_Circ& aCirc);
+                     gp_Circ& aCirc);
 static
   void AdjustToSeam (const gp_Sphere& aQuad,
-                    gp_Circ& aCirc,
-                    const Standard_Real aTolAng);
+                     gp_Circ& aCirc,
+                     const Standard_Real aTolAng);
 static
   void AdjustToSeam (const gp_Cone& aQuad,
-                    gp_Circ& aCirc);
+                     gp_Circ& aCirc);
+static
+  void AdjustToSeam (const gp_Torus& aQuad,
+                     gp_Circ& aCirc);
 //modified by NIZNHY-PKV Thu Sep 15 11:09:13 2011
 
 //=======================================================================
 // Traitement du cas Plan/Plan
 //=======================================================================
 Standard_Boolean IntPP (const IntSurf_Quadric& Quad1,
-                       const IntSurf_Quadric& Quad2,
-                       const Standard_Real Tolang,
-                       const Standard_Real TolTang,
-                       Standard_Boolean& Same,
-                       IntPatch_SequenceOfLine& slin)
+                        const IntSurf_Quadric& Quad2,
+                        const Standard_Real Tolang,
+                        const Standard_Real TolTang,
+                        Standard_Boolean& Same,
+                        IntPatch_SequenceOfLine& slin)
 
 {
   IntSurf_TypeTrans trans1,trans2;
 // Traitement du cas Plan/Cylindre et reciproquement
 //=======================================================================
 Standard_Boolean IntPCy (const IntSurf_Quadric& Quad1,
-                        const IntSurf_Quadric& Quad2,
-                        const Standard_Real Tolang,
-                        const Standard_Real TolTang,
-                        const Standard_Boolean Reversed,
-                        Standard_Boolean& Empty,
-                        IntPatch_SequenceOfLine& slin,
+                         const IntSurf_Quadric& Quad2,
+                         const Standard_Real Tolang,
+                         const Standard_Real TolTang,
+                         const Standard_Boolean Reversed,
+                         Standard_Boolean& Empty,
+                         IntPatch_SequenceOfLine& slin,
                          const Standard_Real H)
 
 {
   Empty = Standard_False;
 
   switch (typint) {
-           
+            
     case IntAna_Empty :    {
       Empty = Standard_True;
     }
       gp_Lin linsol = inter.Line(1);
       gp_Pnt orig(linsol.Location());
       if (NbSol == 1) {                 // ligne de tangence
-       gp_Vec TestCurvature(orig,Cy.Location());
-       gp_Vec Normp,Normcyl;
-       if (!Reversed) {
-         Normp = Quad1.Normale(orig);
-         Normcyl = Quad2.Normale(orig);
-       }
-       else {
-         Normp = Quad2.Normale(orig);
-         Normcyl = Quad1.Normale(orig);
-       }
-       
-       IntSurf_Situation situcyl;
-       IntSurf_Situation situp;
+        gp_Vec TestCurvature(orig,Cy.Location());
+        gp_Vec Normp,Normcyl;
+        if (!Reversed) {
+          Normp = Quad1.Normale(orig);
+          Normcyl = Quad2.Normale(orig);
+        }
+        else {
+          Normp = Quad2.Normale(orig);
+          Normcyl = Quad1.Normale(orig);
+        }
+        
+        IntSurf_Situation situcyl;
+        IntSurf_Situation situp;
 
-       if (Normp.Dot(TestCurvature) > 0.) {
-         situcyl = IntSurf_Outside;
-         if (Normp.Dot(Normcyl) > 0.) {
-           situp = IntSurf_Inside;
-         }
-         else {
-           situp = IntSurf_Outside;
-         }
-       }
-       else {
-         situcyl = IntSurf_Inside;
-         if (Normp.Dot(Normcyl) > 0.) {
-           situp = IntSurf_Outside;
-         }
-         else {
-           situp = IntSurf_Inside;
-         }
-       }
-       Handle(IntPatch_GLine) glig;
-       if (!Reversed) {
-         glig = new IntPatch_GLine(linsol, Standard_True, situp, situcyl);
-       }
-       else {
-         glig = new IntPatch_GLine(linsol, Standard_True, situcyl, situp);
-       }
-       slin.Append(glig);
+        if (Normp.Dot(TestCurvature) > 0.) {
+          situcyl = IntSurf_Outside;
+          if (Normp.Dot(Normcyl) > 0.) {
+            situp = IntSurf_Inside;
+          }
+          else {
+            situp = IntSurf_Outside;
+          }
+        }
+        else {
+          situcyl = IntSurf_Inside;
+          if (Normp.Dot(Normcyl) > 0.) {
+            situp = IntSurf_Outside;
+          }
+          else {
+            situp = IntSurf_Inside;
+          }
+        }
+        Handle(IntPatch_GLine) glig;
+        if (!Reversed) {
+          glig = new IntPatch_GLine(linsol, Standard_True, situp, situcyl);
+        }
+        else {
+          glig = new IntPatch_GLine(linsol, Standard_True, situcyl, situp);
+        }
+        slin.Append(glig);
       }
       else {      
-       // on a 2 droites. Il faut determiner les transitions
-       // de chacune.
-       
-       if (linsol.Direction().DotCross(Quad2.Normale(orig),
-                                       Quad1.Normale(orig)) >0.) {
-         trans1 = IntSurf_Out;
-         trans2 = IntSurf_In;
-       }
-       else {
-         trans1 = IntSurf_In;
-         trans2 = IntSurf_Out;
-       }
-       Handle(IntPatch_GLine) glig = 
-         new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
-       slin.Append(glig);
-       
-       linsol = inter.Line(2);
-       orig = linsol.Location();
+        // on a 2 droites. Il faut determiner les transitions
+        // de chacune.
+        
+        if (linsol.Direction().DotCross(Quad2.Normale(orig),
+                                        Quad1.Normale(orig)) >0.) {
+          trans1 = IntSurf_Out;
+          trans2 = IntSurf_In;
+        }
+        else {
+          trans1 = IntSurf_In;
+          trans2 = IntSurf_Out;
+        }
+        Handle(IntPatch_GLine) glig = 
+          new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
+        slin.Append(glig);
+        
+        linsol = inter.Line(2);
+        orig = linsol.Location();
 
-       if (linsol.Direction().DotCross(Quad2.Normale(orig),
-                                       Quad1.Normale(orig)) >0.) {
-         trans1 = IntSurf_Out;
-         trans2 = IntSurf_In;
-       }
-       else {
-         trans1 = IntSurf_In;
-         trans2 = IntSurf_Out;
-       }
-       glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
-       slin.Append(glig);
+        if (linsol.Direction().DotCross(Quad2.Normale(orig),
+                                        Quad1.Normale(orig)) >0.) {
+          trans1 = IntSurf_Out;
+          trans2 = IntSurf_In;
+        }
+        else {
+          trans1 = IntSurf_In;
+          trans2 = IntSurf_Out;
+        }
+        glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
+        slin.Append(glig);
       }
     }
       break;
       ElCLib::D1(0.,cirsol,ptref,Tgt);
       
       if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
-       trans1 = IntSurf_Out;
-       trans2 = IntSurf_In;
+        trans1 = IntSurf_Out;
+        trans2 = IntSurf_In;
       }
       else {
-       trans1 = IntSurf_In;
-       trans2 = IntSurf_Out;
+        trans1 = IntSurf_In;
+        trans2 = IntSurf_Out;
       }
       Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
       slin.Append(glig);
       ElCLib::D1(0.,elipsol,ptref,Tgt);
       
       if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
-       trans1 = IntSurf_Out;
-       trans2 = IntSurf_In;
+        trans1 = IntSurf_Out;
+        trans2 = IntSurf_In;
       }
       else {
-       trans1 = IntSurf_In;
-       trans2 = IntSurf_Out;
+        trans1 = IntSurf_In;
+        trans2 = IntSurf_Out;
       }
       Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
       slin.Append(glig);
 // Traitement du cas Plan/Sphere et reciproquement
 //=======================================================================
 Standard_Boolean IntPSp (const IntSurf_Quadric& Quad1,
-                        const IntSurf_Quadric& Quad2,
-                        //modified by NIZNHY-PKV Tue Sep 20 08:59:36 2011f
-                        const Standard_Real Tolang,
-                        //modified by NIZNHY-PKV Tue Sep 20 08:59:39 2011t
-                        const Standard_Real TolTang,
-                        const Standard_Boolean Reversed,
-                        Standard_Boolean& Empty,
-                        IntPatch_SequenceOfLine& slin,
-                        IntPatch_SequenceOfPoint& spnt)
+                         const IntSurf_Quadric& Quad2,
+                         //modified by NIZNHY-PKV Tue Sep 20 08:59:36 2011f
+                         const Standard_Real Tolang,
+                         //modified by NIZNHY-PKV Tue Sep 20 08:59:39 2011t
+                         const Standard_Real TolTang,
+                         const Standard_Boolean Reversed,
+                         Standard_Boolean& Empty,
+                         IntPatch_SequenceOfLine& slin,
+                         IntPatch_SequenceOfPoint& spnt)
 
 
 {
       ElCLib::D1(0.,cirsol,ptref,Tgt);
 
       if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) >0.) {
-       trans1 = IntSurf_Out;
-       trans2 = IntSurf_In;
+        trans1 = IntSurf_Out;
+        trans2 = IntSurf_In;
       }
       else {
-       trans1 = IntSurf_In;
-       trans2 = IntSurf_Out;
+        trans1 = IntSurf_In;
+        trans2 = IntSurf_Out;
       }
       Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
       slin.Append(glig);
 // Traitement du cas Plan/Cone et reciproquement
 //=======================================================================
 Standard_Boolean IntPCo (const IntSurf_Quadric& Quad1,
-                        const IntSurf_Quadric& Quad2,
-                        const Standard_Real Tolang,
-                        const Standard_Real TolTang,
-                        const Standard_Boolean Reversed,
-                        Standard_Boolean& Empty,
-                        Standard_Boolean& Multpoint,
-                        IntPatch_SequenceOfLine& slin,
-                        IntPatch_SequenceOfPoint& spnt)
+                         const IntSurf_Quadric& Quad2,
+                         const Standard_Real Tolang,
+                         const Standard_Real TolTang,
+                         const Standard_Boolean Reversed,
+                         Standard_Boolean& Empty,
+                         Standard_Boolean& Multpoint,
+                         IntPatch_SequenceOfLine& slin,
+                         IntPatch_SequenceOfPoint& spnt)
 
 
 {
     case IntAna_Line:    {
       gp_Lin linsol = inter.Line(1);
       if (linsol.Direction().Dot(Co.Axis().Direction()) <0.) {
-       linsol.SetDirection(linsol.Direction().Reversed());
+        linsol.SetDirection(linsol.Direction().Reversed());
       }
       Standard_Real para = ElCLib::Parameter(linsol, apex);
       gp_Pnt ptbid (ElCLib::Value(para+5.,linsol));
       Quad2.Parameters(apex,U2,V2);
       
       if (NbSol == 1) {                 // ligne de tangence
-       IntPatch_Point ptsol;
-       ptsol.SetValue(apex,TolTang,Standard_False);
-       ptsol.SetParameters(U1,V1,U2,V2);
-       ptsol.SetParameter(para);
-       gp_Pnt ptbid2(apex.XYZ() + 5.*Co.Axis().Direction().XYZ());
-       gp_Vec TestCurvature(ptbid,ptbid2);
-       gp_Vec Normp,Normco;
-       if (!Reversed) {
-         Normp = Quad1.Normale(ptbid);
-         Normco = Quad2.Normale(ptbid);
-       }
-       else {
-         Normp = Quad2.Normale(ptbid);
-         Normco = Quad1.Normale(ptbid);
-       }
-       IntSurf_Situation situco,situco_otherside;
-       IntSurf_Situation situp,situp_otherside;
-       
-       if (Normp.Dot(TestCurvature) > 0.) {
-         situco           = IntSurf_Outside;
-         situco_otherside = IntSurf_Inside;
-         if (Normp.Dot(Normco) > 0.) {
-           situp           = IntSurf_Inside;
-           situp_otherside = IntSurf_Outside;
-         }
-         else {
-           situp           = IntSurf_Outside;
-           situp_otherside = IntSurf_Inside;
-         }
-       }
-       else {
-         situco           = IntSurf_Inside;
-         situco_otherside = IntSurf_Outside;
-         if (Normp.Dot(Normco) > 0.) {
-           situp           = IntSurf_Outside;
-           situp_otherside = IntSurf_Inside;
-         }
-         else {
-           situp           = IntSurf_Inside;
-           situp_otherside = IntSurf_Outside;
-         }
-       }
-       //----------------------------------------------------------
-       //--              Apex ---> Cone.Direction
-       //--
-       Handle(IntPatch_GLine) glig;
-       if (!Reversed) {
-         glig = new IntPatch_GLine(linsol, Standard_True, situp, situco);
-       }
-       else {
-         glig = new IntPatch_GLine(linsol, Standard_True, situco, situp);
-       }
-       glig->AddVertex(ptsol);
-       glig->SetFirstPoint(1);
-       slin.Append(glig);
-       //----------------------------------------------------------
-       //--   -Cone.Direction <------- Apex
-       //--
-       linsol.SetDirection(linsol.Direction().Reversed());
-       if (!Reversed) {
-         glig = new IntPatch_GLine(linsol, Standard_True, situp_otherside, situco_otherside);
-       }
-       else {
-         glig = new IntPatch_GLine(linsol, Standard_True, situco_otherside, situp_otherside);
-       }
-       glig->AddVertex(ptsol);
-       glig->SetFirstPoint(1);
-       slin.Append(glig);
+        IntPatch_Point ptsol;
+        ptsol.SetValue(apex,TolTang,Standard_False);
+        ptsol.SetParameters(U1,V1,U2,V2);
+        ptsol.SetParameter(para);
+        gp_Pnt ptbid2(apex.XYZ() + 5.*Co.Axis().Direction().XYZ());
+        gp_Vec TestCurvature(ptbid,ptbid2);
+        gp_Vec Normp,Normco;
+        if (!Reversed) {
+          Normp = Quad1.Normale(ptbid);
+          Normco = Quad2.Normale(ptbid);
+        }
+        else {
+          Normp = Quad2.Normale(ptbid);
+          Normco = Quad1.Normale(ptbid);
+        }
+        IntSurf_Situation situco,situco_otherside;
+        IntSurf_Situation situp,situp_otherside;
+        
+        if (Normp.Dot(TestCurvature) > 0.) {
+          situco           = IntSurf_Outside;
+          situco_otherside = IntSurf_Inside;
+          if (Normp.Dot(Normco) > 0.) {
+            situp           = IntSurf_Inside;
+            situp_otherside = IntSurf_Outside;
+          }
+          else {
+            situp           = IntSurf_Outside;
+            situp_otherside = IntSurf_Inside;
+          }
+        }
+        else {
+          situco           = IntSurf_Inside;
+          situco_otherside = IntSurf_Outside;
+          if (Normp.Dot(Normco) > 0.) {
+            situp           = IntSurf_Outside;
+            situp_otherside = IntSurf_Inside;
+          }
+          else {
+            situp           = IntSurf_Inside;
+            situp_otherside = IntSurf_Outside;
+          }
+        }
+        //----------------------------------------------------------
+        //--              Apex ---> Cone.Direction
+        //--
+        Handle(IntPatch_GLine) glig;
+        if (!Reversed) {
+          glig = new IntPatch_GLine(linsol, Standard_True, situp, situco);
+        }
+        else {
+          glig = new IntPatch_GLine(linsol, Standard_True, situco, situp);
+        }
+        glig->AddVertex(ptsol);
+        glig->SetFirstPoint(1);
+        slin.Append(glig);
+        //----------------------------------------------------------
+        //--   -Cone.Direction <------- Apex
+        //--
+        linsol.SetDirection(linsol.Direction().Reversed());
+        if (!Reversed) {
+          glig = new IntPatch_GLine(linsol, Standard_True, situp_otherside, situco_otherside);
+        }
+        else {
+          glig = new IntPatch_GLine(linsol, Standard_True, situco_otherside, situp_otherside);
+        }
+        glig->AddVertex(ptsol);
+        glig->SetFirstPoint(1);
+        slin.Append(glig);
       }
       else {      
-       // on a 2 droites. Il faut determiner les transitions
-       // de chacune. On oriente chaque ligne dans le sens
-       // de l axe du cone. Les transitions de chaque ligne seront
-       // inverses l une de l autre => on ne fait le calcul que sur
-       // la premiere.
-       if (linsol.Direction().DotCross
-           (Quad2.Normale(ptbid),Quad1.Normale(ptbid)) >0.) {
-         trans1 = IntSurf_Out;
-         trans2 = IntSurf_In;
-       }
-       else {
-         trans1 = IntSurf_In;
-         trans2 = IntSurf_Out;
-       }
+        // on a 2 droites. Il faut determiner les transitions
+        // de chacune. On oriente chaque ligne dans le sens
+        // de l axe du cone. Les transitions de chaque ligne seront
+        // inverses l une de l autre => on ne fait le calcul que sur
+        // la premiere.
+        if (linsol.Direction().DotCross
+            (Quad2.Normale(ptbid),Quad1.Normale(ptbid)) >0.) {
+          trans1 = IntSurf_Out;
+          trans2 = IntSurf_In;
+        }
+        else {
+          trans1 = IntSurf_In;
+          trans2 = IntSurf_Out;
+        }
 
-       Multpoint = Standard_True;
-       //------------------------------------------- Ligne 1 -------
-       IntPatch_Point ptsol;
-       ptsol.SetValue(apex,TolTang,Standard_False);
-       ptsol.SetParameters(U1,V1,U2,V2);
-       ptsol.SetParameter(para);
-       ptsol.SetMultiple(Standard_True);
-       Handle(IntPatch_GLine) glig;
-       glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
-       glig->AddVertex(ptsol);
-       glig->SetFirstPoint(1);
-       slin.Append(glig);
-       //-----------------------------------------------------------
-       //-- Other Side : Les transitions restent les memes
-       //--    linsol -> -linsol   et Quad1(2).N -> -Quad1(2).N
-       //-- 
-       linsol.SetDirection(linsol.Direction().Reversed());
-       glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
-       para = ElCLib::Parameter(linsol, apex);
-       ptsol.SetParameter(para);
-       glig->AddVertex(ptsol);
-       glig->SetFirstPoint(1);
-       slin.Append(glig);
-       
-       //------------------------------------------- Ligne 2 -------
-       linsol = inter.Line(2);
-       if (linsol.Direction().Dot(Co.Axis().Direction()) <0.) {
-         linsol.SetDirection(linsol.Direction().Reversed());
-       }
-       para = ElCLib::Parameter(linsol, apex);
-       ptbid  = ElCLib::Value(para+5.,linsol);
-       if (linsol.Direction().DotCross
-           (Quad2.Normale(ptbid),Quad1.Normale(ptbid)) >0.) {
-         trans1 = IntSurf_Out;
-         trans2 = IntSurf_In;
-       }
-       else {
-         trans1 = IntSurf_In;
-         trans2 = IntSurf_Out;
-       }
-       ptsol.SetParameter(para);
-       glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
-       para = ElCLib::Parameter(linsol, apex);
-       ptsol.SetParameter(para);
-       glig->AddVertex(ptsol);
-       glig->SetFirstPoint(1);
-       slin.Append(glig);
-       //-----------------------------------------------------------
-       //-- Other Side : Les transitions restent les memes
-       //--    linsol -> -linsol   et Quad1(2).N -> -Quad1(2).N
-       //-- 
-       linsol.SetDirection(linsol.Direction().Reversed());
-       glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
-       para = ElCLib::Parameter(linsol, apex);
-       ptsol.SetParameter(para);
-       glig->AddVertex(ptsol);
-       glig->SetFirstPoint(1);
-       slin.Append(glig);
+        Multpoint = Standard_True;
+        //------------------------------------------- Ligne 1 -------
+        IntPatch_Point ptsol;
+        ptsol.SetValue(apex,TolTang,Standard_False);
+        ptsol.SetParameters(U1,V1,U2,V2);
+        ptsol.SetParameter(para);
+        ptsol.SetMultiple(Standard_True);
+        Handle(IntPatch_GLine) glig;
+        glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
+        glig->AddVertex(ptsol);
+        glig->SetFirstPoint(1);
+        slin.Append(glig);
+        //-----------------------------------------------------------
+        //-- Other Side : Les transitions restent les memes
+        //--    linsol -> -linsol   et Quad1(2).N -> -Quad1(2).N
+        //-- 
+        linsol.SetDirection(linsol.Direction().Reversed());
+        glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
+        para = ElCLib::Parameter(linsol, apex);
+        ptsol.SetParameter(para);
+        glig->AddVertex(ptsol);
+        glig->SetFirstPoint(1);
+        slin.Append(glig);
+        
+        //------------------------------------------- Ligne 2 -------
+        linsol = inter.Line(2);
+        if (linsol.Direction().Dot(Co.Axis().Direction()) <0.) {
+          linsol.SetDirection(linsol.Direction().Reversed());
+        }
+        para = ElCLib::Parameter(linsol, apex);
+        ptbid  = ElCLib::Value(para+5.,linsol);
+        if (linsol.Direction().DotCross
+            (Quad2.Normale(ptbid),Quad1.Normale(ptbid)) >0.) {
+          trans1 = IntSurf_Out;
+          trans2 = IntSurf_In;
+        }
+        else {
+          trans1 = IntSurf_In;
+          trans2 = IntSurf_Out;
+        }
+        ptsol.SetParameter(para);
+        glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
+        para = ElCLib::Parameter(linsol, apex);
+        ptsol.SetParameter(para);
+        glig->AddVertex(ptsol);
+        glig->SetFirstPoint(1);
+        slin.Append(glig);
+        //-----------------------------------------------------------
+        //-- Other Side : Les transitions restent les memes
+        //--    linsol -> -linsol   et Quad1(2).N -> -Quad1(2).N
+        //-- 
+        linsol.SetDirection(linsol.Direction().Reversed());
+        glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
+        para = ElCLib::Parameter(linsol, apex);
+        ptsol.SetParameter(para);
+        glig->AddVertex(ptsol);
+        glig->SetFirstPoint(1);
+        slin.Append(glig);
       }
     }
     break;
       ElCLib::D1(0.,cirsol,ptref,Tgt);
       
       if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) >0.) {
-       trans1 = IntSurf_Out;
-       trans2 = IntSurf_In;
+        trans1 = IntSurf_Out;
+        trans2 = IntSurf_In;
       }
       else {
-       trans1 = IntSurf_In;
-       trans2 = IntSurf_Out;
+        trans1 = IntSurf_In;
+        trans2 = IntSurf_Out;
       }
       Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
       slin.Append(glig);
       ElCLib::D1(0.,elipsol,ptref,Tgt);
       
       if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) >0.) {
-       trans1 = IntSurf_Out;
-       trans2 = IntSurf_In;
+        trans1 = IntSurf_Out;
+        trans2 = IntSurf_In;
       }
       else {
-       trans1 = IntSurf_In;
-       trans2 = IntSurf_Out;
+        trans1 = IntSurf_In;
+        trans2 = IntSurf_Out;
       }
       Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
       slin.Append(glig);
       
       gp_Vec Tgtorig(parabsol.YAxis().Direction());
       Standard_Real ptran = Tgtorig.DotCross(Quad2.Normale(parabsol.Location()),
-                                            Quad1.Normale(parabsol.Location()));
+                                             Quad1.Normale(parabsol.Location()));
       if (ptran >0.00000001) {
-       trans1 = IntSurf_Out;
-       trans2 = IntSurf_In;
+        trans1 = IntSurf_Out;
+        trans2 = IntSurf_In;
       }
       else if (ptran <-0.00000001) {
-       trans1 = IntSurf_In;
-       trans2 = IntSurf_Out;
+        trans1 = IntSurf_In;
+        trans2 = IntSurf_Out;
       }
       else { 
-       trans1=trans2=IntSurf_Undecided; 
+        trans1=trans2=IntSurf_Undecided; 
       }
       Handle(IntPatch_GLine) glig = new IntPatch_GLine(parabsol,Standard_False,trans1,trans2);
       slin.Append(glig);
       gp_Vec Tgttop;
       
       for(Standard_Integer i=1; i<=2; i++) { 
-       gp_Hypr hyprsol = inter.Hyperbola(i);
-       tophypr = ElCLib::Value(hyprsol.MajorRadius(), 
-                               hyprsol.XAxis());
-       Tgttop = hyprsol.YAxis().Direction();
-       Standard_Real qwe = Tgttop.DotCross(Quad2.Normale(tophypr),
-                                           Quad1.Normale(tophypr));
-       
-       if (qwe>0.00000001) { 
-         trans1 = IntSurf_Out;
-         trans2 = IntSurf_In;
-       }
-       else if (qwe<-0.00000001){
-         trans1 = IntSurf_In;
-         trans2 = IntSurf_Out;
-       }
-       else { 
-         trans1=trans2=IntSurf_Undecided;
-       }
-       Handle(IntPatch_GLine) glig = new IntPatch_GLine(hyprsol,Standard_False,trans1,trans2);
-       slin.Append(glig);
+        gp_Hypr hyprsol = inter.Hyperbola(i);
+        tophypr = ElCLib::Value(hyprsol.MajorRadius(), 
+                                hyprsol.XAxis());
+        Tgttop = hyprsol.YAxis().Direction();
+        Standard_Real qwe = Tgttop.DotCross(Quad2.Normale(tophypr),
+                                            Quad1.Normale(tophypr));
+        
+        if (qwe>0.00000001) { 
+          trans1 = IntSurf_Out;
+          trans2 = IntSurf_In;
+        }
+        else if (qwe<-0.00000001){
+          trans1 = IntSurf_In;
+          trans2 = IntSurf_Out;
+        }
+        else { 
+          trans1=trans2=IntSurf_Undecided;
+        }
+        Handle(IntPatch_GLine) glig = new IntPatch_GLine(hyprsol,Standard_False,trans1,trans2);
+        slin.Append(glig);
       }
     }
       break;
   }
   return Standard_True;
 }
+//=======================================================================
+//function : IntPTo
+//purpose  : 
+//=======================================================================
+Standard_Boolean IntPTo(const IntSurf_Quadric& theQuad1,
+                        const IntSurf_Quadric& theQuad2,
+                        const Standard_Real theTolTang,
+                        const Standard_Boolean bReversed,
+                        Standard_Boolean& bEmpty,
+                        IntPatch_SequenceOfLine& theSeqLin)
+{
+  const gp_Pln aPln = bReversed ? theQuad2.Plane() : theQuad1.Plane();
+  const gp_Torus aTorus = bReversed ? theQuad1.Torus() : theQuad2.Torus();
+  //
+  IntAna_QuadQuadGeo inter(aPln, aTorus, theTolTang);
+  Standard_Boolean bRet = inter.IsDone();
+  //
+  if (!bRet) {
+    return bRet;
+  }
+  //
+  IntAna_ResultType typint = inter.TypeInter();
+  Standard_Integer NbSol = inter.NbSolutions();
+  bEmpty = Standard_False;
+  //
+  switch (typint) {
+  case IntAna_Empty :
+    bEmpty = Standard_True;
+    break;
+  //
+  case IntAna_Circle : {
+    Standard_Integer i;
+    IntSurf_TypeTrans trans1, trans2;
+    gp_Pnt ptref;
+    gp_Vec Tgt;
+    //
+    for (i = 1; i <= NbSol; ++i) {
+      gp_Circ aC = inter.Circle(i);
+      if (!aPln.Axis().IsNormal(aTorus.Axis(), Precision::Angular())) {
+        AdjustToSeam(aTorus, aC);
+      }
+      ElCLib::D1(0., aC, ptref, Tgt);
+      //
+      if (Tgt.DotCross(theQuad2.Normale(ptref),theQuad1.Normale(ptref)) > 0.0) {
+        trans1 = IntSurf_Out;
+        trans2 = IntSurf_In;
+      }
+      else {
+        trans1 = IntSurf_In;
+        trans2 = IntSurf_Out;
+      }
+      //
+      Handle(IntPatch_GLine) glig = 
+        new IntPatch_GLine(aC, Standard_False, trans1, trans2);
+      theSeqLin.Append(glig);
+    }
+  }
+    break;
+  //
+  case IntAna_NoGeometricSolution:
+  default:
+    bRet = Standard_False;
+    break;
+  }
+  //
+  return bRet;
+}
 //
 //modified by NIZNHY-PKV Thu Sep 15 10:53:39 2011f
 //=======================================================================
 //purpose  : 
 //=======================================================================
 void AdjustToSeam (const gp_Cone& aQuad,
-                  gp_Circ& aCirc)
+                   gp_Circ& aCirc)
 {
    gp_Ax2 aAx2;
    //
 //purpose  : 
 //=======================================================================
 void AdjustToSeam (const gp_Sphere& aQuad,
-                  gp_Circ& aCirc,
-                  const Standard_Real aTolAng)
+                   gp_Circ& aCirc,
+                   const Standard_Real aTolAng)
 {
    gp_Ax2 aAx2;
    //
 //purpose  : 
 //=======================================================================
 void AdjustToSeam (const gp_Cylinder& aQuad,
-                  gp_Circ& aCirc)
+                   gp_Circ& aCirc)
 {
    gp_Ax2 aAx2;
    //
    aCirc.SetPosition(aAx2);
 } 
 //=======================================================================
+//function : AdjustToSeam
+//purpose  : 
+//=======================================================================
+void AdjustToSeam (const gp_Torus& aQuad,
+                   gp_Circ& aCirc)
+{
+  gp_Ax2 aAx2;
+  //
+  const gp_Pnt& aPLoc=aCirc.Location();
+  const gp_Ax3& aAx3=aQuad.Position();
+  SeamPosition(aPLoc, aAx3, aAx2);
+  aCirc.SetPosition(aAx2);
+} 
+//=======================================================================
 //function : SeamPosition
 //purpose  : 
 //=======================================================================
 void SeamPosition(const gp_Pnt& aPLoc, 
-                 const gp_Ax3& aPos, 
-                 gp_Ax2& aSeamPos)
+                  const gp_Ax3& aPos, 
+                  gp_Ax2& aSeamPos)
 {
   const gp_Dir& aDZ=aPos.Direction();
   const gp_Dir& aDX=aPos.XDirection();
   gp_Ax2 aAx2(aPLoc, aDZ, aDX);
   aSeamPos=aAx2;
 }
-                           
+    
 //modified by NIZNHY-PKV Thu Sep 15 10:53:41 2011t
 
--- /dev/null
+// Created on: 1992-05-07
+// Created by: Jacques GOUSSARD
+// Copyright (c) 1992-1999 Matra Datavision
+// Copyright (c) 1999-2012 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and / or modify it
+// under the terms of the GNU Lesser General Public version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+
+static 
+  Standard_Boolean TreatResultTorus(const IntSurf_Quadric& theQuad1,
+                                    const IntSurf_Quadric& theQuad2,
+                                    const IntAna_QuadQuadGeo& anInt,
+                                    Standard_Boolean& bEmpty,
+                                    IntPatch_SequenceOfLine& theSeqLin);
+
+//=======================================================================
+//function : IntCyTo
+//purpose  : 
+//=======================================================================
+Standard_Boolean IntCyTo(const IntSurf_Quadric& theQuad1,
+                         const IntSurf_Quadric& theQuad2,
+                         const Standard_Real theTolTang,
+                         const Standard_Boolean bReversed,
+                         Standard_Boolean& bEmpty,
+                         IntPatch_SequenceOfLine& theSeqLin)
+{
+  const gp_Cylinder aCyl = bReversed ? theQuad2.Cylinder() : theQuad1.Cylinder();
+  const gp_Torus aTorus = bReversed ? theQuad1.Torus() : theQuad2.Torus();
+  //
+  IntAna_QuadQuadGeo anInt(aCyl, aTorus, theTolTang);
+  Standard_Boolean bRet = 
+    TreatResultTorus(theQuad1, theQuad2, anInt, bEmpty, theSeqLin);
+  //
+  return bRet;
+}
+
+//=======================================================================
+//function : IntCoTo
+//purpose  : 
+//=======================================================================
+Standard_Boolean IntCoTo(const IntSurf_Quadric& theQuad1,
+                         const IntSurf_Quadric& theQuad2,
+                         const Standard_Real theTolTang,
+                         const Standard_Boolean bReversed,
+                         Standard_Boolean& bEmpty,
+                         IntPatch_SequenceOfLine& theSeqLin)
+{
+  const gp_Cone aCone = bReversed ? theQuad2.Cone() : theQuad1.Cone();
+  const gp_Torus aTorus = bReversed ? theQuad1.Torus() : theQuad2.Torus();
+  //
+  IntAna_QuadQuadGeo anInt(aCone, aTorus, theTolTang);
+  Standard_Boolean bRet = 
+    TreatResultTorus(theQuad1, theQuad2, anInt, bEmpty, theSeqLin);
+  //
+  return bRet;
+}
+
+//=======================================================================
+//function : IntSpTo
+//purpose  : 
+//=======================================================================
+Standard_Boolean IntSpTo(const IntSurf_Quadric& theQuad1,
+                         const IntSurf_Quadric& theQuad2,
+                         const Standard_Real theTolTang,
+                         const Standard_Boolean bReversed,
+                         Standard_Boolean& bEmpty,
+                         IntPatch_SequenceOfLine& theSeqLin)
+{
+  const gp_Sphere aSphere = bReversed ? theQuad2.Sphere() : theQuad1.Sphere();
+  const gp_Torus aTorus = bReversed ? theQuad1.Torus() : theQuad2.Torus();
+  //
+  IntAna_QuadQuadGeo anInt(aSphere, aTorus, theTolTang);
+  Standard_Boolean bRet = 
+    TreatResultTorus(theQuad1, theQuad2, anInt, bEmpty, theSeqLin);
+  //
+  return bRet;
+}
+
+//=======================================================================
+//function : IntToTo
+//purpose  : 
+//=======================================================================
+Standard_Boolean IntToTo(const IntSurf_Quadric& theQuad1,
+                         const IntSurf_Quadric& theQuad2,
+                         const Standard_Real theTolTang,
+                         Standard_Boolean& bSameSurf,
+                         Standard_Boolean& bEmpty,
+                         IntPatch_SequenceOfLine& theSeqLin)
+{
+  const gp_Torus aTorus1 = theQuad1.Torus();
+  const gp_Torus aTorus2 = theQuad2.Torus();
+  //
+  IntAna_QuadQuadGeo anInt(aTorus1, aTorus2, theTolTang);
+  Standard_Boolean bRet = anInt.IsDone();
+  if (bRet) {
+    if (anInt.TypeInter() == IntAna_Same) {
+      bEmpty = Standard_False;
+      bSameSurf = Standard_True;
+    } else {
+      bRet = TreatResultTorus(theQuad1, theQuad2, anInt, bEmpty, theSeqLin);
+    }
+  }
+  //
+  return bRet;
+}
+
+//=======================================================================
+//function : TreatResultTorus
+//purpose  : 
+//=======================================================================
+static Standard_Boolean TreatResultTorus(const IntSurf_Quadric& theQuad1,
+                                         const IntSurf_Quadric& theQuad2,
+                                         const IntAna_QuadQuadGeo& anInt,
+                                         Standard_Boolean& bEmpty,
+                                         IntPatch_SequenceOfLine& theSeqLin)
+{
+  Standard_Boolean bRet = anInt.IsDone();
+  //
+  if (!bRet) {
+    return bRet;
+  }
+  //
+  IntAna_ResultType typint = anInt.TypeInter();
+  Standard_Integer NbSol = anInt.NbSolutions();
+  bEmpty = Standard_False;
+  //
+  switch (typint) {
+  case IntAna_Empty :
+    bEmpty = Standard_True;
+    break;
+  //
+  case IntAna_Circle : {
+    Standard_Integer i;
+    IntSurf_TypeTrans trans1, trans2;
+    gp_Vec Tgt;
+    gp_Pnt ptref;
+    //
+    for (i = 1; i <= NbSol; ++i) {
+      gp_Circ aC = anInt.Circle(i);
+      if (theQuad1.TypeQuadric() == theQuad2.TypeQuadric()) {
+        AdjustToSeam(theQuad1.Torus(), aC);
+      }
+      ElCLib::D1(0., aC, ptref, Tgt);
+      Standard_Real qwe = Tgt.DotCross(theQuad2.Normale(ptref),
+                                       theQuad1.Normale(ptref));
+      if(qwe> 0.00000001) {
+        trans1 = IntSurf_Out;
+        trans2 = IntSurf_In;
+      }
+      else if(qwe< -0.00000001) {
+        trans1 = IntSurf_In;
+        trans2 = IntSurf_Out;
+      }
+      else { 
+        trans1=trans2=IntSurf_Undecided; 
+      }
+      //
+      Handle(IntPatch_GLine) glig = 
+        new IntPatch_GLine(aC, Standard_False, trans1, trans2);
+      theSeqLin.Append(glig);
+    }
+  }
+    break;
+  //
+  case IntAna_NoGeometricSolution:
+  default:
+    bRet = Standard_False;
+    break;
+  }
+  //
+  return bRet;
+}
 
     case GeomAbs_Cone: ts2 = 1; break;
     default: break;
   }
-
+  //
+  // treatment of the cases with torus and any other geom surface
+  if ((typs1 == GeomAbs_Torus && ts2) ||
+      (typs2 == GeomAbs_Torus && ts1) ||
+      (typs1 == GeomAbs_Torus && typs2 == GeomAbs_Torus)) {
+    // check if axes collinear
+    //
+    const Handle(Adaptor3d_HSurface)& aTorSurf = 
+      (typs1 == GeomAbs_Torus) ? theS1 : theS2;
+    const Handle(Adaptor3d_HSurface)& aGeomSurf = 
+      (typs1 == GeomAbs_Torus) ? theS2 : theS1;
+    //
+    Standard_Boolean bValid = 
+      aTorSurf->Torus().MajorRadius() > aTorSurf->Torus().MinorRadius();
+    if (bValid && (typs1 == typs2)) {
+      bValid = aGeomSurf->Torus().MajorRadius() > aGeomSurf->Torus().MinorRadius();
+    }
+    //
+    if (bValid) {
+      Standard_Boolean bCheck, bImpImp;
+      const gp_Ax1 aTorAx = aTorSurf->Torus().Axis();
+      const gp_Lin aL1(aTorAx);
+      //
+      bCheck = Standard_True;
+      bImpImp = Standard_False;
+      //
+      gp_Ax1 aGeomAx;
+      switch (aGeomSurf->GetType()) {
+      case GeomAbs_Plane: {
+        aGeomAx = aGeomSurf->Plane().Axis();
+        if (aTorAx.IsParallel(aGeomAx, Precision::Angular()) ||
+            (aTorAx.IsNormal(aGeomAx, Precision::Angular()) && 
+             (aGeomSurf->Plane().Distance(aTorAx.Location()) < Precision::Confusion()))) {
+          bImpImp = Standard_True;
+        }
+        bCheck = Standard_False;
+        break;
+      }
+      case GeomAbs_Sphere: {
+        if (aL1.Distance(aGeomSurf->Sphere().Location()) < Precision::Confusion()) {
+          bImpImp = Standard_True;
+        }
+        bCheck = Standard_False;
+        break;
+      }
+      case GeomAbs_Cylinder:
+        aGeomAx = aGeomSurf->Cylinder().Axis();
+        break;
+      case GeomAbs_Cone: 
+        aGeomAx = aGeomSurf->Cone().Axis();
+        break;
+      case GeomAbs_Torus: 
+        aGeomAx = aGeomSurf->Torus().Axis();
+        break;
+      default: 
+        bCheck = Standard_False;
+        break;
+      }
+      //
+      if (bCheck) {
+        if (aTorAx.IsParallel(aGeomAx, Precision::Angular()) &&
+            (aL1.Distance(aGeomAx.Location()) <= Precision::Confusion())) {
+          bImpImp = Standard_True;
+        }
+      }
+      //
+      if (bImpImp) {
+        ts1 = 1;
+        ts2 = 1;
+      }
+    }
+  }
+  //
   // Possible intersection types: 1. ts1 == ts2 == 1 <Geom-Geom>
   //                              2. ts1 != ts2      <Geom-Param>
   //                              3. ts1 == ts2 == 0 <Param-Param>
     case GeomAbs_Cone: ts2 = 1; break;
     default: break;
   }
-
+  //
+  // treatment of the cases with torus and any other geom surface
+  if ((typs1 == GeomAbs_Torus && ts2) ||
+      (typs2 == GeomAbs_Torus && ts1) ||
+      (typs1 == GeomAbs_Torus && typs2 == GeomAbs_Torus)) {
+    // check if axes collinear
+    //
+    const Handle(Adaptor3d_HSurface)& aTorSurf = 
+      (typs1 == GeomAbs_Torus) ? theS1 : theS2;
+    const Handle(Adaptor3d_HSurface)& aGeomSurf = 
+      (typs1 == GeomAbs_Torus) ? theS2 : theS1;
+    //
+    Standard_Boolean bValid = 
+      aTorSurf->Torus().MajorRadius() > aTorSurf->Torus().MinorRadius();
+    if (bValid && (typs1 == typs2)) {
+      bValid = aGeomSurf->Torus().MajorRadius() > aGeomSurf->Torus().MinorRadius();
+    }
+    //
+    if (bValid) {
+      Standard_Boolean bCheck, bImpImp;
+      const gp_Ax1 aTorAx = aTorSurf->Torus().Axis();
+      const gp_Lin aL1(aTorAx);
+      //
+      bCheck = Standard_True;
+      bImpImp = Standard_False;
+      //
+      gp_Ax1 aGeomAx;
+      switch (aGeomSurf->GetType()) {
+      case GeomAbs_Plane: {
+        aGeomAx = aGeomSurf->Plane().Axis();
+        if (aTorAx.IsParallel(aGeomAx, Precision::Angular()) ||
+            (aTorAx.IsNormal(aGeomAx, Precision::Angular()) && 
+             (aGeomSurf->Plane().Distance(aTorAx.Location()) < Precision::Confusion()))) {
+          bImpImp = Standard_True;
+        }
+        bCheck = Standard_False;
+        break;
+      }
+      case GeomAbs_Sphere: {
+        if (aL1.Distance(aGeomSurf->Sphere().Location()) < Precision::Confusion()) {
+          bImpImp = Standard_True;
+        }
+        bCheck = Standard_False;
+        break;
+      }
+      case GeomAbs_Cylinder:
+        aGeomAx = aGeomSurf->Cylinder().Axis();
+        break;
+      case GeomAbs_Cone: 
+        aGeomAx = aGeomSurf->Cone().Axis();
+        break;
+      case GeomAbs_Torus: 
+        aGeomAx = aGeomSurf->Torus().Axis();
+        break;
+      default: 
+        bCheck = Standard_False;
+        break;
+      }
+      //
+      if (bCheck) {
+        if (aTorAx.IsParallel(aGeomAx, Precision::Angular()) &&
+            (aL1.Distance(aGeomAx.Location()) <= Precision::Confusion())) {
+          bImpImp = Standard_True;
+        }
+      }
+      //
+      if (bImpImp) {
+        ts1 = 1;
+        ts2 = 1;
+      }
+    }
+  }
+  //
   // Possible intersection types: 1. ts1 == ts2 == 1 <Geom-Geom>
   //                              2. ts1 != ts2      <Geom-Param>
   //                              3. ts1 == ts2 == 0 <Param-Param>
             Quad1.SetValue(theS1->Cone());
             break;
 
+          case GeomAbs_Torus:
+            Quad1.SetValue(theS1->Torus());
+            break;
+
           default:
             break;
           }
             Quad2.SetValue(theS2->Cone());
             break;
 
+          case GeomAbs_Torus:
+            Quad2.SetValue(theS2->Torus());
+            break;
+
           default:
             break;
           }
 
 //=======================================================================
 
 static void Parameters(const Handle(Adaptor3d_HSurface)& myHS1,
-                      const Handle(Adaptor3d_HSurface)& myHS2,
-                      const gp_Pnt& Ptref,
-                      Standard_Real& U1,
-                      Standard_Real& V1,
-                      Standard_Real& U2,
-                      Standard_Real& V2)
+                       const Handle(Adaptor3d_HSurface)& myHS2,
+                       const gp_Pnt& Ptref,
+                       Standard_Real& U1,
+                       Standard_Real& V1,
+                       Standard_Real& U2,
+                       Standard_Real& V2)
 {
   IntSurf_Quadric quad1,quad2;
   GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
   case GeomAbs_Sphere:
     quad1.SetValue(myHS1->Surface().Sphere());
     break;
+  case GeomAbs_Torus:
+    quad1.SetValue(myHS1->Surface().Torus());
+    break;
   default:
     Standard_ConstructionError::Raise("IntPatch_IntSS::MakeCurve");
   }
   case GeomAbs_Sphere:
     quad2.SetValue(myHS2->Surface().Sphere());
     break;
+  case GeomAbs_Torus:
+    quad2.SetValue(myHS2->Surface().Torus());
+    break;
   default:
     Standard_ConstructionError::Raise("IntPatch_IntSS::MakeCurve");
   }
 
      Pln         from gp,
      Cylinder    from gp,
      Sphere      from gp,
-     Cone        from gp,
+     Cone        from gp, 
+     Torus       from gp,
      SurfaceType from GeomAbs     
 
 is
     Create(C: Cone from gp)
     
        returns Quadric from IntSurf;
+ 
+ 
+    Create(T: Torus from gp)
+    
+       returns Quadric from IntSurf;
 
 
     SetValue(me: in out; P: Pln from gp)
     SetValue(me: in out; C: Cone from gp)
     
        is static;
+ 
+    SetValue(me: in out; T: Torus from gp)
+    
+       is static;
 
 
     Distance(me; P: Pnt from gp)
        ---C++: inline
        
        is static;
+ 
+ 
+    Torus(me)
+    
+       returns Torus from gp
+       ---C++: inline
+       
+       is static;
 
 
 
 
 #include <gp.hxx>
 
 
-
 // ============================================================
 IntSurf_Quadric::IntSurf_Quadric ():typ(GeomAbs_OtherSurface),
    prm1(0.), prm2(0.), prm3(0.), prm4(0.)
   prm4 = 0.0;
 }
 // ============================================================
+IntSurf_Quadric::IntSurf_Quadric (const gp_Torus& T):
+
+       ax3(T.Position()),typ(GeomAbs_Torus)
+{
+  ax3direc = ax3.Direct();
+  lin.SetPosition(ax3.Axis());
+  prm1 = T.MajorRadius();
+  prm2 = T.MinorRadius();
+  prm3 = 0.0;
+  prm4 = 0.0;
+}
+// ============================================================
 void IntSurf_Quadric::SetValue (const gp_Pln& P)
 {
   typ = GeomAbs_Plane;
   prm4 = 0.0;
 }
 // ============================================================
+void IntSurf_Quadric::SetValue (const gp_Torus& T)
+{
+  typ = GeomAbs_Torus;
+  ax3 = T.Position();
+  ax3direc = ax3.Direct();
+  lin.SetPosition(ax3.Axis());
+  prm1 = T.MajorRadius();
+  prm2 = T.MinorRadius();
+  prm3 = 0.0;
+  prm4 = 0.0;
+}
+// ============================================================
 Standard_Real IntSurf_Quadric::Distance (const gp_Pnt& P) const {
   switch (typ) {
   case GeomAbs_Plane:   // plan
       dist = (dist-distp)/prm3;
       return(dist);
     }
+  case GeomAbs_Torus:   // torus
+    {
+      gp_Pnt O, Pp, PT;
+      //
+      O = ax3.Location();
+      gp_Vec OZ (ax3.Direction());
+      Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
+      //
+      gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ? 
+        ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
+      PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
+      //
+      Standard_Real dist = P.Distance(PT) - prm2;
+      return dist;
+    }
   default:
     {
     }
       ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
       grad=D1u.Crossed(D1v);
       if(ax3direc==Standard_False) { 
-       grad.Reverse();
+        grad.Reverse();
       }
       grad.Normalize();
     }
     break;
+  case GeomAbs_Torus:   // torus
+    {
+      gp_Pnt O, Pp, PT;
+      //
+      O = ax3.Location();
+      gp_Vec OZ (ax3.Direction());
+      Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
+      //
+      gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ? 
+        ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
+      PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
+      //
+      grad.SetXYZ(P.XYZ() - PT.XYZ());
+      Standard_Real N = grad.Magnitude();
+      if(N>1e-14) { grad.Divide(N); } 
+      else { grad.SetCoord(0., 0., 0.); }
+    }
+    break;
   default:
     {}
     break;
 }
 // ============================================================
 void IntSurf_Quadric::ValAndGrad (const gp_Pnt& P,
-                                  Standard_Real& Dist,
-                                  gp_Vec& Grad) const
+                                  Standard_Real& Dist,
+                                  gp_Vec& Grad) const
 {
 
   switch (typ) {
       Dist = dist;
       Grad=D1u.Crossed(D1v);
       if(ax3direc==Standard_False) { 
-       Grad.Reverse();
+        Grad.Reverse();
       }
       //-- lbr le 7 mars 96 
       //-- Si le gardient est nul, on est sur l axe 
       //-- et dans ce cas dist vaut 0 
       //-- On peut donc renvoyer une valeur quelconque. 
       if(   Grad.X() > 1e-13 
-        || Grad.Y() > 1e-13 
-        || Grad.Z() > 1e-13) { 
-       Grad.Normalize();
+         || Grad.Y() > 1e-13 
+         || Grad.Z() > 1e-13) { 
+        Grad.Normalize();
       }
     }
     break;
+  case GeomAbs_Torus:   
+    {
+      gp_Pnt O, Pp, PT;
+      //
+      O = ax3.Location();
+      gp_Vec OZ (ax3.Direction());
+      Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
+      //
+      gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ? 
+        ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
+      PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
+      //
+      Dist = P.Distance(PT) - prm2;
+      //
+      Grad.SetXYZ(P.XYZ()-PT.XYZ());
+      Standard_Real N = Grad.Magnitude();
+      if(N>1e-14) { Grad.Divide(N); } 
+      else { Grad.SetCoord(0., 0., 0.); }
+    }
+    break;
   default:
     {}
     break;
 }
 // ============================================================
 gp_Pnt IntSurf_Quadric::Value(const Standard_Real U,
-                              const Standard_Real V) const
+                              const Standard_Real V) const
 {
   switch (typ) {
 
     return ElSLib::SphereValue(U,V,ax3,prm1);
   case GeomAbs_Cone:  
     return ElSLib::ConeValue(U,V,ax3,prm1,prm2);
+  case GeomAbs_Torus:  
+    return ElSLib::TorusValue(U,V,ax3,prm1,prm2);
   default:
     {
       gp_Pnt p(0,0,0);
 }
 // ============================================================
 void IntSurf_Quadric::D1(const Standard_Real U,
-                        const Standard_Real V,
-                        gp_Pnt& P,
-                        gp_Vec& D1U,
-                        gp_Vec& D1V) const
+                         const Standard_Real V,
+                         gp_Pnt& P,
+                         gp_Vec& D1U,
+                         gp_Vec& D1V) const
 {
   switch (typ) {
   case GeomAbs_Plane:  
   case GeomAbs_Cone:  
     ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
     break;
+  case GeomAbs_Torus:  
+    ElSLib::TorusD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
+    break;
   default:
     {
     }
 }
 // ============================================================
 gp_Vec IntSurf_Quadric::DN(const Standard_Real U,
-                           const Standard_Real V,
-                           const Standard_Integer Nu,
-                           const Standard_Integer Nv) const
+                           const Standard_Real V,
+                           const Standard_Integer Nu,
+                           const Standard_Integer Nv) const
 {
   switch (typ) {
   case GeomAbs_Plane: 
     return ElSLib::SphereDN(U,V,ax3,prm1,Nu,Nv);
   case GeomAbs_Cone:  
     return ElSLib::ConeDN(U,V,ax3,prm1,prm2,Nu,Nv);
+  case GeomAbs_Torus:  
+    return ElSLib::TorusDN(U,V,ax3,prm1,prm2,Nu,Nv);
   default:
     {
       gp_Vec v(0,0,0);
 }
 // ============================================================
 gp_Vec IntSurf_Quadric::Normale(const Standard_Real U,
-                                const Standard_Real V) const
+                                const Standard_Real V) const
 {
   switch (typ) {
   case GeomAbs_Plane:  
       gp_Vec D1u,D1v;
       ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1u,D1v);
       if(D1u.Magnitude()<0.0000001) { 
-       gp_Vec Vn(0.0,0.0,0.0); 
-       return(Vn);
+        gp_Vec Vn(0.0,0.0,0.0); 
+        return(Vn);
       }
       return(D1u.Crossed(D1v));
     }
+  case GeomAbs_Torus:
+    return Normale(Value(U,V));
   default:     
     {
       gp_Vec v(0,0,0);
   case GeomAbs_Cylinder:   
     {
       if(ax3direc) { 
-       return lin.Normal(P).Direction();
+        return lin.Normal(P).Direction();
       }
       else { 
-       gp_Dir D(lin.Normal(P).Direction());
-       D.Reverse();
-       return(D);
+        gp_Dir D(lin.Normal(P).Direction());
+        D.Reverse();
+        return(D);
       }
     }
   case GeomAbs_Sphere:   
     {
       if(ax3direc) { 
-       gp_Vec ax3P(ax3.Location(),P);
-       return gp_Dir(ax3P);
+        gp_Vec ax3P(ax3.Location(),P);
+        return gp_Dir(ax3P);
       }
       else { 
-       gp_Vec Pax3(P,ax3.Location());
-       return gp_Dir(Pax3);
+        gp_Vec Pax3(P,ax3.Location());
+        return gp_Dir(Pax3);
       }
     }
   case GeomAbs_Cone:   
     {
       Standard_Real U,V;
       ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
-      return Normale(U,V);;
+      return Normale(U,V);
+    }
+  case GeomAbs_Torus:
+    {
+      gp_Pnt O, Pp, PT;
+      //
+      O = ax3.Location();
+      gp_Vec OZ (ax3.Direction());
+      Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
+      //
+      gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ? 
+        ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
+      PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
+      if (PT.SquareDistance(P) < 1e-14) {
+        return gp_Dir(OZ);
+      }
+      gp_Dir aD(ax3direc ? gp_Vec(PT, P) : gp_Vec(P, PT));
+      return aD;
     }
   default:      
     {
 }
 // ============================================================
 void IntSurf_Quadric::Parameters (const gp_Pnt& P,
-                                  Standard_Real& U,
-                                  Standard_Real& V) const
+                                  Standard_Real& U,
+                                  Standard_Real& V) const
 {
   switch (typ) {
   case GeomAbs_Plane:   
     ElSLib::SphereParameters(ax3,prm1,P,U,V);
     break;
   case GeomAbs_Cone: 
-    {
-      ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
-    }
+    ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
+    break;
+  case GeomAbs_Torus: 
+    ElSLib::TorusParameters(ax3,prm1,prm2,P,U,V);
     break;
   default:
-    {}
     break;
   }
 }
 
   return gp_Cone(ax3,prm2,prm1);
 }
 
+inline gp_Torus IntSurf_Quadric::Torus () const {
+
+  return gp_Torus(ax3, prm1, prm2);
+}
+
 
 
     case GeomAbs_Sphere:
       quad1.SetValue(myHS1->Surface().Sphere());
       break;
+    case GeomAbs_Torus:
+      quad1.SetValue(myHS1->Surface().Torus());
+      break;
     default:
       Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
     }
     case GeomAbs_Sphere:
       quad2.SetValue(myHS2->Surface().Sphere());
       break;
+    case GeomAbs_Torus:
+      quad2.SetValue(myHS2->Surface().Torus());
+      break;
     default:
       Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
     }
   case GeomAbs_Sphere:
     quad1.SetValue(HS1->Surface().Sphere());
     break;
+  case GeomAbs_Torus:
+    quad1.SetValue(HS1->Surface().Torus());
+    break;
   default:
     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
   }
   case GeomAbs_Sphere:
     quad2.SetValue(HS2->Surface().Sphere());
     break;
+  case GeomAbs_Torus:
+    quad2.SetValue(HS2->Surface().Torus());
+    break;
   default:
     Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
   }
 
 //purpose  : 
 //=======================================================================
 void Parameters(const Handle(GeomAdaptor_HSurface)& myHS1,
-               const gp_Pnt& Ptref,
-               Standard_Real& U1,
-               Standard_Real& V1)
+                const gp_Pnt& Ptref,
+                Standard_Real& U1,
+                Standard_Real& V1)
 {
   IntSurf_Quadric quad1;
   //
     case GeomAbs_Sphere:   
       quad1.SetValue(myHS1->Surface().Sphere()); 
       break;
+    case GeomAbs_Torus:
+      quad1.SetValue(myHS1->Surface().Torus()); 
+      break;
     default: 
       Standard_ConstructionError::Raise("IntTools_LineConstructor::Parameters");
   }
 
--- /dev/null
+puts "=========="
+puts "OCC24470"
+puts "=========="
+puts ""
+##################################################
+# Wrong result done by General Fuse algorithm
+##################################################
+
+restore [locate_data_file bug24470_box.brep] b1
+restore [locate_data_file bug24470_cone.brep] b2
+restore [locate_data_file bug24470_cylinder.brep] b3
+restore [locate_data_file bug24470_rotBox.brep] b4
+restore [locate_data_file bug24470_rotCyl.brep] b5
+restore [locate_data_file bug24470_sphere.brep] b6
+restore [locate_data_file bug24470_torus.brep] b7
+
+bclearobjects
+bcleartools
+baddobjects b1 b2 b3 b4 b5 b6 b7
+bfillds
+bbuild result
+checkshape result
+
+set square 595443
+set 2dviewer 1
+
+
+