// Created on: 1992-05-07
// Created by: Jacques GOUSSARD
// Copyright (c) 1992-1999 Matra Datavision
// Copyright (c) 1999-2014 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 License 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.

#include <algorithm>
#include <Bnd_Range.hxx>
#include <IntAna_ListOfCurve.hxx>
#include <math_Matrix.hxx>
#include <NCollection_IncAllocator.hxx>
#include <Standard_DivideByZero.hxx>
#include <math_Vector.hxx>

//If Abs(a) <= aNulValue then it is considered that a = 0.
static const Standard_Real aNulValue = 1.0e-11;

static void ShortCosForm( const Standard_Real theCosFactor,
                          const Standard_Real theSinFactor,
                          Standard_Real& theCoeff,
                          Standard_Real& theAngle);
//
static Standard_Boolean ExploreCurve(const gp_Cone& theCo,
                                     IntAna_Curve& aC,
                                     const Standard_Real aTol,
                                     IntAna_ListOfCurve& aLC);

static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
                                      const Standard_Real theUlTarget,
                                      Standard_Real& theUGiven,
                                      const Standard_Real theTol2D,
                                      const Standard_Real thePeriod,
                                      const Standard_Boolean theFlForce);


class ComputationMethods
{
  //Every cylinder can be represented by the following equation in parametric form:
  //    S(U,V) = L + R*cos(U)*Xd+R*sin(U)*Yd+V*Zd,
  //where location L, directions Xd, Yd and Zd have type gp_XYZ.

  //Intersection points between two cylinders can be found from the following system:
  //    S1(U1, V1) = S2(U2, V2)
  //or
  //    {X01 + R1*cos(U1)*Xx1 + R1*sin(U1)*Yx1 + V1*Zx1 = X02 + R2*cos(U2)*Xx2 + R2*sin(U2)*Yx2 + V2*Zx2
  //    {Y01 + R1*cos(U1)*Xy1 + R1*sin(U1)*Yy1 + V1*Zy1 = Y02 + R2*cos(U2)*Xy2 + R2*sin(U2)*Yy2 + V2*Zy2   (1)
  //    {Z01 + R1*cos(U1)*Xz1 + R1*sin(U1)*Yz1 + V1*Zz1 = Z02 + R2*cos(U2)*Xz2 + R2*sin(U2)*Yz2 + V2*Zz2

  //The formula (1) can be rewritten as follows
  //    {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1
  //    {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2                                   (2)
  //    {C13*V1+C23*V2=A13*cos(U1)+B13*sin(U1)+A23*cos(U2)+B23*sin(U2)+D3

  //Hereafter we consider that in system
  //    {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1                                   (3)
  //    {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2
  //variables V1 and V2 can be found unambiguously, i.e. determinant
  //            |C11 C21|
  //            |       | != 0
  //            |C12 C22|
  //            
  //In this case, variables V1 and V2 can be found as follows:
  //    {V1 = K11*sin(U1)+K21*sin(U2)+L11*cos(U1)+L21*cos(U2)+M1 = K1*cos(U1-FIV1)+L1*cos(U2-PSIV1)+M1      (4)
  //    {V2 = K12*sin(U1)+K22*sin(U2)+L12*cos(U1)+L22*cos(U2)+M2 = K2*cos(U2-FIV2)+L2*cos(U2-PSIV2)+M2

  //Having substituted result of (4) to the 3rd equation of (2), we will obtain equation
  //    cos(U2-FI2) = B*cos(U1-FI1)+C.                                                                      (5)

  //I.e. when U1 is taken different given values (from domain), correspond U2 value can be computed
  //from equation (5). After that, V1 and V2 can be computed from the system (4) (see
  //CylCylComputeParameters(...) methods).

  //It is important to remark that equation (5) (in general) has two solutions: U2=FI2 +/- f(U1).
  //Therefore, we are getting here two intersection lines.

public:
  //Stores equations coefficients
  struct stCoeffsValue
  {
    stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);

    math_Vector mVecA1;
    math_Vector mVecA2;
    math_Vector mVecB1;
    math_Vector mVecB2;
    math_Vector mVecC1;
    math_Vector mVecC2;
    math_Vector mVecD;

    Standard_Real mK21; //sinU2
    Standard_Real mK11; //sinU1
    Standard_Real mL21; //cosU2
    Standard_Real mL11;  //cosU1
    Standard_Real mM1;  //Free member

    Standard_Real mK22; //sinU2
    Standard_Real mK12; //sinU1
    Standard_Real mL22; //cosU2
    Standard_Real mL12; //cosU1
    Standard_Real mM2; //Free member

    Standard_Real mK1;
    Standard_Real mL1;
    Standard_Real mK2;
    Standard_Real mL2;

    Standard_Real mFIV1;
    Standard_Real mPSIV1;
    Standard_Real mFIV2;
    Standard_Real mPSIV2;

    Standard_Real mB;
    Standard_Real mC;
    Standard_Real mFI1;
    Standard_Real mFI2;
  };


  //! Determines, if U2(U1) function is increasing.
  static Standard_Boolean CylCylMonotonicity(const Standard_Real theU1par,
                                             const Standard_Integer theWLIndex,
                                             const stCoeffsValue& theCoeffs,
                                             const Standard_Real thePeriod,
                                             Standard_Boolean& theIsIncreasing);

  //! Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
  //! esimates the tolerance of U2-computing (estimation result is
  //! assigned to *theDelta value).
  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
                                                const Standard_Integer theWLIndex,
                                                const stCoeffsValue& theCoeffs,
                                                Standard_Real& theU2,
                                                Standard_Real* const theDelta = 0);

  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
                                                  const Standard_Real theU2,
                                                  const stCoeffsValue& theCoeffs,
                                                  Standard_Real& theV1,
                                                  Standard_Real& theV2);

  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
                                                  const Standard_Integer theWLIndex,
                                                  const stCoeffsValue& theCoeffs,
                                                  Standard_Real& theU2,
                                                  Standard_Real& theV1,
                                                  Standard_Real& theV2);

};

ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1,
                                                 const gp_Cylinder& theCyl2):
  mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()),
  mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()),
  mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()),
  mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()),
  mVecC1(theCyl1.Axis().Direction().XYZ()),
  mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()),
  mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ())
{
  enum CoupleOfEquation
  {
    COENONE = 0,
    COE12 = 1,
    COE23 = 2,
    COE13 = 3
  }aFoundCouple = COENONE;


  Standard_Real aDetV1V2 = 0.0;

  const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2
  const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3
  const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3
  const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
  const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
  const Standard_Real anAbsD3 = Abs(aDelta3); //1-3

  if(anAbsD1 >= anAbsD2)
  {
    if(anAbsD3 > anAbsD1)
    {
      aFoundCouple = COE13;
      aDetV1V2 = aDelta3;
    }
    else
    {
      aFoundCouple = COE12;
      aDetV1V2 = aDelta1;
    }
  }
  else
  {
    if(anAbsD3 > anAbsD2)
    {
      aFoundCouple = COE13;
      aDetV1V2 = aDelta3;
    }
    else
    {
      aFoundCouple = COE23;
      aDetV1V2 = aDelta2;
    }
  }

  // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is
  // cross-product between directions (i.e. sine of angle).
  // If sine is too small then sine is (approx.) equal to angle itself.
  // Therefore, in this case we should compare sine with angular tolerance. 
  // This constant is used for check if axes are parallel (see constructor
  // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file).
  if(Abs(aDetV1V2) < Precision::Angular())
  {
    throw Standard_Failure("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
  }

  switch(aFoundCouple)
  {
  case COE12:
    break;
  case COE23:
    {
      math_Vector aVTemp(mVecA1);
      mVecA1(1) = aVTemp(2);
      mVecA1(2) = aVTemp(3);
      mVecA1(3) = aVTemp(1);

      aVTemp = mVecA2;
      mVecA2(1) = aVTemp(2);
      mVecA2(2) = aVTemp(3);
      mVecA2(3) = aVTemp(1);

      aVTemp = mVecB1;
      mVecB1(1) = aVTemp(2);
      mVecB1(2) = aVTemp(3);
      mVecB1(3) = aVTemp(1);

      aVTemp = mVecB2;
      mVecB2(1) = aVTemp(2);
      mVecB2(2) = aVTemp(3);
      mVecB2(3) = aVTemp(1);

      aVTemp = mVecC1;
      mVecC1(1) = aVTemp(2);
      mVecC1(2) = aVTemp(3);
      mVecC1(3) = aVTemp(1);

      aVTemp = mVecC2;
      mVecC2(1) = aVTemp(2);
      mVecC2(2) = aVTemp(3);
      mVecC2(3) = aVTemp(1);

      aVTemp = mVecD;
      mVecD(1) = aVTemp(2);
      mVecD(2) = aVTemp(3);
      mVecD(3) = aVTemp(1);

      break;
    }
  case COE13:
    {
      math_Vector aVTemp = mVecA1;
      mVecA1(2) = aVTemp(3);
      mVecA1(3) = aVTemp(2);

      aVTemp = mVecA2;
      mVecA2(2) = aVTemp(3);
      mVecA2(3) = aVTemp(2);

      aVTemp = mVecB1;
      mVecB1(2) = aVTemp(3);
      mVecB1(3) = aVTemp(2);

      aVTemp = mVecB2;
      mVecB2(2) = aVTemp(3);
      mVecB2(3) = aVTemp(2);

      aVTemp = mVecC1;
      mVecC1(2) = aVTemp(3);
      mVecC1(3) = aVTemp(2);

      aVTemp = mVecC2;
      mVecC2(2) = aVTemp(3);
      mVecC2(3) = aVTemp(2);

      aVTemp = mVecD;
      mVecD(2) = aVTemp(3);
      mVecD(3) = aVTemp(2);

      break;
    }
  default:
    break;
  }

  //------- For V1 (begin)
  //sinU2
  mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
  //sinU1
  mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
  //cosU2
  mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
  //cosU1
  mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
  //Free member
  mM1 =  (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
  //------- For V1 (end)

  //------- For V2 (begin)
  //sinU2
  mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
  //sinU1
  mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
  //cosU2
  mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
  //cosU1
  mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
  //Free member
  mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
  //------- For V1 (end)

  ShortCosForm(mL11, mK11, mK1, mFIV1);
  ShortCosForm(mL21, mK21, mL1, mPSIV1);
  ShortCosForm(mL12, mK12, mK2, mFIV2);
  ShortCosForm(mL22, mK22, mL2, mPSIV2);

  const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2
    aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2
    aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1
    aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1

  mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2;  //Free

  Standard_Real aA = 0.0;

  ShortCosForm(aB2,aB1,mB,mFI1);
  ShortCosForm(aA2,aA1,aA,mFI2);

  mB /= aA;
  mC /= aA;
}

class WorkWithBoundaries
{
public:
  enum SearchBoundType
  {
    SearchNONE = 0,
    SearchV1 = 1,
    SearchV2 = 2
  };

  struct StPInfo
  {
    StPInfo()
    {
      mySurfID = 0;
      myU1 = RealLast();
      myV1 = RealLast();
      myU2 = RealLast();
      myV2 = RealLast();
    }

    //Equal to 0 for 1st surface non-zero for 2nd one.
    Standard_Integer mySurfID;

    Standard_Real myU1;
    Standard_Real myV1;
    Standard_Real myU2;
    Standard_Real myV2;

    bool operator>(const StPInfo& theOther) const
    {
      return myU1 > theOther.myU1;
    }

    bool operator<(const StPInfo& theOther) const
    {
      return myU1 < theOther.myU1;
    }

    bool operator==(const StPInfo& theOther) const
    {
      return myU1 == theOther.myU1;
    }
  };

  WorkWithBoundaries(const IntSurf_Quadric& theQuad1,
                     const IntSurf_Quadric& theQuad2,
                     const ComputationMethods::stCoeffsValue& theCoeffs,
                     const Bnd_Box2d& theUVSurf1,
                     const Bnd_Box2d& theUVSurf2,
                     const Standard_Integer theNbWLines,
                     const Standard_Real thePeriod,
                     const Standard_Real theTol3D,
                     const Standard_Real theTol2D,
                     const Standard_Boolean isTheReverse) :
      myQuad1(theQuad1), myQuad2(theQuad2), myCoeffs(theCoeffs),
      myUVSurf1(theUVSurf1), myUVSurf2(theUVSurf2), myNbWLines(theNbWLines),
      myPeriod(thePeriod), myTol3D(theTol3D), myTol2D(theTol2D),
      myIsReverse(isTheReverse)
  {
  };

  // Returns parameters of system solved while finding
  // intersection line
  const ComputationMethods::stCoeffsValue &SICoeffs() const
  {
    return myCoeffs;
  }

  // Returns quadric correspond to the index theIdx.
  const IntSurf_Quadric& GetQSurface(const Standard_Integer theIdx) const
  {
    if (theIdx <= 1)
      return myQuad1;

    return myQuad2;
  }

  // Returns TRUE in case of reverting surfaces
  Standard_Boolean IsReversed() const
  {
    return myIsReverse;
  }

  // Returns 2D-tolerance
  Standard_Real Get2dTolerance() const
  {
    return myTol2D;
  }

  // Returns 3D-tolerance
  Standard_Real Get3dTolerance() const
  {
    return myTol3D;
  }

  // Returns UV-bounds of 1st surface
  const Bnd_Box2d& UVS1() const
  {
    return myUVSurf1;
  }

  // Returns UV-bounds of 2nd surface
  const Bnd_Box2d& UVS2() const
  {
    return myUVSurf2;
  }

  void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
                        const Standard_Real theU1,
                        const Standard_Real theU1Min,
                        const Standard_Real theU2,
                        const Standard_Real theV1,
                        const Standard_Real theV1Prev,
                        const Standard_Real theV2,
                        const Standard_Real theV2Prev,
                        const Standard_Integer theWLIndex,
                        const Standard_Boolean theFlForce,
                        Standard_Boolean& isTheFound1,
                        Standard_Boolean& isTheFound2) const;
  
  static Standard_Boolean BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
                                              const Standard_Real thePeriod,
                                              Bnd_Range theURange[]);
  
  void BoundaryEstimation(const gp_Cylinder& theCy1,
                          const gp_Cylinder& theCy2,
                          Bnd_Range& theOutBoxS1,
                          Bnd_Range& theOutBoxS2) const;

protected:

  //Solves equation (2) (see declaration of ComputationMethods class) in case,
  //when V1 or V2 (is set by theSBType argument) is known (corresponds to the boundary
  //and equal to theVzad) but U1 is unknown. Computation is made by numeric methods and
  //requires initial values (theVInit, theInitU2 and theInitMainVar).
  Standard_Boolean
              SearchOnVBounds(const SearchBoundType theSBType,
                              const Standard_Real theVzad,
                              const Standard_Real theVInit,
                              const Standard_Real theInitU2,
                              const Standard_Real theInitMainVar,
                              Standard_Real& theMainVariableValue) const;

  const WorkWithBoundaries& operator=(const WorkWithBoundaries&);

private:
  friend class ComputationMethods;

  const IntSurf_Quadric& myQuad1;
  const IntSurf_Quadric& myQuad2;
  const ComputationMethods::stCoeffsValue& myCoeffs;
  const Bnd_Box2d& myUVSurf1;
  const Bnd_Box2d& myUVSurf2;
  const Standard_Integer myNbWLines;
  const Standard_Real myPeriod;
  const Standard_Real myTol3D;
  const Standard_Real myTol2D;
  const Standard_Boolean myIsReverse;
};

static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
                                  const IntSurf_Quadric& theQuad2,
                                  const Handle(IntSurf_LineOn2S)& theLine,
                                  const ComputationMethods::stCoeffsValue& theCoeffs,
                                  const Standard_Integer theWLIndex,
                                  const Standard_Integer theMinNbPoints,
                                  const Standard_Integer theStartPointOnLine,
                                  const Standard_Integer theEndPointOnLine,
                                  const Standard_Real theTol2D,
                                  const Standard_Real thePeriodOfSurf2,
                                  const Standard_Boolean isTheReverse);

//=======================================================================
//function : MinMax
//purpose  : Replaces theParMIN = MIN(theParMIN, theParMAX),
//                    theParMAX = MAX(theParMIN, theParMAX).
//=======================================================================
static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
{
  if(theParMIN > theParMAX)
  {
    const Standard_Real aux = theParMAX;
    theParMAX = theParMIN;
    theParMIN = aux;
  }
}

//=======================================================================
//function : ExtremaLineLine
//purpose  : Computes extrema between the given lines. Returns parameters
//          on correspond curve (see correspond method for Extrema_ExtElC class). 
//=======================================================================
static inline void ExtremaLineLine(const gp_Ax1& theC1,
                                   const gp_Ax1& theC2,
                                   const Standard_Real theCosA,
                                   const Standard_Real theSqSinA,
                                   Standard_Real& thePar1,
                                   Standard_Real& thePar2)
{
  const gp_Dir &aD1 = theC1.Direction(),
               &aD2 = theC2.Direction();

  const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ();
  const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2),
                      aD2L = aD2.XYZ().Dot(aL1L2);

  thePar1 = (aD1L - theCosA * aD2L) / theSqSinA;
  thePar2 = (theCosA * aD1L - aD2L) / theSqSinA;
}

//=======================================================================
//function : VBoundaryPrecise
//purpose  : By default, we shall consider, that V1 and V2 will be increased
//            if U1 is increased. But if it is not, new V1set and/or V2set
//            must be computed as [V1current - DeltaV1] (analogically
//            for V2). This function processes this case.
//=======================================================================
static void VBoundaryPrecise( const math_Matrix& theMatr,
                              const Standard_Real theV1AfterDecrByDelta,
                              const Standard_Real theV2AfterDecrByDelta,
                              Standard_Real& theV1Set,
                              Standard_Real& theV2Set)
{
  //Now we are going to define if V1 (and V2) increases
  //(or decreases) when U1 will increase.
  const Standard_Integer aNbDim = 3;
  math_Matrix aSyst(1, aNbDim, 1, aNbDim);

  aSyst.SetCol(1, theMatr.Col(1));
  aSyst.SetCol(2, theMatr.Col(2));
  aSyst.SetCol(3, theMatr.Col(4));

  //We have the system (see comment to StepComputing(...) function)
  //    {a11*dV1 + a12*dV2 + a14*dU2 = -a13*dU1
  //    {a21*dV1 + a22*dV2 + a24*dU2 = -a23*dU1
  //    {a31*dV1 + a32*dV2 + a34*dU2 = -a33*dU1

  const Standard_Real aDet = aSyst.Determinant();

  aSyst.SetCol(1, theMatr.Col(3));
  const Standard_Real aDet1 = aSyst.Determinant();

  aSyst.SetCol(1, theMatr.Col(1));
  aSyst.SetCol(2, theMatr.Col(3));

  const Standard_Real aDet2 = aSyst.Determinant();

  //Now,
  //    dV1 = -dU1*aDet1/aDet
  //    dV2 = -dU1*aDet2/aDet

  //If U1 is increased then dU1 > 0.
  //If (aDet1/aDet > 0) then dV1 < 0 and 
  //V1 will be decreased after increasing U1.

  //We have analogical situation with V2-parameter. 

  if(aDet*aDet1 > 0.0)
  {
    theV1Set = theV1AfterDecrByDelta;
  }

  if(aDet*aDet2 > 0.0)
  {
    theV2Set = theV2AfterDecrByDelta;
  }
}

//=======================================================================
//function : DeltaU1Computing
//purpose  : Computes new step for U1 parameter.
//=======================================================================
static inline 
        Standard_Boolean DeltaU1Computing(const math_Matrix& theSyst,
                                          const math_Vector& theFree,
                                          Standard_Real& theDeltaU1Found)
{
  Standard_Real aDet = theSyst.Determinant();

  if(Abs(aDet) > aNulValue)
  {
    math_Matrix aSyst1(theSyst);
    aSyst1.SetCol(2, theFree);

    theDeltaU1Found = Abs(aSyst1.Determinant()/aDet);
    return Standard_True;
  }

  return Standard_False;
}

//=======================================================================
//function : StepComputing
//purpose  : 
//
//Attention!!!:
//            theMatr must have 3*5-dimension strictly.
//            For system
//                {a11*V1+a12*V2+a13*dU1+a14*dU2=b1; 
//                {a21*V1+a22*V2+a23*dU1+a24*dU2=b2; 
//                {a31*V1+a32*V2+a33*dU1+a34*dU2=b3; 
//            theMatr must be following:
//                (a11 a12 a13 a14 b1) 
//                (a21 a22 a23 a24 b2) 
//                (a31 a32 a33 a34 b3) 
//=======================================================================
static Standard_Boolean StepComputing(const math_Matrix& theMatr,
                                      const Standard_Real theV1Cur,
                                      const Standard_Real theV2Cur,
                                      const Standard_Real theDeltaV1,
                                      const Standard_Real theDeltaV2,
                                      Standard_Real& theDeltaU1Found/*,
                                      Standard_Real& theDeltaU2Found,
                                      Standard_Real& theV1Found,
                                      Standard_Real& theV2Found*/)
{
#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
  bool flShow = false;

  if(flShow)
  {
    printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
              theMatr(1,1), theMatr(1,2), theMatr(1,3), theMatr(1,4), theMatr(1,5));
    printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
              theMatr(2,1), theMatr(2,2), theMatr(2,3), theMatr(2,4), theMatr(2,5));
    printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
              theMatr(3,1), theMatr(3,2), theMatr(3,3), theMatr(3,4), theMatr(3,5));
  }
#endif

  Standard_Boolean isSuccess = Standard_False;
  theDeltaU1Found/* = theDeltaU2Found*/ = RealLast();
  //theV1Found = theV1set;
  //theV2Found = theV2Set;
  const Standard_Integer aNbDim = 3;

  math_Matrix aSyst(1, aNbDim, 1, aNbDim);
  math_Vector aFree(1, aNbDim);

  //By default, increasing V1(U1) and V2(U1) functions is
  //considered
  Standard_Real aV1Set = theV1Cur + theDeltaV1,
                aV2Set = theV2Cur + theDeltaV2;

  //However, what is indeed?
  VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1,
                    theV2Cur - theDeltaV2, aV1Set, aV2Set);

  aSyst.SetCol(2, theMatr.Col(3));
  aSyst.SetCol(3, theMatr.Col(4));

  for(Standard_Integer i = 0; i < 2; i++)
  {
    if(i == 0)
    {//V1 is known
      aSyst.SetCol(1, theMatr.Col(2));
      aFree.Set(1, aNbDim, theMatr.Col(5)-aV1Set*theMatr.Col(1));
    }
    else
    {//i==1 => V2 is known
      aSyst.SetCol(1, theMatr.Col(1));
      aFree.Set(1, aNbDim, theMatr.Col(5)-aV2Set*theMatr.Col(2));
    }

    Standard_Real aNewDU = theDeltaU1Found;
    if(DeltaU1Computing(aSyst, aFree, aNewDU))
    {
      isSuccess = Standard_True;
      if(aNewDU < theDeltaU1Found)
      {
        theDeltaU1Found = aNewDU;
      }
    }
  }

  if(!isSuccess)
  {
    aFree = theMatr.Col(5) - aV1Set*theMatr.Col(1) - aV2Set*theMatr.Col(2);
    math_Matrix aSyst1(1, aNbDim, 1, 2);
    aSyst1.SetCol(1, aSyst.Col(2));
    aSyst1.SetCol(2, aSyst.Col(3));

    //Now we have overdetermined system.

    const Standard_Real aDet1 = theMatr(1,3)*theMatr(2,4) - theMatr(2,3)*theMatr(1,4);
    const Standard_Real aDet2 = theMatr(1,3)*theMatr(3,4) - theMatr(3,3)*theMatr(1,4);
    const Standard_Real aDet3 = theMatr(2,3)*theMatr(3,4) - theMatr(3,3)*theMatr(2,4);
    const Standard_Real anAbsD1 = Abs(aDet1);
    const Standard_Real anAbsD2 = Abs(aDet2);
    const Standard_Real anAbsD3 = Abs(aDet3);

    if(anAbsD1 >= anAbsD2)
    {
      if(anAbsD1 >= anAbsD3)
      {
        //Det1
        if(anAbsD1 <= aNulValue)
          return isSuccess;

        theDeltaU1Found = Abs(aFree(1)*theMatr(2,4) - aFree(2)*theMatr(1,4))/anAbsD1;
        isSuccess = Standard_True;
      }
      else
      {
        //Det3
        if(anAbsD3 <= aNulValue)
          return isSuccess;

        theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
        isSuccess = Standard_True;
      }
    }
    else
    {
      if(anAbsD2 >= anAbsD3)
      {
        //Det2
        if(anAbsD2 <= aNulValue)
          return isSuccess;

        theDeltaU1Found = Abs(aFree(1)*theMatr(3,4) - aFree(3)*theMatr(1,4))/anAbsD2;
        isSuccess = Standard_True;
      }
      else
      {
        //Det3
        if(anAbsD3 <= aNulValue)
          return isSuccess;

        theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
        isSuccess = Standard_True;
      }
    }
  }

  return isSuccess;
}

//=======================================================================
//function : ProcessBounds
//purpose  : 
//=======================================================================
void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne courante
		   const IntPatch_SequenceOfLine& slin,
		   const IntSurf_Quadric& Quad1,
		   const IntSurf_Quadric& Quad2,
		   Standard_Boolean& procf,
		   const gp_Pnt& ptf,                     //-- Debut Ligne Courante
		   const Standard_Real first,             //-- Paramf
		   Standard_Boolean& procl,
		   const gp_Pnt& ptl,                     //-- Fin Ligne courante
		   const Standard_Real last,              //-- Paraml
		   Standard_Boolean& Multpoint,
		   const Standard_Real Tol)
{  
  Standard_Integer j,k;
  Standard_Real U1,V1,U2,V2;
  IntPatch_Point ptsol;
  Standard_Real d;
  
  if (procf && procl) {
    j = slin.Length() + 1;
  }
  else {
    j = 1;
  }


  //-- On prend les lignes deja enregistrees

  while (j <= slin.Length()) {  
    if(slin.Value(j)->ArcType() == IntPatch_Analytic) { 
      const Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(j));
      k = 1;

      //-- On prend les vertex des lignes deja enregistrees

      while (k <= aligold->NbVertex()) {
	ptsol = aligold->Vertex(k);            
	if (!procf) {
	  d=ptf.Distance(ptsol.Value());
	  if (d <= Tol) {
            ptsol.SetTolerance(Tol);
	    if (!ptsol.IsMultiple()) {
	      //-- le point ptsol (de aligold) est declare multiple sur aligold
	      Multpoint = Standard_True;
	      ptsol.SetMultiple(Standard_True);
	      aligold->Replace(k,ptsol);
	    }
	    ptsol.SetParameter(first);
	    alig->AddVertex(ptsol);
	    alig->SetFirstPoint(alig->NbVertex());
	    procf = Standard_True;

	    //-- On restore le point avec son parametre sur aligold
	    ptsol = aligold->Vertex(k); 
                                        
	  }
	}
	if (!procl) {
	  if (ptl.Distance(ptsol.Value()) <= Tol) {
            ptsol.SetTolerance(Tol);
	    if (!ptsol.IsMultiple()) {
	      Multpoint = Standard_True;
	      ptsol.SetMultiple(Standard_True);
	      aligold->Replace(k,ptsol);
	    }
	    ptsol.SetParameter(last);
	    alig->AddVertex(ptsol);
	    alig->SetLastPoint(alig->NbVertex());
	    procl = Standard_True;

	    //-- On restore le point avec son parametre sur aligold
	    ptsol = aligold->Vertex(k); 
             
	  }
	}
	if (procf && procl) {
	  k = aligold->NbVertex()+1;
	}
	else {
	  k = k+1;
	}
      }
      if (procf && procl) {
	j = slin.Length()+1;
      }
      else {
	j = j+1;
      }
    }
  }
  
  ptsol.SetTolerance(Tol);
  if (!procf && !procl) {
    Quad1.Parameters(ptf,U1,V1);
    Quad2.Parameters(ptf,U2,V2);
    ptsol.SetValue(ptf,Tol,Standard_False);
    ptsol.SetParameters(U1,V1,U2,V2);
    ptsol.SetParameter(first);
    if (ptf.Distance(ptl) <= Tol) {
      ptsol.SetMultiple(Standard_True); // a voir
      Multpoint = Standard_True;        // a voir de meme
      alig->AddVertex(ptsol);
      alig->SetFirstPoint(alig->NbVertex());
      
      ptsol.SetParameter(last);
      alig->AddVertex(ptsol);
      alig->SetLastPoint(alig->NbVertex());
    }
    else { 
      alig->AddVertex(ptsol);
      alig->SetFirstPoint(alig->NbVertex());
      Quad1.Parameters(ptl,U1,V1);
      Quad2.Parameters(ptl,U2,V2);
      ptsol.SetValue(ptl,Tol,Standard_False);
      ptsol.SetParameters(U1,V1,U2,V2);
      ptsol.SetParameter(last);
      alig->AddVertex(ptsol);
      alig->SetLastPoint(alig->NbVertex());
    }
  }
  else if (!procf) {
    Quad1.Parameters(ptf,U1,V1);
    Quad2.Parameters(ptf,U2,V2);
    ptsol.SetValue(ptf,Tol,Standard_False);
    ptsol.SetParameters(U1,V1,U2,V2);
    ptsol.SetParameter(first);
    alig->AddVertex(ptsol);
    alig->SetFirstPoint(alig->NbVertex());
  }
  else if (!procl) {
    Quad1.Parameters(ptl,U1,V1);
    Quad2.Parameters(ptl,U2,V2);
    ptsol.SetValue(ptl,Tol,Standard_False);
    ptsol.SetParameters(U1,V1,U2,V2);
    ptsol.SetParameter(last);
    alig->AddVertex(ptsol);
    alig->SetLastPoint(alig->NbVertex());
  }
}

//=======================================================================
//function : CyCyAnalyticalIntersect
//purpose  : Checks if intersection curve is analytical (line, circle, ellipse)
//            and returns these curves.
//=======================================================================
Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1,
                                          const IntSurf_Quadric& Quad2,
                                          const IntAna_QuadQuadGeo& theInter,
                                          const Standard_Real Tol,
                                          Standard_Boolean& Empty,
                                          Standard_Boolean& Same,
                                          Standard_Boolean& Multpoint,
                                          IntPatch_SequenceOfLine& slin,
                                          IntPatch_SequenceOfPoint& spnt)
{
  IntPatch_Point ptsol;
  
  Standard_Integer i;

  IntSurf_TypeTrans trans1,trans2;
  IntAna_ResultType typint;

  gp_Elips elipsol;
  gp_Lin linsol;

  gp_Cylinder Cy1(Quad1.Cylinder());
  gp_Cylinder Cy2(Quad2.Cylinder());

  typint = theInter.TypeInter();
  Standard_Integer NbSol = theInter.NbSolutions();
  Empty = Standard_False;
  Same  = Standard_False;

  switch (typint)
      {
  case IntAna_Empty:
    {
      Empty = Standard_True;
      }
    break;

  case IntAna_Same:
      {
      Same  = Standard_True;
      }
    break;

  case IntAna_Point:
    {
      gp_Pnt psol(theInter.Point(1));
      ptsol.SetValue(psol,Tol,Standard_True);

      Standard_Real U1,V1,U2,V2;
      Quad1.Parameters(psol, U1, V1);
      Quad2.Parameters(psol, U2, V2);

      ptsol.SetParameters(U1,V1,U2,V2);
      spnt.Append(ptsol);
    }
    break;

  case IntAna_Line:
    {
      gp_Pnt ptref;
      if (NbSol == 1)
      { // Cylinders are tangent to each other by line
        linsol = theInter.Line(1);
        ptref = linsol.Location();
        
        //Radius-vectors
        gp_Dir crb1(gp_Vec(ptref,Cy1.Location()));
        gp_Dir crb2(gp_Vec(ptref,Cy2.Location()));

        //outer normal lines
        gp_Vec norm1(Quad1.Normale(ptref));
        gp_Vec norm2(Quad2.Normale(ptref));
        IntSurf_Situation situcyl1;
        IntSurf_Situation situcyl2;

        if (crb1.Dot(crb2) < 0.)
        { // centre de courbures "opposes"
            //ATTENTION!!! 
            //        Normal and Radius-vector of the 1st(!) cylinder
            //        is used for judging what the situation of the 2nd(!)
            //        cylinder is.

          if (norm1.Dot(crb1) > 0.)
          {
            situcyl2 = IntSurf_Inside;
          }
          else
          {
            situcyl2 = IntSurf_Outside;
          }

          if (norm2.Dot(crb2) > 0.)
          {
            situcyl1 = IntSurf_Inside;
          }
          else
          {
            situcyl1 = IntSurf_Outside;
          }
        }
        else
          {
          if (Cy1.Radius() < Cy2.Radius())
          {
            if (norm1.Dot(crb1) > 0.)
            {
              situcyl2 = IntSurf_Inside;
            }
            else
            {
              situcyl2 = IntSurf_Outside;
            }

            if (norm2.Dot(crb2) > 0.)
            {
              situcyl1 = IntSurf_Outside;
            }
            else
            {
              situcyl1 = IntSurf_Inside;
            }
          }
          else
          {
            if (norm1.Dot(crb1) > 0.)
            {
              situcyl2 = IntSurf_Outside;
            }
            else
            {
              situcyl2 = IntSurf_Inside;
            }

            if (norm2.Dot(crb2) > 0.)
            {
              situcyl1 = IntSurf_Inside;
            }
            else
            {
              situcyl1 = IntSurf_Outside;
            }
          }
        }

        Handle(IntPatch_GLine) glig =  new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2);
        slin.Append(glig);
      }
      else
      {
        for (i=1; i <= NbSol; i++)
        {
          linsol = theInter.Line(i);
          ptref = linsol.Location();
          gp_Vec lsd = linsol.Direction();

          //Theoretically, qwe = +/- 1.0.
          Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.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(linsol, Standard_False,trans1,trans2);
          slin.Append(glig);
        }
      }
    }
    break;
    
  case IntAna_Ellipse:
    {
      gp_Vec Tgt;
      gp_Pnt ptref;
      IntPatch_Point pmult1, pmult2;

      elipsol = theInter.Ellipse(1);

      gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
      gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));

      Multpoint = Standard_True;
      pmult1.SetValue(pttang1,Tol,Standard_True);
      pmult2.SetValue(pttang2,Tol,Standard_True);
      pmult1.SetMultiple(Standard_True);
      pmult2.SetMultiple(Standard_True);

      Standard_Real oU1,oV1,oU2,oV2;
      Quad1.Parameters(pttang1, oU1, oV1);
      Quad2.Parameters(pttang1, oU2, oV2);

      pmult1.SetParameters(oU1,oV1,oU2,oV2);
      Quad1.Parameters(pttang2,oU1,oV1); 
      Quad2.Parameters(pttang2,oU2,oV2);

      pmult2.SetParameters(oU1,oV1,oU2,oV2);

      // on traite la premiere ellipse
        
      //-- Calcul de la Transition de la ligne 
      ElCLib::D1(0.,elipsol,ptref,Tgt);

      //Theoretically, qwe = +/- |Tgt|.
      Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.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; 
      }

      //-- Transition calculee au point 0 -> Trans2 , Trans1 
      //-- car ici, on devarit calculer en PI
      Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans2,trans1);
      //
      {
        Standard_Real aU1, aV1, aU2, aV2;
        IntPatch_Point aIP;
        gp_Pnt aP (ElCLib::Value(0., elipsol));
        //
        aIP.SetValue(aP,Tol,Standard_False);
        aIP.SetMultiple(Standard_False);
        //
        Quad1.Parameters(aP, aU1, aV1); 
        Quad2.Parameters(aP, aU2, aV2);

        aIP.SetParameters(aU1, aV1, aU2, aV2);
        //
        aIP.SetParameter(0.);
        glig->AddVertex(aIP);
        glig->SetFirstPoint(1);
        //
        aIP.SetParameter(2.*M_PI);
        glig->AddVertex(aIP);
        glig->SetLastPoint(2);
      }
      //
      pmult1.SetParameter(0.5*M_PI);
      glig->AddVertex(pmult1);
      //
      pmult2.SetParameter(1.5*M_PI);
      glig->AddVertex(pmult2);
     
      //
      slin.Append(glig);
      
      //-- Transitions calculee au point 0    OK
      //
      // on traite la deuxieme ellipse
      elipsol = theInter.Ellipse(2);

      Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
      Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
      Standard_Real parampourtransition = 0.0;
      if (param1 < param2)
      {
        pmult1.SetParameter(0.5*M_PI);
        pmult2.SetParameter(1.5*M_PI);
        parampourtransition = M_PI;
      }
      else {
        pmult1.SetParameter(1.5*M_PI);
        pmult2.SetParameter(0.5*M_PI);
        parampourtransition = 0.0;
      }
      
      //-- Calcul des transitions de ligne pour la premiere ligne 
      ElCLib::D1(parampourtransition,elipsol,ptref,Tgt);      

      //Theoretically, qwe = +/- |Tgt|.
      qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.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;
      }

      //-- La transition a ete calculee sur un point de cette ligne 
      glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
      //
      {
        Standard_Real aU1, aV1, aU2, aV2;
        IntPatch_Point aIP;
        gp_Pnt aP (ElCLib::Value(0., elipsol));
        //
        aIP.SetValue(aP,Tol,Standard_False);
        aIP.SetMultiple(Standard_False);
        //

        Quad1.Parameters(aP, aU1, aV1);
        Quad2.Parameters(aP, aU2, aV2);

        aIP.SetParameters(aU1, aV1, aU2, aV2);
        //
        aIP.SetParameter(0.);
        glig->AddVertex(aIP);
        glig->SetFirstPoint(1);
        //
        aIP.SetParameter(2.*M_PI);
        glig->AddVertex(aIP);
        glig->SetLastPoint(2);
      }
      //
      glig->AddVertex(pmult1);
      glig->AddVertex(pmult2);
      //
      slin.Append(glig);
    }
    break;

  case IntAna_Parabola:
  case IntAna_Hyperbola:
    throw Standard_Failure("IntCyCy(): Wrong intersection type!");

  case IntAna_Circle:
    // Circle is useful when we will work with trimmed surfaces
    // (two cylinders can be tangent by their basises, e.g. circle)
  case IntAna_NoGeometricSolution:
  default:
    return Standard_False;
  }

  return Standard_True;
}

//=======================================================================
//function : ShortCosForm
//purpose  : Represents theCosFactor*cosA+theSinFactor*sinA as
//            theCoeff*cos(A-theAngle) if it is possibly (all angles are
//            in radians).
//=======================================================================
static void ShortCosForm( const Standard_Real theCosFactor,
                          const Standard_Real theSinFactor,
                          Standard_Real& theCoeff,
                          Standard_Real& theAngle)
{
  theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
  theAngle = 0.0;
  if(IsEqual(theCoeff, 0.0))
  {
    theAngle = 0.0;
    return;
  }

  theAngle = acos(Abs(theCosFactor/theCoeff));

  if(theSinFactor > 0.0)
  {
    if(IsEqual(theCosFactor, 0.0))
    {
      theAngle = M_PI/2.0;
    }
    else if(theCosFactor < 0.0)
    {
      theAngle = M_PI-theAngle;
    }
  }
  else if(IsEqual(theSinFactor, 0.0))
  {
    if(theCosFactor < 0.0)
    {
      theAngle = M_PI;
    }
  }
  if(theSinFactor < 0.0)
  {
    if(theCosFactor > 0.0)
    {
      theAngle = 2.0*M_PI-theAngle;
    }
    else if(IsEqual(theCosFactor, 0.0))
    {
      theAngle = 3.0*M_PI/2.0;
    }
    else if(theCosFactor < 0.0)
    {
      theAngle = M_PI+theAngle;
    }
  }
}

//=======================================================================
//function : CylCylMonotonicity
//purpose  : Determines, if U2(U1) function is increasing.
//=======================================================================
Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par,
                                                        const Standard_Integer theWLIndex,
                                                        const stCoeffsValue& theCoeffs,
                                                        const Standard_Real thePeriod,
                                                        Standard_Boolean& theIsIncreasing)
{
  // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
  //Accordingly,
  //Func. U2(X1) = FI2 + X1;
  //Func. X1(X2) = anArccosFactor*X2;
  //Func. X2(X3) = acos(X3);
  //Func. X3(X4) = B*X4 + C;
  //Func. X4(U1) = cos(U1 - FI1).
  //
  //Consequently,
  //U2(X1) is always increasing.
  //X1(X2) is increasing if anArccosFactor > 0.0 and
  //is decreasing otherwise.
  //X2(X3) is always decreasing.
  //Therefore, U2(X3) is decreasing if anArccosFactor > 0.0 and
  //is increasing otherwise.
  //X3(X4) is increasing if B > 0 and is decreasing otherwise.
  //X4(U1) is increasing if U1 - FI1 in [PI, 2*PI) and
  //is decreasing U1 - FI1 in [0, PI) or if (U1 - FI1 == 2PI).
  //After that, we can predict behaviour of U2(U1) function:
  //if it is increasing or decreasing.

  //For "+/-" sign. If isPlus == TRUE, "+" is chosen, otherwise, we choose "-".
  Standard_Boolean isPlus = Standard_False;

  switch(theWLIndex)
  {
  case 0: 
    isPlus = Standard_True;
    break;
  case 1:
    isPlus = Standard_False;
    break;
  default:
    //throw Standard_Failure("Error. Range Error!!!!");
    return Standard_False;
  }

  Standard_Real aU1Temp = theU1par - theCoeffs.mFI1;
  InscribePoint(0, thePeriod, aU1Temp, 0.0, thePeriod, Standard_False);

  theIsIncreasing = Standard_True;

  if(((M_PI - aU1Temp) < RealSmall()) && (aU1Temp < thePeriod))
  {
    theIsIncreasing = Standard_False;
  }

  if(theCoeffs.mB < 0.0)
  {
    theIsIncreasing = !theIsIncreasing;
  }

  if(!isPlus)
  {
    theIsIncreasing = !theIsIncreasing;
  }  

  return Standard_True;
}

//=======================================================================
//function : CylCylComputeParameters
//purpose  : Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
//            estimates the tolerance of U2-computing (estimation result is
//            assigned to *theDelta value).
//=======================================================================
Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
                                                const Standard_Integer theWLIndex,
                                                const stCoeffsValue& theCoeffs,
                                                Standard_Real& theU2,
                                                Standard_Real* const theDelta)
{
  //This formula is got from some experience and can be changed.
  const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue);
  const Standard_Real aTol = 1.0 - aTol0;

  if(theWLIndex < 0 || theWLIndex > 1)
    return Standard_False;

  const Standard_Real aSign = theWLIndex ? -1.0 : 1.0;

  Standard_Real anArg = cos(theU1par - theCoeffs.mFI1);
  anArg = theCoeffs.mB*anArg + theCoeffs.mC;

  if(anArg >= aTol)
  {
    if(theDelta)
      *theDelta = 0.0;

    anArg = 1.0;
  }
  else if(anArg <= -aTol)
  {
    if(theDelta)
      *theDelta = 0.0;

    anArg = -1.0;
  }
  else if(theDelta)
  {
    //There is a case, when
    //  const double aPar = 0.99999999999721167;
    //  const double aFI2 = 2.3593296083566181e-006;

    //Then
    //  aPar - cos(aFI2) == -5.10703e-015 ==> cos(aFI2) == aPar. 
    //Theoretically, in this case
    //  aFI2 == +/- acos(aPar).
    //However,
    //  acos(aPar) - aFI2 == 2.16362e-009.
    //Error is quite big.

    //This error should be estimated. Let use following way, which is based
    //on expanding to Taylor series.

    //  acos(p)-acos(p+x) = x/sqrt(1-p*p).

    //If p == (1-d) (when p > 0) or p == (-1+d) (when p < 0) then
    //  acos(p)-acos(p+x) = x/sqrt(d*(2-d)).

    //Here always aTol0 <= d <= 1. Max(x) is considered (!) to be equal to aTol0.
    //In this case
    //  8*aTol0 <= acos(p)-acos(p+x) <= sqrt(2/(2-aTol0)-1),
    //                                              because 0 < aTol0 < 1.
    //E.g. when aTol0 = 1.0e-11,
    //   8.0e-11 <= acos(p)-acos(p+x) < 2.24e-6.

    const Standard_Real aDelta = Min(1.0-anArg, 1.0+anArg);
    Standard_DivideByZero_Raise_if((aDelta*aDelta < RealSmall()) || (aDelta >= 2.0), 
          "IntPatch_ImpImpIntersection_4.gxx, CylCylComputeParameters()");
    *theDelta = aTol0/sqrt(aDelta*(2.0-aDelta));
  }

  theU2 = acos(anArg);
  theU2 = theCoeffs.mFI2 + aSign*theU2;

  return Standard_True;
}

//=======================================================================
//function : CylCylComputeParameters
//purpose  : Computes V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
//=======================================================================
Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1,
                                                const Standard_Real theU2,
                                                const stCoeffsValue& theCoeffs,
                                                Standard_Real& theV1,
                                                Standard_Real& theV2)
{
  theV1 = theCoeffs.mK21 * sin(theU2) + 
          theCoeffs.mK11 * sin(theU1) +
          theCoeffs.mL21 * cos(theU2) +
          theCoeffs.mL11 * cos(theU1) + theCoeffs.mM1;

  theV2 = theCoeffs.mK22 * sin(theU2) +
          theCoeffs.mK12 * sin(theU1) +
          theCoeffs.mL22 * cos(theU2) +
          theCoeffs.mL12 * cos(theU1) + theCoeffs.mM2;

  return Standard_True;
}

//=======================================================================
//function : CylCylComputeParameters
//purpose  : Computes U2 (U-parameter of the 2nd cylinder),
//            V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
//=======================================================================
Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
                                                const Standard_Integer theWLIndex,
                                                const stCoeffsValue& theCoeffs,
                                                Standard_Real& theU2,
                                                Standard_Real& theV1,
                                                Standard_Real& theV2)
{
  if(!CylCylComputeParameters(theU1par, theWLIndex, theCoeffs, theU2))
    return Standard_False;

  if(!CylCylComputeParameters(theU1par, theU2, theCoeffs, theV1, theV2))
    return Standard_False;

  return Standard_True;
}

//=======================================================================
//function : SearchOnVBounds
//purpose  : 
//=======================================================================
Standard_Boolean WorkWithBoundaries::
                        SearchOnVBounds(const SearchBoundType theSBType,
                                        const Standard_Real theVzad,
                                        const Standard_Real theVInit,
                                        const Standard_Real theInitU2,
                                        const Standard_Real theInitMainVar,
                                        Standard_Real& theMainVariableValue) const
{
  const Standard_Integer aNbDim = 3;
  const Standard_Real aMaxError = 4.0*M_PI; // two periods
  
  theMainVariableValue = theInitMainVar;
  const Standard_Real aTol2 = 1.0e-18;
  Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit;
  
  //Structure of aMatr:
  //  C_{1}*U_{1} & C_{2}*U_{2} & C_{3}*V_{*},
  //where C_{1}, C_{2} and C_{3} are math_Vector.
  math_Matrix aMatr(1, aNbDim, 1, aNbDim);

  Standard_Real anError = RealLast();
  Standard_Real anErrorPrev = anError;
  Standard_Integer aNbIter = 0;
  do
  {
    if(++aNbIter > 1000)
      return Standard_False;

    const Standard_Real aSinU1 = sin(aMainVarPrev),
                        aCosU1 = cos(aMainVarPrev),
                        aSinU2 = sin(aU2Prev),
                        aCosU2 = cos(aU2Prev);

    math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev +
                                              myCoeffs.mVecB2) * aSinU2 -
                              (myCoeffs.mVecB2 * aU2Prev -
                                              myCoeffs.mVecA2) * aCosU2 +
                              (myCoeffs.mVecA1 * aMainVarPrev +
                                              myCoeffs.mVecB1) * aSinU1 -
                              (myCoeffs.mVecB1 * aMainVarPrev -
                                              myCoeffs.mVecA1) * aCosU1 +
                                                            myCoeffs.mVecD;

    math_Vector aMSum(1, 3);

    switch(theSBType)
    {
    case SearchV1:
      aMatr.SetCol(3, myCoeffs.mVecC2);
      aMSum = myCoeffs.mVecC1 * theVzad;
      aVecFreeMem -= aMSum;
      aMSum += myCoeffs.mVecC2*anOtherVar;
      break;

    case SearchV2:
      aMatr.SetCol(3, myCoeffs.mVecC1);
      aMSum = myCoeffs.mVecC2 * theVzad;
      aVecFreeMem -= aMSum;
      aMSum += myCoeffs.mVecC1*anOtherVar;
      break;

    default:
      return Standard_False;
    }

    aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1);
    aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2);

    Standard_Real aDetMainSyst = aMatr.Determinant();

    if(Abs(aDetMainSyst) < aNulValue)
    {
      return Standard_False;
    }

    math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr);
    aM1.SetCol(1, aVecFreeMem);
    aM2.SetCol(2, aVecFreeMem);
    aM3.SetCol(3, aVecFreeMem);

    const Standard_Real aDetMainVar = aM1.Determinant();
    const Standard_Real aDetVar1    = aM2.Determinant();
    const Standard_Real aDetVar2    = aM3.Determinant();

    Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;

    if(Abs(aDelta) > aMaxError)
      return Standard_False;

    anError = aDelta*aDelta;
    aMainVarPrev += aDelta;
        
    ///
    aDelta = aDetVar1/aDetMainSyst-aU2Prev;

    if(Abs(aDelta) > aMaxError)
      return Standard_False;

    anError += aDelta*aDelta;
    aU2Prev += aDelta;

    ///
    aDelta = aDetVar2/aDetMainSyst-anOtherVar;
    anError += aDelta*aDelta;
    anOtherVar += aDelta;

    if(anError > anErrorPrev)
    {//Method diverges. Keep the best result
      const Standard_Real aSinU1Last = sin(aMainVarPrev),
                          aCosU1Last = cos(aMainVarPrev),
                          aSinU2Last = sin(aU2Prev),
                          aCosU2Last = cos(aU2Prev);
      aMSum -= (myCoeffs.mVecA1*aCosU1Last + 
                myCoeffs.mVecB1*aSinU1Last +
                myCoeffs.mVecA2*aCosU2Last +
                myCoeffs.mVecB2*aSinU2Last +
                myCoeffs.mVecD);
      const Standard_Real aSQNorm = aMSum.Norm2();
      return (aSQNorm < aTol2);
    }
    else
    {
      theMainVariableValue = aMainVarPrev;
    }

    anErrorPrev = anError;
  }
  while(anError > aTol2);

  theMainVariableValue = aMainVarPrev;

  return Standard_True;
}

//=======================================================================
//function : InscribePoint
//purpose  : If theFlForce==TRUE theUGiven will be changed forcefully
//            even if theUGiven is already inscibed in the boundary
//            (if it is possible; i.e. if new theUGiven is inscribed
//            in the boundary, too).
//=======================================================================
Standard_Boolean InscribePoint( const Standard_Real theUfTarget,
                                const Standard_Real theUlTarget,
                                Standard_Real& theUGiven,
                                const Standard_Real theTol2D,
                                const Standard_Real thePeriod,
                                const Standard_Boolean theFlForce)
{
  if(Precision::IsInfinite(theUGiven))
  {
    return Standard_False;
  }

  if((theUfTarget - theUGiven <= theTol2D) &&
              (theUGiven - theUlTarget <= theTol2D))
  {//it has already inscribed

    /*
             Utf    U      Utl
              +     *       +
    */
    
    if(theFlForce)
    {
      Standard_Real anUtemp = theUGiven + thePeriod;
      if((theUfTarget - anUtemp <= theTol2D) &&
                (anUtemp - theUlTarget <= theTol2D))
      {
        theUGiven = anUtemp;
        return Standard_True;
      }
      
      anUtemp = theUGiven - thePeriod;
      if((theUfTarget - anUtemp <= theTol2D) &&
                (anUtemp - theUlTarget <= theTol2D))
      {
        theUGiven = anUtemp;
      }
    }

    return Standard_True;
  }

  const Standard_Real aUf = theUfTarget - theTol2D;
  const Standard_Real aUl = aUf + thePeriod;

  theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl);
  
  return ((theUfTarget - theUGiven <= theTol2D) &&
          (theUGiven - theUlTarget <= theTol2D));
}

//=======================================================================
//function : InscribeInterval
//purpose  : Shifts theRange in order to make at least one of its 
//            boundary in the range [theUfTarget, theUlTarget]
//=======================================================================
static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
                                         const Standard_Real theUlTarget,
                                         Bnd_Range &theRange,
                                         const Standard_Real theTol2D,
                                         const Standard_Real thePeriod)
{
  Standard_Real anUpar = 0.0;
  if (!theRange.GetMin(anUpar))
  {
    return Standard_False;
  }

  const Standard_Real aDelta = theRange.Delta();
  if(InscribePoint(theUfTarget, theUlTarget, anUpar, 
          theTol2D, thePeriod, (Abs(theUlTarget-anUpar) < theTol2D)))
  {
    theRange.SetVoid();
    theRange.Add(anUpar);
    theRange.Add(anUpar + aDelta);
  }
  else 
  {
    if (!theRange.GetMax (anUpar))
    {
      return Standard_False;
    }

    if(InscribePoint(theUfTarget, theUlTarget, anUpar,
          theTol2D, thePeriod, (Abs(theUfTarget-anUpar) < theTol2D)))
    {
      theRange.SetVoid();
      theRange.Add(anUpar);
      theRange.Add(anUpar - aDelta);
    }
    else
    {
      return Standard_False;
    }
  }

  return Standard_True;
}

//=======================================================================
//function : ExcludeNearElements
//purpose  : Checks if theArr contains two almost equal elements.
//            If it is true then one of equal elements will be excluded
//            (made infinite).
//           Returns TRUE if at least one element of theArr has been changed.
//ATTENTION!!!
//           1. Every not infinite element of theArr is considered to be 
//            in [0, T] interval (where T is the period);
//           2. theArr must be sorted in ascending order.
//=======================================================================
static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
                                            const Standard_Integer theNOfMembers,
                                            const Standard_Real theUSurf1f,
                                            const Standard_Real theUSurf1l,
                                            const Standard_Real theTol)
{
  Standard_Boolean aRetVal = Standard_False;
  for(Standard_Integer i = 1; i < theNOfMembers; i++)
  {
    Standard_Real &anA = theArr[i],
                  &anB = theArr[i-1];

    //Here, anA >= anB

    if(Precision::IsInfinite(anA))
      break;

    if((anA-anB) < theTol)
    {
      if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l)) 
      anA = (anA + anB)/2.0;
      else
        anA = anB;

      //Make this element infinite an forget it
      //(we will not use it in next iterations).
      anB = Precision::Infinite();
      aRetVal = Standard_True;
    }
  }

  return aRetVal;
}

//=======================================================================
//function : AddPointIntoWL
//purpose  : Surf1 is a surface, whose U-par is variable.
//           If theFlBefore == TRUE then we enable the U1-parameter
//            of the added point to be less than U1-parameter of
//           previously added point (in general U1-parameter is
//           always increased; therefore, this situation is abnormal).
//           If theOnlyCheck==TRUE then no point will be added to theLine.
//=======================================================================
static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
                                        const IntSurf_Quadric& theQuad2,
                                        const ComputationMethods::stCoeffsValue& theCoeffs,
                                        const Standard_Boolean isTheReverse,
                                        const Standard_Boolean isThePrecise,
                                        const gp_Pnt2d& thePntOnSurf1,
                                        const gp_Pnt2d& thePntOnSurf2,
                                        const Standard_Real theUfSurf1,
                                        const Standard_Real theUlSurf1,
                                        const Standard_Real theUfSurf2,
                                        const Standard_Real theUlSurf2,
                                        const Standard_Real theVfSurf1,
                                        const Standard_Real theVlSurf1,
                                        const Standard_Real theVfSurf2,
                                        const Standard_Real theVlSurf2,
                                        const Standard_Real thePeriodOfSurf1,
                                        const Handle(IntSurf_LineOn2S)& theLine,
                                        const Standard_Integer theWLIndex,
                                        const Standard_Real theTol3D,
                                        const Standard_Real theTol2D,
                                        const Standard_Boolean theFlBefore = Standard_False,
                                        const Standard_Boolean theOnlyCheck = Standard_False)
{
  //Check if the point is in the domain or can be inscribed in the domain after adjusting.
  
  gp_Pnt  aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
          aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));

  Standard_Real aU1par = thePntOnSurf1.X();

  // aU1par always increases. Therefore, we must reduce its
  // value in order to continue creation of WLine.
  if(!InscribePoint(theUfSurf1, theUlSurf1, aU1par, theTol2D,
                  thePeriodOfSurf1, aU1par > 0.5*(theUfSurf1+theUlSurf1)))
    return Standard_False;

  if ((theLine->NbPoints() > 0) &&
      ((theUlSurf1 - theUfSurf1) >= (thePeriodOfSurf1 - theTol2D)) &&      
      (((aU1par + thePeriodOfSurf1 - theUlSurf1) <= theTol2D) ||
       ((aU1par - thePeriodOfSurf1 - theUfSurf1) >= theTol2D)))
  {
    // aU1par can be adjusted to both theUlSurf1 and theUfSurf1
    // with equal possibilities. This code fragment allows choosing
    // correct parameter from these two variants.

    Standard_Real aU1 = 0.0, aV1 = 0.0;
    if (isTheReverse)
    {
      theLine->Value(theLine->NbPoints()).ParametersOnS2(aU1, aV1);
    }
    else
    {
      theLine->Value(theLine->NbPoints()).ParametersOnS1(aU1, aV1);
    }

    const Standard_Real aDelta = aU1 - aU1par;
    if (2.0*Abs(aDelta) > thePeriodOfSurf1)
    {
      aU1par += Sign(thePeriodOfSurf1, aDelta);
    }
  }

  Standard_Real aU2par = thePntOnSurf2.X();
  if(!InscribePoint(theUfSurf2, theUlSurf2, aU2par, theTol2D,
                                    thePeriodOfSurf1, Standard_False))
    return Standard_False;

  Standard_Real aV1par = thePntOnSurf1.Y();
  if((aV1par - theVlSurf1 > theTol2D) || (theVfSurf1 - aV1par > theTol2D))
    return Standard_False;

  Standard_Real aV2par = thePntOnSurf2.Y();
  if((aV2par -  theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
    return Standard_False;
  
  //Get intersection point and add it in the WL
  IntSurf_PntOn2S aPnt;
  
  if(isTheReverse)
  {
    aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
                        aU2par, aV2par,
                        aU1par, aV1par);
  }
  else
  {
    aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
                        aU1par, aV1par,
                        aU2par, aV2par);
  }

  Standard_Integer aNbPnts = theLine->NbPoints();
  if(aNbPnts > 0)
  {
    Standard_Real aUl = 0.0, aVl = 0.0;
    const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts);
    if(isTheReverse)
      aPlast.ParametersOnS2(aUl, aVl);
    else
      aPlast.ParametersOnS1(aUl, aVl);

    if(!theFlBefore && (aU1par <= aUl))
    {
      //Parameter value must be increased if theFlBefore == FALSE.

      aU1par += thePeriodOfSurf1;

      //The condition is as same as in
      //InscribePoint(...) function
      if((theUfSurf1 - aU1par > theTol2D) ||
         (aU1par - theUlSurf1 > theTol2D))
      {
        //New aU1par is out of target interval.
        //Go back to old value.

        return Standard_False;
      }
    }

    if (theOnlyCheck)
      return Standard_True;

    //theTol2D is minimal step along parameter changed. 
    //Therefore, if we apply this minimal step two 
    //neighbour points will be always "same". Consequently,
    //we should reduce tolerance for IsSame checking.
    const Standard_Real aDTol = 1.0-Epsilon(1.0);
    if(aPnt.IsSame(aPlast, theTol3D*aDTol, theTol2D*aDTol))
    {
      theLine->RemovePoint(aNbPnts);
    }
  }

  if (theOnlyCheck)
    return Standard_True;

  theLine->Add(aPnt);

  if(!isThePrecise)
    return Standard_True;

  //Try to precise existing WLine
  aNbPnts = theLine->NbPoints();
  if(aNbPnts >= 3)
  {
    Standard_Real aU1 = 0.0, aU2 = 0.0, aU3 = 0.0, aV = 0.0;
    if(isTheReverse)
    {
      theLine->Value(aNbPnts).ParametersOnS2(aU3, aV);
      theLine->Value(aNbPnts-1).ParametersOnS2(aU2, aV);
      theLine->Value(aNbPnts-2).ParametersOnS2(aU1, aV);
    }
    else
    {
      theLine->Value(aNbPnts).ParametersOnS1(aU3, aV);
      theLine->Value(aNbPnts-1).ParametersOnS1(aU2, aV);
      theLine->Value(aNbPnts-2).ParametersOnS1(aU1, aV);
    }

    const Standard_Real aStepPrev = aU2-aU1;
    const Standard_Real aStep = aU3-aU2;

    const Standard_Integer aDeltaStep = RealToInt(aStepPrev/aStep);

    if((1 < aDeltaStep) && (aDeltaStep < 2000))
    {
      //Add new points in case of non-uniform distribution of existing points
      SeekAdditionalPoints( theQuad1, theQuad2, theLine, 
                            theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
                            aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse);
    }
  }

  return Standard_True;
}

//=======================================================================
//function : AddBoundaryPoint
//purpose  : Find intersection point on V-boundary.
//=======================================================================
void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
                                          const Standard_Real theU1,
                                          const Standard_Real theU1Min,
                                          const Standard_Real theU2,
                                          const Standard_Real theV1,
                                          const Standard_Real theV1Prev,
                                          const Standard_Real theV2,
                                          const Standard_Real theV2Prev,
                                          const Standard_Integer theWLIndex,
                                          const Standard_Boolean theFlForce,
                                          Standard_Boolean& isTheFound1,
                                          Standard_Boolean& isTheFound2) const
{
  Standard_Real aUSurf1f = 0.0, //const
                aUSurf1l = 0.0,
                aVSurf1f = 0.0,
                aVSurf1l = 0.0;
  Standard_Real aUSurf2f = 0.0, //const
                aUSurf2l = 0.0,
                aVSurf2f = 0.0,
                aVSurf2l = 0.0;

  myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
  myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
  
  const Standard_Integer aSize = 4;
  const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l};

  StPInfo aUVPoint[aSize];

  for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2)
  {
    const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2,
                        aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev;

    const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2;

    for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++)
    {
      const Standard_Integer anIndex = anIDSurf+anIDBound;

      aUVPoint[anIndex].mySurfID = anIDSurf;

      if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) &&
          ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0))
      {
        continue;
      }

      //Segment [aVf, aVl] intersects at least one V-boundary (first or last)
      // (in general, case is possible, when aVf > aVl).

      // Precise intersection point
      const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex],
                                                    (anIDSurf == 0) ? theV2 : theV1,
                                                    theU2, theU1,
                                                    aUVPoint[anIndex].myU1);

      // aUVPoint[anIndex].myU1 is considered to be nearer to theU1 than
      // to theU1+/-Period
      if (!aRes || (aUVPoint[anIndex].myU1 >= theU1) ||
                              (aUVPoint[anIndex].myU1 < theU1Min))
      {
        //Intersection point is not found or out of the domain
        aUVPoint[anIndex].myU1 = RealLast();
        continue;
      }
      else
      {
        //intersection point is found

        Standard_Real &aU1 = aUVPoint[anIndex].myU1,
                      &aU2 = aUVPoint[anIndex].myU2,
                      &aV1 = aUVPoint[anIndex].myV1,
                      &aV2 = aUVPoint[anIndex].myV2;
        aU2 = theU2;
        aV1 = theV1;
        aV2 = theV2;

        if(!ComputationMethods::
                  CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2))
        {
          // Found point is wrong
          aU1 = RealLast();
          continue;
        }

        //Point on true V-boundary.
        if(aTS == SearchV1)
          aV1 = anArrVzad[anIndex];
        else //if(aTS[anIndex] == SearchV2)
          aV2 = anArrVzad[anIndex];
      }
    }
  }

  // Sort with acceding U1-parameter.
  std::sort(aUVPoint, aUVPoint+aSize);
    
  isTheFound1 = isTheFound2 = Standard_False;

  //Adding found points on boundary in the WLine.
  for(Standard_Integer i = 0; i < aSize; i++)
  {
    if(aUVPoint[i].myU1 == RealLast())
      break;

    if(!AddPointIntoWL(myQuad1, myQuad2, myCoeffs, myIsReverse, Standard_False,
                        gp_Pnt2d(aUVPoint[i].myU1, aUVPoint[i].myV1),
                        gp_Pnt2d(aUVPoint[i].myU2, aUVPoint[i].myV2),
                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, myPeriod,
                        theWL->Curve(), theWLIndex, myTol3D, myTol2D, theFlForce))
    {
      continue;
    }

    if(aUVPoint[i].mySurfID == 0)
    {
      isTheFound1 = Standard_True;
    }
    else
    {
      isTheFound2 = Standard_True;
    }
  }
}

//=======================================================================
//function : SeekAdditionalPoints
//purpose  : Inserts additional intersection points between neighbor points.
//            This action is bone with several iterations. During every iteration,
//          new point is inserted in middle of every interval.
//            The process will be finished if NbPoints >= theMinNbPoints.
//=======================================================================
static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
                                  const IntSurf_Quadric& theQuad2,
                                  const Handle(IntSurf_LineOn2S)& theLine,
                                  const ComputationMethods::stCoeffsValue& theCoeffs,
                                  const Standard_Integer theWLIndex,
                                  const Standard_Integer theMinNbPoints,
                                  const Standard_Integer theStartPointOnLine,
                                  const Standard_Integer theEndPointOnLine,
                                  const Standard_Real theTol2D,
                                  const Standard_Real thePeriodOfSurf2,
                                  const Standard_Boolean isTheReverse)
{
  if(theLine.IsNull())
    return;
  
  Standard_Integer aNbPoints = theEndPointOnLine - theStartPointOnLine + 1;

  Standard_Real aMinDeltaParam = theTol2D;

  {
    Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;

    if(isTheReverse)
    {
      theLine->Value(theStartPointOnLine).ParametersOnS2(u1, v1);
      theLine->Value(theEndPointOnLine).ParametersOnS2(u2, v2);
    }
    else
    {
      theLine->Value(theStartPointOnLine).ParametersOnS1(u1, v1);
      theLine->Value(theEndPointOnLine).ParametersOnS1(u2, v2);
    }
    
    aMinDeltaParam = Max(Abs(u2 - u1)/IntToReal(theMinNbPoints), aMinDeltaParam);
  }

  Standard_Integer aLastPointIndex = theEndPointOnLine;
  Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;

  Standard_Integer aNbPointsPrev = 0;
  do
  {
    aNbPointsPrev = aNbPoints;
    for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
    {
      Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
      Standard_Real U1l = 0.0, V1l = 0.0; //last  point in 1st suraface

      Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface
      Standard_Real U2l = 0.0, V2l = 0.0; //last  point in 2nd suraface

      lp = fp+1;
      
      if(isTheReverse)
      {
        theLine->Value(fp).ParametersOnS2(U1f, V1f);
        theLine->Value(lp).ParametersOnS2(U1l, V1l);

        theLine->Value(fp).ParametersOnS1(U2f, V2f);
        theLine->Value(lp).ParametersOnS1(U2l, V2l);
      }
      else
      {
        theLine->Value(fp).ParametersOnS1(U1f, V1f);
        theLine->Value(lp).ParametersOnS1(U1l, V1l);

        theLine->Value(fp).ParametersOnS2(U2f, V2f);
        theLine->Value(lp).ParametersOnS2(U2l, V2l);
      }

      if(Abs(U1l - U1f) <= aMinDeltaParam)
      {
        //Step is minimal. It is not necessary to divide it.
        continue;
      }

      U1prec = 0.5*(U1f+U1l);
      
      if(!ComputationMethods::
            CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
      {
        continue;
      }

      MinMax(U2f, U2l);
      if(!InscribePoint(U2f, U2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False))
      {
        continue;
      }

      const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec));
      const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));

#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
      std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl;
#endif

      IntSurf_PntOn2S anIP;
      if(isTheReverse)
      {
        anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
      }
      else
      {
        anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
      }

      theLine->InsertBefore(lp, anIP);

      aNbPoints++;
      aLastPointIndex++;
    }

    if(aNbPoints >= theMinNbPoints)
    {
      return;
    }
  } while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev));
}

//=======================================================================
//function : BoundariesComputing
//purpose  : Computes true domain of future intersection curve (allows
//            avoiding excess iterations)
//=======================================================================
Standard_Boolean WorkWithBoundaries::
            BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
                                const Standard_Real thePeriod,
                                Bnd_Range theURange[])
{
  //All comments to this method is related to the comment
  //to ComputationMethods class
  
  //So, we have the equation
  //    cos(U2-FI2)=B*cos(U1-FI1)+C
  //Evidently,
  //    -1 <= B*cos(U1-FI1)+C <= 1

  if (theCoeffs.mB > 0.0)
  {
    // -(1+C)/B <= cos(U1-FI1) <= (1-C)/B

    if (theCoeffs.mB + Abs(theCoeffs.mC) < -1.0)
    {
      //(1-C)/B < -1 or -(1+C)/B > 1  ==> No solution
      
      return Standard_False;
    }
    else if (theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
    {
      //(1-C)/B >= 1 and -(1+C)/B <= -1 ==> U=[0;2*PI]+aFI1
      theURange[0].Add(theCoeffs.mFI1);
      theURange[0].Add(thePeriod + theCoeffs.mFI1);
    }
    else if ((1 + theCoeffs.mC <= theCoeffs.mB) &&
             (theCoeffs.mB <= 1 - theCoeffs.mC))
    {
      //(1-C)/B >= 1 and -(1+C)/B >= -1 ==> 
      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
      //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB)

      Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
      if(anArg > 1.0)
        anArg = 1.0;
      if(anArg < -1.0)
        anArg = -1.0;

      const Standard_Real aDAngle = acos(anArg);
      theURange[0].Add(theCoeffs.mFI1);
      theURange[0].Add(aDAngle + theCoeffs.mFI1);
      theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
      theURange[1].Add(thePeriod + theCoeffs.mFI1);
    }
    else if ((1 - theCoeffs.mC <= theCoeffs.mB) &&
             (theCoeffs.mB <= 1 + theCoeffs.mC))
    {
      //(1-C)/B <= 1 and -(1+C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1
      //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)

      Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
      if(anArg > 1.0)
        anArg = 1.0;
      if(anArg < -1.0)
        anArg = -1.0;

      const Standard_Real aDAngle = acos(anArg);
      theURange[0].Add(aDAngle + theCoeffs.mFI1);
      theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
    }
    else if (theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
    {
      //(1-C)/B <= 1 and -(1+C)/B >= -1 ==>
      //(U=[aDAngle1;aDAngle2]+aFI1) ||
      //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
      //where aDAngle1 = acos((1 - myCoeffs.mC) / myCoeffs.mB),
      //      aDAngle2 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).

      Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB,
                    anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB;
      if(anArg1 > 1.0)
        anArg1 = 1.0;
      if(anArg1 < -1.0)
        anArg1 = -1.0;

      if(anArg2 > 1.0)
        anArg2 = 1.0;
      if(anArg2 < -1.0)
        anArg2 = -1.0;

      const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
      //(U=[aDAngle1;aDAngle2]+aFI1) ||
      //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
      theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
      theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
      theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
      theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
    }
    else
    {
      return Standard_False;
    }
  }
  else if (theCoeffs.mB < 0.0)
  {
    // (1-C)/B <= cos(U1-FI1) <= -(1+C)/B

    if (theCoeffs.mB + Abs(theCoeffs.mC) > 1.0)
    {
      // -(1+C)/B < -1 or (1-C)/B > 1 ==> No solutions
      return Standard_False;
    }
    else if (-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
    {
      //  -(1+C)/B >= 1 and (1-C)/B <= -1 ==> U=[0;2*PI]+aFI1
      theURange[0].Add(theCoeffs.mFI1);
      theURange[0].Add(thePeriod + theCoeffs.mFI1);
    }
    else if ((-theCoeffs.mC - 1 <= theCoeffs.mB) &&
             (theCoeffs.mB <= theCoeffs.mC - 1))
    {
      //  -(1+C)/B >= 1 and (1-C)/B >= -1 ==> 
      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
      //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)

      Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
      if(anArg > 1.0)
        anArg = 1.0;
      if(anArg < -1.0)
        anArg = -1.0;

      const Standard_Real aDAngle = acos(anArg);
      theURange[0].Add(theCoeffs.mFI1);
      theURange[0].Add(aDAngle + theCoeffs.mFI1);
      theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
      theURange[1].Add(thePeriod + theCoeffs.mFI1);
    }
    else if ((theCoeffs.mC - 1 <= theCoeffs.mB) &&
             (theCoeffs.mB <= -theCoeffs.mB - 1))
    {
      //  -(1+C)/B <= 1 and (1-C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1,
      //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).

      Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
      if(anArg > 1.0)
        anArg = 1.0;
      if(anArg < -1.0)
        anArg = -1.0;

      const Standard_Real aDAngle = acos(anArg);
      theURange[0].Add(aDAngle + theCoeffs.mFI1);
      theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
    }
    else if (-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
    {
      //  -(1+C)/B <= 1 and (1-C)/B >= -1 ==>
      //(U=[aDAngle1;aDAngle2]+aFI1) || (U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1),
      //where aDAngle1 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB),
      //      aDAngle2 = acos((1 - myCoeffs.mC) / myCoeffs.mB)

      Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB,
                    anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB;
      if(anArg1 > 1.0)
        anArg1 = 1.0;
      if(anArg1 < -1.0)
        anArg1 = -1.0;

      if(anArg2 > 1.0)
        anArg2 = 1.0;
      if(anArg2 < -1.0)
        anArg2 = -1.0;

      const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
      theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
      theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
      theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
      theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
    }
    else
    {
      return Standard_False;
    }
  }
  else
  {
    return Standard_False;
  }

  return Standard_True;
}

//=======================================================================
//function : CriticalPointsComputing
//purpose  : theNbCritPointsMax contains true number of critical points.
//            It must be initialized correctly before function calling
//=======================================================================
static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs,
                                    const Standard_Real theUSurf1f,
                                    const Standard_Real theUSurf1l,
                                    const Standard_Real theUSurf2f,
                                    const Standard_Real theUSurf2l,
                                    const Standard_Real thePeriod,
                                    const Standard_Real theTol2D,
                                    Standard_Integer& theNbCritPointsMax,
                                    Standard_Real theU1crit[])
{
  //[0...1] - in these points parameter U1 goes through
  //          the seam-edge of the first cylinder.
  //[2...3] - First and last U1 parameter.
  //[4...5] - in these points parameter U2 goes through
  //          the seam-edge of the second cylinder.
  //[6...9] - in these points an intersection line goes through
  //          U-boundaries of the second surface.
  //[10...11] - Boundary of monotonicity interval of U2(U1) function
  //            (see CylCylMonotonicity() function)

  theU1crit[0] = 0.0;
  theU1crit[1] = thePeriod;
  theU1crit[2] = theUSurf1f;
  theU1crit[3] = theUSurf1l;

  const Standard_Real aCOS = cos(theCoeffs.mFI2);
  const Standard_Real aBSB = Abs(theCoeffs.mB);
  if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
  {
    Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
    if(anArg > 1.0)
      anArg = 1.0;
    if(anArg < -1.0)
      anArg = -1.0;

    theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
    theU1crit[5] =  acos(anArg) + theCoeffs.mFI1;
  }

  Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
  Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
  MinMax(aSf, aSl);

  //In accorance with pure mathematic, theU1crit[6] and [8]
  //must be -Precision::Infinite() instead of used +Precision::Infinite()
  theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
                  -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
                  Precision::Infinite();
  theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
                  -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
                  Precision::Infinite();
  theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
                  acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
                  Precision::Infinite();
  theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
                  acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
                  Precision::Infinite();

  theU1crit[10] = theCoeffs.mFI1;
  theU1crit[11] = M_PI+theCoeffs.mFI1;

  //preparative treatment of array. This array must have faled to contain negative
  //infinity number

  for(Standard_Integer i = 0; i < theNbCritPointsMax; i++)
  {
    if(Precision::IsInfinite(theU1crit[i]))
    {
      continue;
    }

    theU1crit[i] = fmod(theU1crit[i], thePeriod);
    if(theU1crit[i] < 0.0)
      theU1crit[i] += thePeriod;
  }

  //Here all not infinite elements of theU1crit are in [0, thePeriod) range

  do
  {
    std::sort(theU1crit, theU1crit + theNbCritPointsMax);
  }
  while(ExcludeNearElements(theU1crit, theNbCritPointsMax,
                            theUSurf1f, theUSurf1l, theTol2D));

  //Here all not infinite elements in theU1crit are different and sorted.

  while(theNbCritPointsMax > 0)
  {
    Standard_Real &anB = theU1crit[theNbCritPointsMax-1];
    if(Precision::IsInfinite(anB))
    {
      theNbCritPointsMax--;
      continue;
    }

    //1st not infinte element is found

    if(theNbCritPointsMax == 1)
      break;

    //Here theNbCritPointsMax > 1

    Standard_Real &anA = theU1crit[0];

    //Compare 1st and last significant elements of theU1crit
    //They may still differs by period.

    if (Abs(anB - anA - thePeriod) < theTol2D)
    {//E.g. anA == 2.0e-17, anB == (thePeriod-1.0e-18)
      anA = (anA + anB - thePeriod)/2.0;
      anB = Precision::Infinite();
      theNbCritPointsMax--;
    }

    //Out of "while(theNbCritPointsMax > 0)" cycle.
    break;
  }

  //Attention! Here theU1crit may be unsorted.
}

//=======================================================================
//function : BoundaryEstimation
//purpose  : Rough estimation of the parameter range.
//=======================================================================
void WorkWithBoundaries::BoundaryEstimation(const gp_Cylinder& theCy1,
                                            const gp_Cylinder& theCy2,
                                            Bnd_Range& theOutBoxS1,
                                            Bnd_Range& theOutBoxS2) const
{
  const gp_Dir &aD1 = theCy1.Axis().Direction(),
               &aD2 = theCy2.Axis().Direction();
  const Standard_Real aR1 = theCy1.Radius(),
                      aR2 = theCy2.Radius();

  //Let consider a parallelogram. Its edges are parallel to aD1 and aD2.
  //Its altitudes are equal to 2*aR1 and 2*aR2 (diameters of the cylinders).
  //In fact, this parallelogram is a projection of the cylinders to the plane
  //created by the intersected axes aD1 and aD2 (if the axes are skewed then
  //one axis can be translated by parallel shifting till intersection).

  const Standard_Real aCosA = aD1.Dot(aD2);
  const Standard_Real aSqSinA = aD1.XYZ().CrossSquareMagnitude(aD2.XYZ());

  //If sine is small then it can be compared with angle.
  if (aSqSinA < Precision::Angular()*Precision::Angular())
    return;

  //Half of delta V. Delta V is a distance between 
  //projections of two opposite parallelogram vertices
  //(joined by the maximal diagonal) to the cylinder axis.
  const Standard_Real aSinA = sqrt(aSqSinA);
  const Standard_Real anAbsCosA = Abs(aCosA);
  const Standard_Real aHDV1 = (aR1 * anAbsCosA + aR2) / aSinA,
                      aHDV2 = (aR2 * anAbsCosA + aR1) / aSinA;

#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
  //The code in this block is created for test only.It is stupidly to create
  //OCCT-test for the method, which will be changed possibly never.
  std::cout << "Reference values: aHDV1 = " << aHDV1 << "; aHDV2 = " << aHDV2 << std::endl;
#endif

  //V-parameters of intersection point of the axes (in case of skewed axes,
  //see comment above).
  Standard_Real aV01 = 0.0, aV02 = 0.0;
  ExtremaLineLine(theCy1.Axis(), theCy2.Axis(), aCosA, aSqSinA, aV01, aV02);

  theOutBoxS1.Add(aV01 - aHDV1);
  theOutBoxS1.Add(aV01 + aHDV1);

  theOutBoxS2.Add(aV02 - aHDV2);
  theOutBoxS2.Add(aV02 + aHDV2);

  theOutBoxS1.Enlarge(Precision::Confusion());
  theOutBoxS2.Enlarge(Precision::Confusion());

  Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
  myUVSurf1.Get(aU1, aV1, aU2, aV2);
  theOutBoxS1.Common(Bnd_Range(aV1, aV2));

  myUVSurf2.Get(aU1, aV1, aU2, aV2);
  theOutBoxS2.Common(Bnd_Range(aV1, aV2));
}

//=======================================================================
//function : CyCyNoGeometric
//purpose  : 
//=======================================================================
static IntPatch_ImpImpIntersection::IntStatus
                    CyCyNoGeometric(const gp_Cylinder &theCyl1,
                                    const gp_Cylinder &theCyl2,
                                    const WorkWithBoundaries &theBW,
                                    Bnd_Range theRange[],
                                    const Standard_Integer theNbOfRanges /*=2*/,
                                    Standard_Boolean& isTheEmpty,
                                    IntPatch_SequenceOfLine& theSlin,
                                    IntPatch_SequenceOfPoint& theSPnt)
{
  Standard_Real aUSurf1f = 0.0, aUSurf1l = 0.0,
                aUSurf2f = 0.0, aUSurf2l = 0.0,
                aVSurf1f = 0.0, aVSurf1l = 0.0,
                aVSurf2f = 0.0, aVSurf2l = 0.0;

  theBW.UVS1().Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
  theBW.UVS2().Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);

  Bnd_Range aRangeS1, aRangeS2;
  theBW.BoundaryEstimation(theCyl1, theCyl2, aRangeS1, aRangeS2);
  if (aRangeS1.IsVoid() || aRangeS2.IsVoid())
    return IntPatch_ImpImpIntersection::IntStatus_OK;

  {
    //Quotation of the message from issue #26894 (author MSV):
    //"We should return fail status from intersector if the result should be an
    //infinite curve of non-analytical type... I propose to define the limit for the
    //extent as the radius divided by 1e+2 and multiplied by 1e+7.
    //Thus, taking into account the number of valuable digits (15), we provide reliable
    //computations with an error not exceeding R/100."
    const Standard_Real aF = 1.0e+5;
    const Standard_Real aMaxV1Range = aF*theCyl1.Radius(), aMaxV2Range = aF*theCyl2.Radius();
    if ((aRangeS1.Delta() > aMaxV1Range) || (aRangeS2.Delta() > aMaxV2Range))
      return IntPatch_ImpImpIntersection::IntStatus_InfiniteSectionCurve;
  }
  //
  Standard_Boolean isGoodIntersection = Standard_False;
  Standard_Real anOptdu = 0.;
  for (;;)
  {
    //Checking parameters of cylinders in order to define "good intersection"
    //"Good intersection" means that axes of cylinders are almost perpendicular and
    // one radius is much smaller than the other and small cylinder is "inside" big one.
    const Standard_Real aToMuchCoeff = 3.;
    const Standard_Real aCritAngle = M_PI / 18.; // 10 degree
    Standard_Real anR1 = theCyl1.Radius();
    Standard_Real anR2 = theCyl2.Radius();
    Standard_Real anRmin = 0., anRmax = 0.;
    //Radius criterion
    if (anR1 > aToMuchCoeff * anR2)
    {
      anRmax = anR1; anRmin = anR2;
    }
    else if (anR2 > aToMuchCoeff * anR1)
    {
      anRmax = anR2; anRmin = anR1;
    }
    else
    {
      break;
    }
    //Angle criterion
    const gp_Ax1& anAx1 = theCyl1.Axis();
    const gp_Ax1& anAx2 = theCyl2.Axis();
    if (!anAx1.IsNormal(anAx2, aCritAngle))
    {
      break;
    }
    //Placement criterion
    gp_Lin anL1(anAx1), anL2(anAx2);
    Standard_Real aDist = anL1.Distance(anL2);
    if (aDist > anRmax / 2.)
    {
      break;
    }

    isGoodIntersection = Standard_True;
    //Estimation of "optimal" du
    //Relative deflection, absolut deflection is Rmin*aDeflection
    Standard_Real aDeflection = 0.001;
    Standard_Integer aNbP = 3;
    if (anRmin * aDeflection > 1.e-3)
    {
      Standard_Real anAngle = 1.0e0 - aDeflection;
      anAngle = 2.0e0 * ACos(anAngle);
      aNbP = (Standard_Integer)(2. * M_PI / anAngle) + 1;
    }
    anOptdu = 2. * M_PI_2 / (Standard_Real)(aNbP - 1);
    break;
  } 
//
  const ComputationMethods::stCoeffsValue &anEquationCoeffs = theBW.SICoeffs();
  const IntSurf_Quadric& aQuad1 = theBW.GetQSurface(1);
  const IntSurf_Quadric& aQuad2 = theBW.GetQSurface(2);
  const Standard_Boolean isReversed = theBW.IsReversed();
  const Standard_Real aTol2D = theBW.Get2dTolerance();
  const Standard_Real aTol3D = theBW.Get3dTolerance();
  const Standard_Real aPeriod = 2.0*M_PI;
  Standard_Integer aNbMaxPoints = 1000;
  Standard_Integer aNbMinPoints = 200;
  Standard_Real du;
  if (isGoodIntersection)
  {
    du = anOptdu;
    aNbMaxPoints = 200;
    aNbMinPoints = 50;
  }
  else
  {
    du = 2. * M_PI / aNbMaxPoints;
  }
  Standard_Integer aNbPts = Min(RealToInt((aUSurf1l - aUSurf1f) / du) + 1, 
                                RealToInt(20.0*theCyl1.Radius()));
  const Standard_Integer aNbPoints = Min(Max(aNbMinPoints, aNbPts), aNbMaxPoints);
  const Standard_Real aStepMin = Max(aTol2D, Precision::PConfusion()), 
    aStepMax = (aUSurf1l - aUSurf1f > M_PI / 100.0) ?
    (aUSurf1l - aUSurf1f) / IntToReal(aNbPoints) : aUSurf1l - aUSurf1f;

 
  //The main idea of the algorithm is to change U1-parameter
  //(U-parameter of theCyl1) from aU1f to aU1l with some step
  //(step is adaptive) and to obtain set of intersection points.

  for (Standard_Integer i = 0; i < theNbOfRanges; i++)
  {
    if (theRange[i].IsVoid())
      continue;

    InscribeInterval(aUSurf1f, aUSurf1l, theRange[i], aTol2D, aPeriod);
  }

  if (theRange[0].Union(theRange[1]))
  {
    // Works only if (theNbOfRanges == 2).
    theRange[1].SetVoid();
  }

  //Critical points are the value of U1-parameter in the points
  //where WL must be decomposed.

  //When U1 goes through critical points its value is set up to this
  //parameter forcefully and the intersection point is added in the line.
  //After that, the WL is broken (next U1 value will be correspond to the new WL).

  //See CriticalPointsComputing(...) function to get detail information about this array.
  const Standard_Integer aNbCritPointsMax = 12;
  Standard_Real anU1crit[aNbCritPointsMax] = { Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite(),
                                               Precision::Infinite() };

  //This list of critical points is not full because it does not contain any points
  //which intersection line goes through V-bounds of cylinders in.
  //They are computed by numerical methods on - line (during algorithm working).
  //The moment is caught, when intersection line goes through V-bounds of any cylinder.

  Standard_Integer aNbCritPoints = aNbCritPointsMax;
  CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
                          aPeriod, aTol2D, aNbCritPoints, anU1crit);

  //Getting Walking-line

  enum WLFStatus
  {
    // No points have been added in WL
    WLFStatus_Absent = 0,
    // WL contains at least one point
    WLFStatus_Exist = 1,
    // WL has been finished in some critical point
    // We should start new line
    WLFStatus_Broken = 2
  };

  const Standard_Integer aNbWLines = 2;
  for (Standard_Integer aCurInterval = 0; aCurInterval < theNbOfRanges; aCurInterval++)
  {
    //Process every continuous region
    Standard_Boolean isAddedIntoWL[aNbWLines];
    for (Standard_Integer i = 0; i < aNbWLines; i++)
      isAddedIntoWL[i] = Standard_False;

    Standard_Real anUf = 1.0, anUl = 0.0;
    if (!theRange[aCurInterval].GetBounds(anUf, anUl))
      continue;

    const Standard_Boolean isDeltaPeriod = IsEqual(anUl - anUf, aPeriod);

    //Inscribe and sort critical points
    for (Standard_Integer i = 0; i < aNbCritPoints; i++)
    {
      InscribePoint(anUf, anUl, anU1crit[i], 0.0, aPeriod, Standard_False);
    }

    std::sort(anU1crit, anU1crit + aNbCritPoints);

    while (anUf < anUl)
    {
      //Change value of U-parameter on the 1st surface from anUf to anUl
      //(anUf will be modified in the cycle body).
      //Step is computed adaptively (see comments below).

      Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
      WLFStatus aWLFindStatus[aNbWLines];
      Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
      Standard_Real anUexpect[aNbWLines];
      Standard_Boolean isAddingWLEnabled[aNbWLines];

      Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
      Handle(IntPatch_WLine) aWLine[aNbWLines];
      for (Standard_Integer i = 0; i < aNbWLines; i++)
      {
        aL2S[i] = new IntSurf_LineOn2S();
        aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
        aWLine[i]->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);
        aWLFindStatus[i] = WLFStatus_Absent;
        isAddingWLEnabled[i] = Standard_True;
        aU2[i] = aV1[i] = aV2[i] = 0.0;
        aV1Prev[i] = aV2Prev[i] = 0.0;
        anUexpect[i] = anUf;
      }

      Standard_Real aCriticalDelta[aNbCritPointsMax] = { 0 };
      for (Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++)
      { //We are not interested in elements of aCriticalDelta array
        //if their index is greater than or equal to aNbCritPoints

        aCriticalDelta[aCritPID] = anUf - anU1crit[aCritPID];
      }

      Standard_Real anU1 = anUf, aMinCriticalParam = anUf;
      Standard_Boolean isFirst = Standard_True;

      while (anU1 <= anUl)
      {
        //Change value of U-parameter on the 1st surface from anUf to anUl
        //(anUf will be modified in the cycle body). However, this cycle
        //can be broken if WL goes though some critical point.
        //Step is computed adaptively (see comments below).

        for (Standard_Integer i = 0; i < aNbCritPoints; i++)
        {
          if ((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
          {
            //WL has gone through i-th critical point
            anU1 = anU1crit[i];

            for (Standard_Integer j = 0; j < aNbWLines; j++)
            {
              aWLFindStatus[j] = WLFStatus_Broken;
              anUexpect[j] = anU1;
            }

            break;
          }
        }

        if (IsEqual(anU1, anUl))
        {
          for (Standard_Integer i = 0; i < aNbWLines; i++)
          {
            aWLFindStatus[i] = WLFStatus_Broken;
            anUexpect[i] = anU1;

            if (isDeltaPeriod)
            {
              //if isAddedIntoWL[i] == TRUE WLine contains only one point
              //(which was end point of previous WLine). If we will
              //add point found on the current step WLine will contain only
              //two points. At that both these points will be equal to the
              //points found earlier. Therefore, new WLine will repeat 
              //already existing WLine. Consequently, it is necessary 
              //to forbid building new line in this case.

              isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
            }
            else
            {
              isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
                                      (aWLFindStatus[i] == WLFStatus_Absent));
            }
          }//for(Standard_Integer i = 0; i < aNbWLines; i++)
        }
        else
        {
          for (Standard_Integer i = 0; i < aNbWLines; i++)
          {
            isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
                                    (aWLFindStatus[i] == WLFStatus_Absent));
          }//for(Standard_Integer i = 0; i < aNbWLines; i++)
        }

        for (Standard_Integer i = 0; i < aNbWLines; i++)
        {
          const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
            aWLine[i]->Curve()->NbPoints();

          if ((aWLFindStatus[i] == WLFStatus_Broken) ||
            (aWLFindStatus[i] == WLFStatus_Absent))
          {//Begin and end of WLine must be on boundary point
           //or on seam-edge strictly (if it is possible).

            Standard_Real aTol = aTol2D;
            ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs,
                                                        aU2[i], &aTol);
            InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);

            aTol = Max(aTol, aTol2D);

            if (Abs(aU2[i]) <= aTol)
              aU2[i] = 0.0;
            else if (Abs(aU2[i] - aPeriod) <= aTol)
              aU2[i] = aPeriod;
            else if (Abs(aU2[i] - aUSurf2f) <= aTol)
              aU2[i] = aUSurf2f;
            else if (Abs(aU2[i] - aUSurf2l) <= aTol)
              aU2[i] = aUSurf2l;
          }
          else
          {
            ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
            InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
          }

          if (aNbPntsWL == 0)
          {//the line has not contained any points yet
            if (((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*aTol2D) &&
                ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
                  (Abs(aU2[i] - aUSurf2l) < aTol2D)))
            {
              //In this case aU2[i] can have two values: current aU2[i] or
              //aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
              //correct value.

              Standard_Boolean isIncreasing = Standard_True;
              ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs,
                                                      aPeriod, isIncreasing);

              //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
              //then after the next step (when U1 will be increased) U2 will be
              //increased too. And we will go out of surface boundary.
              //Therefore, If U2(U1) is increasing then U2 must be equal aUSurf2f.
              //Analogically, if U2(U1) is decreasing.
              if (isIncreasing)
              {
                aU2[i] = aUSurf2f;
              }
              else
              {
                aU2[i] = aUSurf2l;
              }
            }
          }
          else
          {//aNbPntsWL > 0
            if (((aUSurf2l - aUSurf2f) >= aPeriod) &&
                ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
                  (Abs(aU2[i] - aUSurf2l) < aTol2D)))
            {//end of the line
              Standard_Real aU2prev = 0.0, aV2prev = 0.0;
              if (isReversed)
                aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
              else
                aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);

              if (2.0*Abs(aU2prev - aU2[i]) > aPeriod)
              {
                if (aU2prev > aU2[i])
                  aU2[i] += aPeriod;
                else
                  aU2[i] -= aPeriod;
              }
            }
          }

          ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
                                                              aV1[i], aV2[i]);

          if (isFirst)
          {
            aV1Prev[i] = aV1[i];
            aV2Prev[i] = aV2[i];
          }
        }//for(Standard_Integer i = 0; i < aNbWLines; i++)

        isFirst = Standard_False;

        //Looking for points into WLine
        Standard_Boolean isBroken = Standard_False;
        for (Standard_Integer i = 0; i < aNbWLines; i++)
        {
          if (!isAddingWLEnabled[i])
          {
            Standard_Boolean isBoundIntersect = Standard_False;
            if ((Abs(aV1[i] - aVSurf1f) <= aTol2D) ||
                ((aV1[i] - aVSurf1f)*(aV1Prev[i] - aVSurf1f) < 0.0))
            {
              isBoundIntersect = Standard_True;
            }
            else if ((Abs(aV1[i] - aVSurf1l) <= aTol2D) ||
                    ((aV1[i] - aVSurf1l)*(aV1Prev[i] - aVSurf1l) < 0.0))
            {
              isBoundIntersect = Standard_True;
            }
            else if ((Abs(aV2[i] - aVSurf2f) <= aTol2D) ||
                    ((aV2[i] - aVSurf2f)*(aV2Prev[i] - aVSurf2f) < 0.0))
            {
              isBoundIntersect = Standard_True;
            }
            else if ((Abs(aV2[i] - aVSurf2l) <= aTol2D) ||
                    ((aV2[i] - aVSurf2l)*(aV2Prev[i] - aVSurf2l) < 0.0))
            {
              isBoundIntersect = Standard_True;
            }

            if (aWLFindStatus[i] == WLFStatus_Broken)
              isBroken = Standard_True;

            if (!isBoundIntersect)
            {
              continue;
            }
            else
            {
              anUexpect[i] = anU1;
            }
          }

          // True if the current point already in the domain
          const Standard_Boolean isInscribe =
              ((aUSurf2f - aU2[i]) <= aTol2D) && ((aU2[i] - aUSurf2l) <= aTol2D) &&
              ((aVSurf1f - aV1[i]) <= aTol2D) && ((aV1[i] - aVSurf1l) <= aTol2D) &&
              ((aVSurf2f - aV2[i]) <= aTol2D) && ((aV2[i] - aVSurf2l) <= aTol2D);

          //isVIntersect == TRUE if intersection line intersects two (!)
          //V-bounds of cylinder (1st or 2nd - no matter)
          const Standard_Boolean isVIntersect =
              (((aVSurf1f - aV1[i])*(aVSurf1f - aV1Prev[i]) < RealSmall()) &&
                ((aVSurf1l - aV1[i])*(aVSurf1l - aV1Prev[i]) < RealSmall())) ||
              (((aVSurf2f - aV2[i])*(aVSurf2f - aV2Prev[i]) < RealSmall()) &&
                ((aVSurf2l - aV2[i])*(aVSurf2l - aV2Prev[i]) < RealSmall()));

          //isFound1 == TRUE if intersection line intersects V-bounds
          //  (First or Last - no matter) of the 1st cylynder
          //isFound2 == TRUE if intersection line intersects V-bounds
          //  (First or Last - no matter) of the 2nd cylynder
          Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
          Standard_Boolean isForce = Standard_False;

          if (aWLFindStatus[i] == WLFStatus_Absent)
          {
            if (((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1 - aUSurf1l) < aTol2D))
            {
              isForce = Standard_True;
            }
          }

          theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
                                 aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i, isForce,
                                 isFound1, isFound2);

          const Standard_Boolean isPrevVBound = !isVIntersect &&
                                                ((Abs(aV1Prev[i] - aVSurf1f) <= aTol2D) ||
                                                 (Abs(aV1Prev[i] - aVSurf1l) <= aTol2D) ||
                                                 (Abs(aV2Prev[i] - aVSurf2f) <= aTol2D) ||
                                                 (Abs(aV2Prev[i] - aVSurf2l) <= aTol2D));

          aV1Prev[i] = aV1[i];
          aV2Prev[i] = aV2[i];

          if ((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
          {
            aWLFindStatus[i] = WLFStatus_Broken; //start a new line
          }
          else if (isInscribe)
          {
            if ((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
            {
              aWLFindStatus[i] = WLFStatus_Exist;
            }

            if ((aWLFindStatus[i] != WLFStatus_Broken) ||
              (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
            {
              if (aWLine[i]->NbPnts() > 0)
              {
                Standard_Real aU2p = 0.0, aV2p = 0.0;
                if (isReversed)
                  aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
                else
                  aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);

                const Standard_Real aDelta = aU2[i] - aU2p;

                if (2.0 * Abs(aDelta) > aPeriod)
                {
                  if (aDelta > 0.0)
                  {
                    aU2[i] -= aPeriod;
                  }
                  else
                  {
                    aU2[i] += aPeriod;
                  }
                }
              }

              if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
                                gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
                                aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
                                aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
                                aWLine[i]->Curve(), i, aTol3D, aTol2D, isForce))
              {
                if (aWLFindStatus[i] == WLFStatus_Absent)
                {
                  aWLFindStatus[i] = WLFStatus_Exist;
                }
              }
              else if (!isFound1 && !isFound2)
              {//We do not add any point while doing this iteration
                if (aWLFindStatus[i] == WLFStatus_Exist)
                {
                  aWLFindStatus[i] = WLFStatus_Broken;
                }
              }
            }
          }
          else
          {//We do not add any point while doing this iteration
            if (aWLFindStatus[i] == WLFStatus_Exist)
            {
              aWLFindStatus[i] = WLFStatus_Broken;
            }
          }

          if (aWLFindStatus[i] == WLFStatus_Broken)
            isBroken = Standard_True;
        }//for(Standard_Integer i = 0; i < aNbWLines; i++)

        if (isBroken)
        {//current lines are filled. Go to the next lines
          anUf = anU1;

          Standard_Boolean isAdded = Standard_True;

          for (Standard_Integer i = 0; i < aNbWLines; i++)
          {
            if (isAddingWLEnabled[i])
            {
              continue;
            }

            isAdded = Standard_False;

            Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;

            theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
                                   aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i,
                                   Standard_False, isFound1, isFound2);

            if (isFound1 || isFound2)
            {
              isAdded = Standard_True;
            }

            if (aWLine[i]->NbPnts() > 0)
            {
              Standard_Real aU2p = 0.0, aV2p = 0.0;
              if (isReversed)
                aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
              else
                aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);

              const Standard_Real aDelta = aU2[i] - aU2p;

              if (2 * Abs(aDelta) > aPeriod)
              {
                if (aDelta > 0.0)
                {
                  aU2[i] -= aPeriod;
                }
                else
                {
                  aU2[i] += aPeriod;
                }
              }
            }

            if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
                              Standard_True, gp_Pnt2d(anU1, aV1[i]),
                              gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
                              aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
                              aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
                              i, aTol3D, aTol2D, Standard_False))
            {
              isAdded = Standard_True;
            }
          }

          if (!isAdded)
          {
            //Before breaking WL, we must complete it correctly
            //(e.g. to prolong to the surface boundary).
            //Therefore, we take the point last added in some WL
            //(have maximal U1-parameter) and try to add it in
            //the current WL.
            Standard_Real anUmaxAdded = RealFirst();

            {
              Standard_Boolean isChanged = Standard_False;
              for (Standard_Integer i = 0; i < aNbWLines; i++)
              {
                if ((aWLFindStatus[i] == WLFStatus_Absent) || (aWLine[i]->NbPnts() == 0))
                  continue;

                Standard_Real aU1c = 0.0, aV1c = 0.0;
                if (isReversed)
                  aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
                else
                  aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);

                anUmaxAdded = Max(anUmaxAdded, aU1c);
                isChanged = Standard_True;
              }

              if (!isChanged)
              { //If anUmaxAdded were not changed in previous cycle then
                //we would break existing WLines.
                break;
              }
            }

            for (Standard_Integer i = 0; i < aNbWLines; i++)
            {
              if (isAddingWLEnabled[i])
              {
                continue;
              }

              ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
                aU2[i], aV1[i], aV2[i]);

              AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
                             Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
                             gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
                             aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
                             aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
                             i, aTol3D, aTol2D, Standard_False);
            }
          }

          break;
        }

        //Step computing

        {
          //Step of aU1-parameter is computed adaptively. The algorithm 
          //aims to provide given aDeltaV1 and aDeltaV2 values (if it is 
          //possible because the intersection line can go along V-isoline)
          //in every iteration. It allows avoiding "flying" intersection
          //points too far each from other (see issue #24915).

          const Standard_Real aDeltaV1 = aRangeS1.Delta() / IntToReal(aNbPoints),
                              aDeltaV2 = aRangeS2.Delta() / IntToReal(aNbPoints);

          math_Matrix aMatr(1, 3, 1, 5);

          Standard_Real aMinUexp = RealLast();

          for (Standard_Integer i = 0; i < aNbWLines; i++)
          {
            if (aTol2D < (anUexpect[i] - anU1))
            {
              continue;
            }

            if (aWLFindStatus[i] == WLFStatus_Absent)
            {
              anUexpect[i] += aStepMax;
              aMinUexp = Min(aMinUexp, anUexpect[i]);
              continue;
            }
            //
            if (isGoodIntersection)
            {
              //Use constant step
              anUexpect[i] += aStepMax;
              aMinUexp = Min(aMinUexp, anUexpect[i]);

              continue;
            }
            //

            Standard_Real aStepTmp = aStepMax;

            const Standard_Real aSinU1 = sin(anU1),
                                aCosU1 = cos(anU1),
                                aSinU2 = sin(aU2[i]),
                                aCosU2 = cos(aU2[i]);

            aMatr.SetCol(1, anEquationCoeffs.mVecC1);
            aMatr.SetCol(2, anEquationCoeffs.mVecC2);
            aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
            aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
            aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
                            anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
                            anEquationCoeffs.mVecD);

            //The main idea is in solving of linearized system (2)
            //(see description to ComputationMethods class) in order to find new U1-value
            //to provide new value V1 or V2, which differs from current one by aDeltaV1 or
            //aDeltaV2 respectively. 

            //While linearizing, following Taylor formulas are used:
            //    cos(x0+dx) = cos(x0) - sin(x0)*dx
            //    sin(x0+dx) = sin(x0) + cos(x0)*dx

            //Consequently cos(U1), cos(U2), sin(U1) and sin(U2) in the system (2)
            //must be substituted by corresponding values.

            //ATTENTION!!!
            //The solution is approximate. More over, all requirements to the
            //linearization must be satisfied in order to obtain quality result.

            if (!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aStepTmp))
            {
              //To avoid cycling-up
              anUexpect[i] += aStepMax;
              aMinUexp = Min(aMinUexp, anUexpect[i]);

              continue;
            }

            if (aStepTmp < aStepMin)
              aStepTmp = aStepMin;

            if (aStepTmp > aStepMax)
              aStepTmp = aStepMax;

            anUexpect[i] = anU1 + aStepTmp;
            aMinUexp = Min(aMinUexp, anUexpect[i]);
          }

          anU1 = aMinUexp;
        }

        if (Precision::PConfusion() >= (anUl - anU1))
          anU1 = anUl;

        anUf = anU1;

        for (Standard_Integer i = 0; i < aNbWLines; i++)
        {
          if (aWLine[i]->NbPnts() != 1)
            isAddedIntoWL[i] = Standard_False;

          if (anU1 == anUl)
          {//strictly equal. Tolerance is considered above.
            anUexpect[i] = anUl;
          }
        }
      }

      for (Standard_Integer i = 0; i < aNbWLines; i++)
      {
        if ((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
        {
          isTheEmpty = Standard_False;
          Standard_Real u1, v1, u2, v2;
          aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
          IntPatch_Point aP;
          aP.SetParameter(u1);
          aP.SetParameters(u1, v1, u2, v2);
          aP.SetTolerance(aTol3D);
          aP.SetValue(aWLine[i]->Point(1).Value());

          //Check whether the added point exists.
          //It is enough to check the last point.
          if (theSPnt.IsEmpty() ||
              !theSPnt.Last().PntOn2S().IsSame(aP.PntOn2S(), Precision::Confusion()))
          {
            theSPnt.Append(aP);
          }
        }
        else if (aWLine[i]->NbPnts() > 1)
        {
          Standard_Boolean isGood = Standard_True;

          if (aWLine[i]->NbPnts() == 2)
          {
            const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
            const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);

            if (aPf.IsSame(aPl, Precision::Confusion()))
              isGood = Standard_False;
          }
          else if (aWLine[i]->NbPnts() > 2)
          {
            // Sometimes points of the WLine are distributed
            // linearly and uniformly. However, such position
            // of the points does not always describe the real intersection
            // curve. I.e. real tangents at the ends of the intersection
            // curve can significantly deviate from this "line" direction.
            // Here we are processing this case by inserting additional points
            // to the beginning/end of the WLine to make it more precise.
            // See description to the issue #30082.

            const Standard_Real aSqTol3D = aTol3D*aTol3D;
            for (Standard_Integer j = 0; j < 2; j++)
            {
              // If j == 0 ==> add point at begin of WLine.
              // If j == 1 ==> add point at end of WLine.

              for (;;)
              {
                if (aWLine[i]->NbPnts() >= aNbMaxPoints)
                {
                  break;
                }

                // Take 1st and 2nd point to compute the "line" direction.
                // For our convenience, we make 2nd point be the ends of the WLine
                // because it will be used for computation of the normals 
                // to the surfaces.
                const Standard_Integer anIdx1 = j ? aWLine[i]->NbPnts() - 1 : 2;
                const Standard_Integer anIdx2 = j ? aWLine[i]->NbPnts() : 1;

                const gp_Pnt &aP1 = aWLine[i]->Point(anIdx1).Value();
                const gp_Pnt &aP2 = aWLine[i]->Point(anIdx2).Value();

                const gp_Vec aDir(aP1, aP2);

                if (aDir.SquareMagnitude() < aSqTol3D)
                {
                  break;
                }

                // Compute tangent in first/last point of the WLine.
                // We do not take into account the flag "isReversed"
                // because strict direction of the tangent is not
                // important here (we are interested in the tangent
                // line itself and nothing to fear if its direction
                // is reversed).
                const gp_Vec aN1 = aQuad1.Normale(aP2);
                const gp_Vec aN2 = aQuad2.Normale(aP2);
                const gp_Vec aTg(aN1.Crossed(aN2));

                if (aTg.SquareMagnitude() < Precision::SquareConfusion())
                {
                  // Tangent zone
                  break;
                }

                // Check of the bending
                Standard_Real anAngle = aDir.Angle(aTg);

                if (anAngle > M_PI_2)
                  anAngle -= M_PI;

                if (Abs(anAngle) > 0.25) // ~ 14deg.
                {
                  const Standard_Integer aNbPntsPrev = aWLine[i]->NbPnts();
                  SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
                                       anEquationCoeffs, i, 3, anIdx1, anIdx2,
                                       aTol2D, aPeriod, isReversed);

                  if (aWLine[i]->NbPnts() == aNbPntsPrev)
                  {
                    // No points have been added. ==> Exit from a loop.
                    break;
                  }
                }
                else
                {
                  // Good result has been achieved. ==> Exit from a loop.
                  break;
                }
              } // for (;;)
            }
          }

          if (isGood)
          {
            isTheEmpty = Standard_False;
            isAddedIntoWL[i] = Standard_True;
            SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
                                 anEquationCoeffs, i, aNbPoints, 1,
                                 aWLine[i]->NbPnts(), aTol2D, aPeriod,
                                 isReversed);

            aWLine[i]->ComputeVertexParameters(aTol3D);
            theSlin.Append(aWLine[i]);
          }
        }
        else
        {
          isAddedIntoWL[i] = Standard_False;
        }

#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
        aWLine[i]->Dump(0);
#endif
      }
    }
  }


  //Delete the points in theSPnt, which
  //lie at least in one of the line in theSlin.
  for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
  {
    for (Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
    {
      Handle(IntPatch_WLine) aWLine1(Handle(IntPatch_WLine)::
        DownCast(theSlin.Value(aNbLin)));

      const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
      const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());

      const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
      if (aPntCur.IsSame(aPntFWL1, aTol3D) ||
        aPntCur.IsSame(aPntLWL1, aTol3D))
      {
        theSPnt.Remove(aNbPnt);
        aNbPnt--;
        break;
      }
    }
  }

  //Try to add new points in the neighborhood of existing point
  for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
  {
    // Standard algorithm (implemented above) could not find any
    // continuous curve in neighborhood of aPnt2S (e.g. because
    // this curve is too small; see tests\bugs\modalg_5\bug25292_35 and _36).
    // Here, we will try to find several new points nearer to aPnt2S.

    // The algorithm below tries to find two points in every
    // intervals [u1 - aStepMax, u1] and [u1, u1 + aStepMax]
    // and every new point will be in maximal distance from
    // u1. If these two points exist they will be joined
    // by the intersection curve.

    const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);

    Standard_Real u1, v1, u2, v2;
    aPnt2S.Parameters(u1, v1, u2, v2);

    Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
    Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
    aWLine->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);

    //Define the index of WLine, which lies the point aPnt2S in.
    Standard_Integer anIndex = 0;

    Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
    if (isReversed)
    {
      anUf = Max(u2 - aStepMax, aUSurf1f);
      anUl = Min(u2 + aStepMax, aUSurf1l);
      aCurU2 = u1;
    }
    else
    {
      anUf = Max(u1 - aStepMax, aUSurf1f);
      anUl = Min(u1 + aStepMax, aUSurf1l);
      aCurU2 = u2;
    }

    const Standard_Real anUinf = anUf, anUsup = anUl, anUmid = 0.5*(anUf + anUl);

    {
      //Find the value of anIndex variable.
      Standard_Real aDelta = RealLast();
      for (Standard_Integer i = 0; i < aNbWLines; i++)
      {
        Standard_Real anU2t = 0.0;
        if (!ComputationMethods::CylCylComputeParameters(anUmid, i, anEquationCoeffs, anU2t))
          continue;

        Standard_Real aDU2 = fmod(Abs(anU2t - aCurU2), aPeriod);
        aDU2 = Min(aDU2, Abs(aDU2 - aPeriod));
        if (aDU2 < aDelta)
        {
          aDelta = aDU2;
          anIndex = i;
        }
      }
    }

    // Bisection method is used in order to find every new point.
    // I.e. if we need to find intersection point in the interval [anUinf, anUmid]
    // we check the point anUC = 0.5*(anUinf+anUmid). If it is an intersection point
    // we try to find another point in the interval [anUinf, anUC] (because we find the point in
    // maximal distance from anUmid). If it is not then we try to find another point in the
    // interval [anUC, anUmid]. Next iterations will be made analogically.
    // When we find intersection point in the interval [anUmid, anUsup] we try to find
    // another point in the interval [anUC, anUsup] if anUC is intersection point and
    // in the interval [anUmid, anUC], otherwise.

    Standard_Real anAddedPar[2] = {isReversed ? u2 : u1, isReversed ? u2 : u1};

    for (Standard_Integer aParID = 0; aParID < 2; aParID++)
    {
      if (aParID == 0)
      {
        anUf = anUinf;
        anUl = anUmid;
      }
      else // if(aParID == 1)
      {
        anUf = anUmid;
        anUl = anUsup;
      }

      Standard_Real &aPar1 = (aParID == 0) ? anUf : anUl,
                    &aPar2 = (aParID == 0) ? anUl : anUf;

      while (Abs(aPar2 - aPar1) > aStepMin)
      {
        Standard_Real anUC = 0.5*(anUf + anUl);
        Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
        Standard_Boolean isDone = ComputationMethods::
                CylCylComputeParameters(anUC, anIndex, anEquationCoeffs, aU2, aV1, aV2);

        if (isDone)
        {
          if (Abs(aV1 - aVSurf1f) <= aTol2D)
            aV1 = aVSurf1f;

          if (Abs(aV1 - aVSurf1l) <= aTol2D)
            aV1 = aVSurf1l;

          if (Abs(aV2 - aVSurf2f) <= aTol2D)
            aV2 = aVSurf2f;

          if (Abs(aV2 - aVSurf2l) <= aTol2D)
            aV2 = aVSurf2l;

          isDone = AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
                                  Standard_True, gp_Pnt2d(anUC, aV1), gp_Pnt2d(aU2, aV2),
                                  aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
                                  aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
                                  aPeriod, aWLine->Curve(), anIndex, aTol3D,
                                  aTol2D, Standard_False, Standard_True);
        }

        if (isDone)
        {
          anAddedPar[0] = Min(anAddedPar[0], anUC);
          anAddedPar[1] = Max(anAddedPar[1], anUC);
          aPar2 = anUC;
        }
        else
        {
          aPar1 = anUC;
        }
      }
    }

    //Fill aWLine by additional points
    if (anAddedPar[1] - anAddedPar[0] > aStepMin)
    {
      for (Standard_Integer aParID = 0; aParID < 2; aParID++)
      {
        Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
        ComputationMethods::CylCylComputeParameters(anAddedPar[aParID], anIndex,
                                                  anEquationCoeffs, aU2, aV1, aV2);

        AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
                        gp_Pnt2d(anAddedPar[aParID], aV1), gp_Pnt2d(aU2, aV2),
                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod, aWLine->Curve(),
                        anIndex, aTol3D, aTol2D, Standard_False, Standard_False);
      }

      SeekAdditionalPoints(aQuad1, aQuad2, aWLine->Curve(),
                            anEquationCoeffs, anIndex, aNbMinPoints,
                            1, aWLine->NbPnts(), aTol2D, aPeriod,
                            isReversed);

      aWLine->ComputeVertexParameters(aTol3D);
      theSlin.Append(aWLine);

      theSPnt.Remove(aNbPnt);
      aNbPnt--;
    }
  }

  return IntPatch_ImpImpIntersection::IntStatus_OK;
}

//=======================================================================
//function : IntCyCy
//purpose  : 
//=======================================================================
IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1,
                                               const IntSurf_Quadric& theQuad2,
                                               const Standard_Real theTol3D,
                                               const Standard_Real theTol2D,
                                               const Bnd_Box2d& theUVSurf1,
                                               const Bnd_Box2d& theUVSurf2,
                                               Standard_Boolean& isTheEmpty,
                                               Standard_Boolean& isTheSameSurface,
                                               Standard_Boolean& isTheMultiplePoint,
                                               IntPatch_SequenceOfLine& theSlin,
                                               IntPatch_SequenceOfPoint& theSPnt)
{
  isTheEmpty = Standard_True;
  isTheSameSurface = Standard_False;
  isTheMultiplePoint = Standard_False;
  theSlin.Clear();
  theSPnt.Clear();

  const gp_Cylinder aCyl1 = theQuad1.Cylinder(),
                    aCyl2 = theQuad2.Cylinder();

  IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);

  if (!anInter.IsDone())
  {
    return IntPatch_ImpImpIntersection::IntStatus_Fail;
  }

  if(anInter.TypeInter() != IntAna_NoGeometricSolution)
  {
    if (CyCyAnalyticalIntersect(theQuad1, theQuad2, anInter,
                                theTol3D, isTheEmpty,
                                isTheSameSurface, isTheMultiplePoint,
                                theSlin, theSPnt))
    {
      return IntPatch_ImpImpIntersection::IntStatus_OK;
    }
  }
  
  //Here, intersection line is not an analytical curve(line, circle, ellipsis etc.)

  Standard_Real aUSBou[2][2], aVSBou[2][2]; //const

  theUVSurf1.Get(aUSBou[0][0], aVSBou[0][0], aUSBou[0][1], aVSBou[0][1]);
  theUVSurf2.Get(aUSBou[1][0], aVSBou[1][0], aUSBou[1][1], aVSBou[1][1]);

  const Standard_Real aPeriod = 2.0*M_PI;
  const Standard_Integer aNbWLines = 2;

  const ComputationMethods::stCoeffsValue anEquationCoeffs1(aCyl1, aCyl2);
  const ComputationMethods::stCoeffsValue anEquationCoeffs2(aCyl2, aCyl1);

  //Boundaries. 
  //Intersection result can include two non-connected regions
  //(see WorkWithBoundaries::BoundariesComputing(...) method).
  const Standard_Integer aNbOfBoundaries = 2;
  Bnd_Range anURange[2][aNbOfBoundaries];   //const

  if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs1, aPeriod, anURange[0]))
    return IntPatch_ImpImpIntersection::IntStatus_OK;

  if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs2, aPeriod, anURange[1]))
    return IntPatch_ImpImpIntersection::IntStatus_OK;

  //anURange[*] can be in different periodic regions in
  //compare with First-Last surface. E.g. the surface
  //is full cylinder [0, 2*PI] but anURange is [5, 7].
  //Trivial common range computation returns [5, 2*PI] and
  //its summary length is 2*PI-5 == 1.28... only. That is wrong.
  //This problem can be solved by the following
  //algorithm:
  // 1. split anURange[*] by the surface boundary;
  // 2. shift every new range in order to inscribe it
  //      in [Ufirst, Ulast] of cylinder;
  // 3. consider only common ranges between [Ufirst, Ulast]
  //    and new ranges.
  //
  // In above example, we obtain following:
  // 1. two ranges: [5, 2*PI] and [2*PI, 7];
  // 2. after shifting: [5, 2*PI] and [0, 7-2*PI];
  // 3. Common ranges: ([5, 2*PI] and [0, 2*PI]) == [5, 2*PI],
  //                   ([0, 7-2*PI] and [0, 2*PI]) == [0, 7-2*PI];
  // 4. Their summary length is (2*PI-5)+(7-2*PI-0)==7-5==2 ==> GOOD.

  Standard_Real aSumRange[2] = { 0.0, 0.0 };
  Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
  for (Standard_Integer aCID = 0; aCID < 2; aCID++)
  {
    anAlloc->Reset(false);
    NCollection_List<Bnd_Range> aListOfRng(anAlloc);
    
    aListOfRng.Append(anURange[aCID][0]);
    aListOfRng.Append(anURange[aCID][1]);

    const Standard_Real aSplitArr[3] = {aUSBou[aCID][0], aUSBou[aCID][1], 0.0};

    NCollection_List<Bnd_Range>::Iterator anITrRng;
    for (Standard_Integer aSInd = 0; aSInd < 3; aSInd++)
    {
      NCollection_List<Bnd_Range> aLstTemp(aListOfRng);
      aListOfRng.Clear();
      for (anITrRng.Init(aLstTemp); anITrRng.More(); anITrRng.Next())
      {
        Bnd_Range& aRng = anITrRng.ChangeValue();
        aRng.Split(aSplitArr[aSInd], aListOfRng, aPeriod);
      }
    }

    anITrRng.Init(aListOfRng);
    for (; anITrRng.More(); anITrRng.Next())
    {
      Bnd_Range& aCurrRange = anITrRng.ChangeValue();

      Bnd_Range aBoundR;
      aBoundR.Add(aUSBou[aCID][0]);
      aBoundR.Add(aUSBou[aCID][1]);

      if (!InscribeInterval(aUSBou[aCID][0], aUSBou[aCID][1],
                                          aCurrRange, theTol2D, aPeriod))
      {
        //If aCurrRange does not have common block with
        //[Ufirst, Ulast] of cylinder then we will try
        //to inscribe [Ufirst, Ulast] in the boundaries of aCurrRange.
        Standard_Real aF = 1.0, aL = 0.0;
        if (!aCurrRange.GetBounds(aF, aL))
          continue;

        if ((aL < aUSBou[aCID][0]))
        {
          aCurrRange.Shift(aPeriod);
        }
        else if (aF > aUSBou[aCID][1])
        {
          aCurrRange.Shift(-aPeriod);
        }
      }

      aBoundR.Common(aCurrRange);

      const Standard_Real aDelta = aBoundR.Delta();

      if (aDelta > 0.0)
      {
        aSumRange[aCID] += aDelta;
      }
    }
  }

  //The bigger range the bigger number of points in Walking-line (WLine)
  //we will be able to add and consequently we will obtain more
  //precise intersection line.
  //Every point of WLine is determined as function from U1-parameter,
  //where U1 is U-parameter on 1st quadric.
  //Therefore, we should use quadric with bigger range as 1st parameter
  //in IntCyCy() function.
  //On the other hand, there is no point in reversing in case of
  //analytical intersection (when result is line, ellipse, point...).
  //This result is independent of the arguments order.
  const Standard_Boolean isToReverse = (aSumRange[1] > aSumRange[0]);

  if (isToReverse)
  {
    const WorkWithBoundaries aBoundWork(theQuad2, theQuad1, anEquationCoeffs2,
                                        theUVSurf2, theUVSurf1, aNbWLines,
                                        aPeriod, theTol3D, theTol2D, Standard_True);

    return CyCyNoGeometric(aCyl2, aCyl1, aBoundWork, anURange[1], aNbOfBoundaries,
                              isTheEmpty, theSlin, theSPnt);
  }
  else
  {
    const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs1,
                                        theUVSurf1, theUVSurf2, aNbWLines,
                                        aPeriod, theTol3D, theTol2D, Standard_False);

    return CyCyNoGeometric(aCyl1, aCyl2, aBoundWork, anURange[0], aNbOfBoundaries,
                              isTheEmpty, theSlin, theSPnt);
  }
}

//=======================================================================
//function : IntCySp
//purpose  : 
//=======================================================================
Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
			 const IntSurf_Quadric& Quad2,
			 const Standard_Real Tol,
			 const Standard_Boolean Reversed,
			 Standard_Boolean& Empty,
			 Standard_Boolean& Multpoint,
			 IntPatch_SequenceOfLine& slin,
			 IntPatch_SequenceOfPoint& spnt)

{
  Standard_Integer i;

  IntSurf_TypeTrans trans1,trans2;
  IntAna_ResultType typint;
  IntPatch_Point ptsol;
  gp_Circ cirsol;

  gp_Cylinder Cy;
  gp_Sphere Sp;

  if (!Reversed) {
    Cy = Quad1.Cylinder();
    Sp = Quad2.Sphere();
  }
  else {
    Cy = Quad2.Cylinder();
    Sp = Quad1.Sphere();
  }
  IntAna_QuadQuadGeo inter(Cy,Sp,Tol);

  if (!inter.IsDone()) {return Standard_False;}

  typint = inter.TypeInter();
  Standard_Integer NbSol = inter.NbSolutions();
  Empty = Standard_False;

  switch (typint) {

  case IntAna_Empty :
    {
      Empty = Standard_True;
    }
    break;

  case IntAna_Point :
    {
      gp_Pnt psol(inter.Point(1));
      Standard_Real U1,V1,U2,V2;
      Quad1.Parameters(psol,U1,V1);
      Quad2.Parameters(psol,U2,V2);
      ptsol.SetValue(psol,Tol,Standard_True);
      ptsol.SetParameters(U1,V1,U2,V2);
      spnt.Append(ptsol);
    }
    break;

  case IntAna_Circle:
    {
      cirsol = inter.Circle(1);
      gp_Vec Tgt;
      gp_Pnt ptref;
      ElCLib::D1(0.,cirsol,ptref,Tgt);

      if (NbSol == 1) {
	gp_Vec TestCurvature(ptref,Sp.Location());
	gp_Vec Normsp,Normcyl;
	if (!Reversed) {
	  Normcyl = Quad1.Normale(ptref);
	  Normsp  = Quad2.Normale(ptref);
	}
	else {
	  Normcyl = Quad2.Normale(ptref);
	  Normsp  = Quad1.Normale(ptref);
	}
	
	IntSurf_Situation situcyl;
	IntSurf_Situation situsp;
	
	if (Normcyl.Dot(TestCurvature) > 0.) {
	  situsp = IntSurf_Outside;
	  if (Normsp.Dot(Normcyl) > 0.) {
	    situcyl = IntSurf_Inside;
	  }
	  else {
	    situcyl = IntSurf_Outside;
	  }
	}
	else {
	  situsp = IntSurf_Inside;
	  if (Normsp.Dot(Normcyl) > 0.) {
	    situcyl = IntSurf_Outside;
	  }
	  else {
	    situcyl = IntSurf_Inside;
	  }
	}
	Handle(IntPatch_GLine) glig;
	if (!Reversed) {
	  glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
	}
	else {
	  glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
	}
	slin.Append(glig);
      }
      else {
	if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
	  trans1 = IntSurf_Out;
	  trans2 = IntSurf_In;
	}
	else {
	  trans1 = IntSurf_In;
	  trans2 = IntSurf_Out;
	}
	Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
	slin.Append(glig);

	cirsol = inter.Circle(2);
	ElCLib::D1(0.,cirsol,ptref,Tgt);
	Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
	if(qwe> 0.0000001) {
	  trans1 = IntSurf_Out;
	  trans2 = IntSurf_In;
	}
	else if(qwe<-0.0000001) {
	  trans1 = IntSurf_In;
	  trans2 = IntSurf_Out;
	}
	else { 
	  trans1=trans2=IntSurf_Undecided; 
	}
	glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
	slin.Append(glig);
      }
    }
    break;
    
  case IntAna_NoGeometricSolution:
    {
      gp_Pnt psol;
      Standard_Real U1,V1,U2,V2;
      IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
      if (!anaint.IsDone()) {
	return Standard_False;
      }

      if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
	Empty = Standard_True;
      }
      else {

	NbSol = anaint.NbPnt();
	for (i = 1; i <= NbSol; i++) {
	  psol = anaint.Point(i);
	  Quad1.Parameters(psol,U1,V1);
	  Quad2.Parameters(psol,U2,V2);
	  ptsol.SetValue(psol,Tol,Standard_True);
	  ptsol.SetParameters(U1,V1,U2,V2);
	  spnt.Append(ptsol);
	}
	
	gp_Pnt ptvalid,ptf,ptl;
	gp_Vec tgvalid;
	Standard_Real first,last,para;
	IntAna_Curve curvsol;
	Standard_Boolean tgfound;
	Standard_Integer kount;
	
	NbSol = anaint.NbCurve();
	for (i = 1; i <= NbSol; i++) {
	  curvsol = anaint.Curve(i);
	  curvsol.Domain(first,last);
	  ptf = curvsol.Value(first);
	  ptl = curvsol.Value(last);

	  para = last;
	  kount = 1;
	  tgfound = Standard_False;
	  
	  while (!tgfound) {
	    para = (1.123*first + para)/2.123;
	    tgfound = curvsol.D1u(para,ptvalid,tgvalid);
	    if (!tgfound) {
	      kount ++;
	      tgfound = kount > 5;
	    }
	  }
	  Handle(IntPatch_ALine) alig;
	  if (kount <= 5) {
	    Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
						 Quad1.Normale(ptvalid));
	    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; 
	    }
	    alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
	  }
	  else {
	    alig = new IntPatch_ALine(curvsol,Standard_False);
	  }
	  Standard_Boolean TempFalse1a = Standard_False;
	  Standard_Boolean TempFalse2a = Standard_False;

	  //-- ptf et ptl : points debut et fin de alig 
	  
	  ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
			TempFalse2a,ptl,last,Multpoint,Tol);
	  slin.Append(alig);
	} //-- boucle sur les lignes 
      } //-- solution avec au moins une lihne 
    }
    break;

  default:
    {
      return Standard_False;
    }
  }
  return Standard_True;
}
//=======================================================================
//function : IntCyCo
//purpose  : 
//=======================================================================
Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
			 const IntSurf_Quadric& Quad2,
			 const Standard_Real Tol,
			 const Standard_Boolean Reversed,
			 Standard_Boolean& Empty,
			 Standard_Boolean& Multpoint,
			 IntPatch_SequenceOfLine& slin,
			 IntPatch_SequenceOfPoint& spnt)

{
  IntPatch_Point ptsol;

  Standard_Integer i;

  IntSurf_TypeTrans trans1,trans2;
  IntAna_ResultType typint;
  gp_Circ cirsol;

  gp_Cylinder Cy;
  gp_Cone     Co;

  if (!Reversed) {
    Cy = Quad1.Cylinder();
    Co = Quad2.Cone();
  }
  else {
    Cy = Quad2.Cylinder();
    Co = Quad1.Cone();
  }
  IntAna_QuadQuadGeo inter(Cy,Co,Tol);

  if (!inter.IsDone()) {return Standard_False;}

  typint = inter.TypeInter();
  Standard_Integer NbSol = inter.NbSolutions();
  Empty = Standard_False;

  switch (typint) {

  case IntAna_Empty : {
    Empty = Standard_True;
  }
    break;

  case IntAna_Point :{
    gp_Pnt psol(inter.Point(1));
    Standard_Real U1,V1,U2,V2;
    Quad1.Parameters(psol,U1,V1);
    Quad1.Parameters(psol,U2,V2);
    ptsol.SetValue(psol,Tol,Standard_True);
    ptsol.SetParameters(U1,V1,U2,V2);
    spnt.Append(ptsol);
  }
    break;
    
  case IntAna_Circle:  {
    gp_Vec Tgt;
    gp_Pnt ptref;
    Standard_Integer j;
    Standard_Real qwe;
    //
    for(j=1; j<=2; ++j) { 
      cirsol = inter.Circle(j);
      ElCLib::D1(0.,cirsol,ptref,Tgt);
      qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.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(cirsol,Standard_False,trans1,trans2);
      slin.Append(glig);
    }
  }
    break;
    
  case IntAna_NoGeometricSolution:    {
    gp_Pnt psol;
    Standard_Real U1,V1,U2,V2;
    IntAna_IntQuadQuad anaint(Cy,Co,Tol);
    if (!anaint.IsDone()) {
      return Standard_False;
    }
    
    if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
      Empty = Standard_True;
    }
    else {
      NbSol = anaint.NbPnt();
      for (i = 1; i <= NbSol; i++) {
	psol = anaint.Point(i);
	Quad1.Parameters(psol,U1,V1);
	Quad2.Parameters(psol,U2,V2);
	ptsol.SetValue(psol,Tol,Standard_True);
	ptsol.SetParameters(U1,V1,U2,V2);
	spnt.Append(ptsol);
      }
      
      gp_Pnt ptvalid, ptf, ptl;
      gp_Vec tgvalid;
      //
      Standard_Real first,last,para;
      Standard_Boolean tgfound,firstp,lastp,kept;
      Standard_Integer kount;
      //
      //
      //IntAna_Curve curvsol;
      IntAna_Curve aC;
      IntAna_ListOfCurve aLC;
      IntAna_ListIteratorOfListOfCurve aIt;
      
      //
      NbSol = anaint.NbCurve();
      for (i = 1; i <= NbSol; ++i) {
	kept = Standard_False;
	//curvsol = anaint.Curve(i);
	aC=anaint.Curve(i);
	aLC.Clear();
	ExploreCurve(Co, aC, 10.*Tol, aLC);
	//
	aIt.Initialize(aLC);
	for (; aIt.More(); aIt.Next()) {
	  IntAna_Curve& curvsol=aIt.ChangeValue();
	  //
	  curvsol.Domain(first, last);
	  firstp = !curvsol.IsFirstOpen();
	  lastp  = !curvsol.IsLastOpen();
	  if (firstp) {
	    ptf = curvsol.Value(first);
	  }
	  if (lastp) {
	    ptl = curvsol.Value(last);
	  }
	  para = last;
	  kount = 1;
	  tgfound = Standard_False;
	  
	  while (!tgfound) {
	    para = (1.123*first + para)/2.123;
	    tgfound = curvsol.D1u(para,ptvalid,tgvalid);
	    if (!tgfound) {
	      kount ++;
	      tgfound = kount > 5;
	    }
	  }
	  Handle(IntPatch_ALine) alig;
	  if (kount <= 5) {
	    Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
						 Quad1.Normale(ptvalid));
	    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; 
	    }
	    alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
	    kept = Standard_True;
	  }
	  else {
	    ptvalid = curvsol.Value(para);
	    alig = new IntPatch_ALine(curvsol,Standard_False);
	    kept = Standard_True;
	    //-- std::cout << "Transition indeterminee" << std::endl;
	  }
	  if (kept) {
	    Standard_Boolean Nfirstp = !firstp;
	    Standard_Boolean Nlastp  = !lastp;
	    ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
			  Nlastp,ptl,last,Multpoint,Tol);
	    slin.Append(alig);
	  }
	} // for (; aIt.More(); aIt.Next())
      } // for (i = 1; i <= NbSol; ++i) 
    }
  }
    break;
    
  default:
    return Standard_False;
    
  } // switch (typint)
  
  return Standard_True;
}
//=======================================================================
//function : ExploreCurve
//purpose  : Splits aC on several curves in the cone apex points.
//=======================================================================
Standard_Boolean ExploreCurve(const gp_Cone& theCo,
                              IntAna_Curve& theCrv,
                              const Standard_Real theTol,
                              IntAna_ListOfCurve& theLC)
{
  const Standard_Real aSqTol = theTol*theTol;
  const gp_Pnt aPapx(theCo.Apex());

  Standard_Real aT1, aT2;
  theCrv.Domain(aT1, aT2);

  theLC.Clear();
  //
  TColStd_ListOfReal aLParams;
  theCrv.FindParameter(aPapx, aLParams);
  if (aLParams.IsEmpty())
  {
    theLC.Append(theCrv);
    return Standard_False;
  }

  for (TColStd_ListIteratorOfListOfReal anItr(aLParams); anItr.More(); anItr.Next())
  {
    Standard_Real aPrm = anItr.Value();

    if ((aPrm - aT1) < Precision::PConfusion())
      continue;

    Standard_Boolean isLast = Standard_False;
    if ((aT2 - aPrm) < Precision::PConfusion())
    {
      aPrm = aT2;
      isLast = Standard_True;
    }

    const gp_Pnt aP = theCrv.Value(aPrm);
    const Standard_Real aSqD = aP.SquareDistance(aPapx);
    if (aSqD < aSqTol)
    {
      IntAna_Curve aC1 = theCrv;
      aC1.SetDomain(aT1, aPrm);
      aT1 = aPrm;
      theLC.Append(aC1);
    }

    if (isLast)
      break;
  }

  if (theLC.IsEmpty())
  {
    theLC.Append(theCrv);
    return Standard_False;
  }

  if ((aT2 - aT1) > Precision::PConfusion())
  {
    IntAna_Curve aC1 = theCrv;
    aC1.SetDomain(aT1, aT2);
    theLC.Append(aC1);
  }

  return Standard_True;
}