src/IntPatch/IntPatch_ImpImpIntersection_4.gxx

   1 // Created on: 1992-05-07
   2 // Created by: Jacques GOUSSARD
   3 // Copyright (c) 1992-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <algorithm>
  18 #include <Bnd_Range.hxx>
  19 #include <IntAna_ListOfCurve.hxx>
  20 #include <math_Matrix.hxx>
  21 #include <NCollection_IncAllocator.hxx>
  22 #include <Standard_DivideByZero.hxx>
  23
  24 //If Abs(a) <= aNulValue then it is considered that a = 0.
  25 static const Standard_Real aNulValue = 1.0e-11;
  26
  27 static void ShortCosForm( const Standard_Real theCosFactor,
  28                           const Standard_Real theSinFactor,
  29                           Standard_Real& theCoeff,
  30                           Standard_Real& theAngle);
  31 //
  32 static
  33   Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
  34                                 const gp_Cone& aCo,
  35                                 IntAna_Curve& aC,
  36                                 const Standard_Real aTol,
  37                                 IntAna_ListOfCurve& aLC);
  38
  39 static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
  40                                       const Standard_Real theUlTarget,
  41                                       Standard_Real& theUGiven,
  42                                       const Standard_Real theTol2D,
  43                                       const Standard_Real thePeriod,
  44                                       const Standard_Boolean theFlForce);
  45
  46
  47 class ComputationMethods
  48 {
  49   //Every cylinder can be represented by the following equation in parametric form:
  50   //    S(U,V) = L + R*cos(U)*Xd+R*sin(U)*Yd+V*Zd,
  51   //where location L, directions Xd, Yd and Zd have type gp_XYZ.
  52
  53   //Intersection points between two cylinders can be found from the following system:
  54   //    S1(U1, V1) = S2(U2, V2)
  55   //or
  56   //    {X01 + R1*cos(U1)*Xx1 + R1*sin(U1)*Yx1 + V1*Zx1 = X02 + R2*cos(U2)*Xx2 + R2*sin(U2)*Yx2 + V2*Zx2
  57   //    {Y01 + R1*cos(U1)*Xy1 + R1*sin(U1)*Yy1 + V1*Zy1 = Y02 + R2*cos(U2)*Xy2 + R2*sin(U2)*Yy2 + V2*Zy2   (1)
  58   //    {Z01 + R1*cos(U1)*Xz1 + R1*sin(U1)*Yz1 + V1*Zz1 = Z02 + R2*cos(U2)*Xz2 + R2*sin(U2)*Yz2 + V2*Zz2
  59
  60   //The formula (1) can be rewritten as follows
  61   //    {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1
  62   //    {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2                                   (2)
  63   //    {C13*V1+C23*V2=A13*cos(U1)+B13*sin(U1)+A23*cos(U2)+B23*sin(U2)+D3
  64
  65   //Hereafter we consider that in system
  66   //    {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1                                   (3)
  67   //    {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2
  68   //variables V1 and V2 can be found unambiguously, i.e. determinant
  69   //            |C11 C21|
  70   //            |       | != 0
  71   //            |C12 C22|
  72   //
  73   //In this case, variables V1 and V2 can be found as follows:
  74   //    {V1 = K11*sin(U1)+K21*sin(U2)+L11*cos(U1)+L21*cos(U2)+M1 = K1*cos(U1-FIV1)+L1*cos(U2-PSIV1)+M1      (4)
  75   //    {V2 = K12*sin(U1)+K22*sin(U2)+L12*cos(U1)+L22*cos(U2)+M2 = K2*cos(U2-FIV2)+L2*cos(U2-PSIV2)+M2
  76
  77   //Having substituted result of (4) to the 3rd equation of (2), we will obtain equation
  78   //    cos(U2-FI2) = B*cos(U1-FI1)+C.                                                                      (5)
  79
  80   //I.e. when U1 is taken different given values (from domain), correspond U2 value can be computed
  81   //from equation (5). After that, V1 and V2 can be computed from the system (4) (see
  82   //CylCylComputeParameters(...) methods).
  83
  84   //It is important to remark that equation (5) (in general) has two solutions: U2=FI2 +/- f(U1).
  85   //Therefore, we are getting here two intersection lines.
  86
  87 public:
  88   //Stores equations coefficients
  89   struct stCoeffsValue
  90   {
  91     stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
  92
  93     math_Vector mVecA1;
  94     math_Vector mVecA2;
  95     math_Vector mVecB1;
  96     math_Vector mVecB2;
  97     math_Vector mVecC1;
  98     math_Vector mVecC2;
  99     math_Vector mVecD;
 100
 101     Standard_Real mK21; //sinU2
 102     Standard_Real mK11; //sinU1
 103     Standard_Real mL21; //cosU2
 104     Standard_Real mL11;  //cosU1
 105     Standard_Real mM1;  //Free member
 106
 107     Standard_Real mK22; //sinU2
 108     Standard_Real mK12; //sinU1
 109     Standard_Real mL22; //cosU2
 110     Standard_Real mL12; //cosU1
 111     Standard_Real mM2; //Free member
 112
 113     Standard_Real mK1;
 114     Standard_Real mL1;
 115     Standard_Real mK2;
 116     Standard_Real mL2;
 117
 118     Standard_Real mFIV1;
 119     Standard_Real mPSIV1;
 120     Standard_Real mFIV2;
 121     Standard_Real mPSIV2;
 122
 123     Standard_Real mB;
 124     Standard_Real mC;
 125     Standard_Real mFI1;
 126     Standard_Real mFI2;
 127   };
 128
 129
 130   //! Determines, if U2(U1) function is increasing.
 131   static Standard_Boolean CylCylMonotonicity(const Standard_Real theU1par,
 132                                              const Standard_Integer theWLIndex,
 133                                              const stCoeffsValue& theCoeffs,
 134                                              const Standard_Real thePeriod,
 135                                              Standard_Boolean& theIsIncreasing);
 136
 137   //! Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
 138   //! esimates the tolerance of U2-computing (estimation result is
 139   //! assigned to *theDelta value).
 140   static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
 141                                                 const Standard_Integer theWLIndex,
 142                                                 const stCoeffsValue& theCoeffs,
 143                                                 Standard_Real& theU2,
 144                                                 Standard_Real* const theDelta = 0);
 145
 146   static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
 147                                                   const Standard_Real theU2,
 148                                                   const stCoeffsValue& theCoeffs,
 149                                                   Standard_Real& theV1,
 150                                                   Standard_Real& theV2);
 151
 152   static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
 153                                                   const Standard_Integer theWLIndex,
 154                                                   const stCoeffsValue& theCoeffs,
 155                                                   Standard_Real& theU2,
 156                                                   Standard_Real& theV1,
 157                                                   Standard_Real& theV2);
 158
 159 };
 160
 161 ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1,
 162                                                  const gp_Cylinder& theCyl2):
 163   mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()),
 164   mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()),
 165   mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()),
 166   mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()),
 167   mVecC1(theCyl1.Axis().Direction().XYZ()),
 168   mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()),
 169   mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ())
 170 {
 171   enum CoupleOfEquation
 172   {
 173     COENONE = 0,
 174     COE12 = 1,
 175     COE23 = 2,
 176     COE13 = 3
 177   }aFoundCouple = COENONE;
 178
 179
 180   Standard_Real aDetV1V2 = 0.0;
 181
 182   const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2
 183   const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3
 184   const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3
 185   const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
 186   const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
 187   const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
 188
 189   if(anAbsD1 >= anAbsD2)
 190   {
 191     if(anAbsD3 > anAbsD1)
 192     {
 193       aFoundCouple = COE13;
 194       aDetV1V2 = aDelta3;
 195     }
 196     else
 197     {
 198       aFoundCouple = COE12;
 199       aDetV1V2 = aDelta1;
 200     }
 201   }
 202   else
 203   {
 204     if(anAbsD3 > anAbsD2)
 205     {
 206       aFoundCouple = COE13;
 207       aDetV1V2 = aDelta3;
 208     }
 209     else
 210     {
 211       aFoundCouple = COE23;
 212       aDetV1V2 = aDelta2;
 213     }
 214   }
 215
 216   // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is
 217   // cross-product between directions (i.e. sine of angle).
 218   // If sine is too small then sine is (approx.) equal to angle itself.
 219   // Therefore, in this case we should compare sine with angular tolerance.
 220   // This constant is used for check if axes are parallel (see constructor
 221   // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file).
 222   if(Abs(aDetV1V2) < Precision::Angular())
 223   {
 224     throw Standard_Failure("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
 225   }
 226
 227   switch(aFoundCouple)
 228   {
 229   case COE12:
 230     break;
 231   case COE23:
 232     {
 233       math_Vector aVTemp(mVecA1);
 234       mVecA1(1) = aVTemp(2);
 235       mVecA1(2) = aVTemp(3);
 236       mVecA1(3) = aVTemp(1);
 237
 238       aVTemp = mVecA2;
 239       mVecA2(1) = aVTemp(2);
 240       mVecA2(2) = aVTemp(3);
 241       mVecA2(3) = aVTemp(1);
 242
 243       aVTemp = mVecB1;
 244       mVecB1(1) = aVTemp(2);
 245       mVecB1(2) = aVTemp(3);
 246       mVecB1(3) = aVTemp(1);
 247
 248       aVTemp = mVecB2;
 249       mVecB2(1) = aVTemp(2);
 250       mVecB2(2) = aVTemp(3);
 251       mVecB2(3) = aVTemp(1);
 252
 253       aVTemp = mVecC1;
 254       mVecC1(1) = aVTemp(2);
 255       mVecC1(2) = aVTemp(3);
 256       mVecC1(3) = aVTemp(1);
 257
 258       aVTemp = mVecC2;
 259       mVecC2(1) = aVTemp(2);
 260       mVecC2(2) = aVTemp(3);
 261       mVecC2(3) = aVTemp(1);
 262
 263       aVTemp = mVecD;
 264       mVecD(1) = aVTemp(2);
 265       mVecD(2) = aVTemp(3);
 266       mVecD(3) = aVTemp(1);
 267
 268       break;
 269     }
 270   case COE13:
 271     {
 272       math_Vector aVTemp = mVecA1;
 273       mVecA1(2) = aVTemp(3);
 274       mVecA1(3) = aVTemp(2);
 275
 276       aVTemp = mVecA2;
 277       mVecA2(2) = aVTemp(3);
 278       mVecA2(3) = aVTemp(2);
 279
 280       aVTemp = mVecB1;
 281       mVecB1(2) = aVTemp(3);
 282       mVecB1(3) = aVTemp(2);
 283
 284       aVTemp = mVecB2;
 285       mVecB2(2) = aVTemp(3);
 286       mVecB2(3) = aVTemp(2);
 287
 288       aVTemp = mVecC1;
 289       mVecC1(2) = aVTemp(3);
 290       mVecC1(3) = aVTemp(2);
 291
 292       aVTemp = mVecC2;
 293       mVecC2(2) = aVTemp(3);
 294       mVecC2(3) = aVTemp(2);
 295
 296       aVTemp = mVecD;
 297       mVecD(2) = aVTemp(3);
 298       mVecD(3) = aVTemp(2);
 299
 300       break;
 301     }
 302   default:
 303     break;
 304   }
 305
 306   //------- For V1 (begin)
 307   //sinU2
 308   mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
 309   //sinU1
 310   mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
 311   //cosU2
 312   mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
 313   //cosU1
 314   mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
 315   //Free member
 316   mM1 =  (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
 317   //------- For V1 (end)
 318
 319   //------- For V2 (begin)
 320   //sinU2
 321   mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
 322   //sinU1
 323   mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
 324   //cosU2
 325   mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
 326   //cosU1
 327   mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
 328   //Free member
 329   mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
 330   //------- For V1 (end)
 331
 332   ShortCosForm(mL11, mK11, mK1, mFIV1);
 333   ShortCosForm(mL21, mK21, mL1, mPSIV1);
 334   ShortCosForm(mL12, mK12, mK2, mFIV2);
 335   ShortCosForm(mL22, mK22, mL2, mPSIV2);
 336
 337   const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2
 338     aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2
 339     aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1
 340     aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1
 341
 342   mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2;  //Free
 343
 344   Standard_Real aA = 0.0;
 345
 346   ShortCosForm(aB2,aB1,mB,mFI1);
 347   ShortCosForm(aA2,aA1,aA,mFI2);
 348
 349   mB /= aA;
 350   mC /= aA;
 351 }
 352
 353 class WorkWithBoundaries
 354 {
 355 public:
 356   enum SearchBoundType
 357   {
 358     SearchNONE = 0,
 359     SearchV1 = 1,
 360     SearchV2 = 2
 361   };
 362
 363   struct StPInfo
 364   {
 365     StPInfo()
 366     {
 367       mySurfID = 0;
 368       myU1 = RealLast();
 369       myV1 = RealLast();
 370       myU2 = RealLast();
 371       myV2 = RealLast();
 372     }
 373
 374     //Equal to 0 for 1st surface non-zero for 2nd one.
 375     Standard_Integer mySurfID;
 376
 377     Standard_Real myU1;
 378     Standard_Real myV1;
 379     Standard_Real myU2;
 380     Standard_Real myV2;
 381
 382     bool operator>(const StPInfo& theOther) const
 383     {
 384       return myU1 > theOther.myU1;
 385     }
 386
 387     bool operator<(const StPInfo& theOther) const
 388     {
 389       return myU1 < theOther.myU1;
 390     }
 391
 392     bool operator==(const StPInfo& theOther) const
 393     {
 394       return myU1 == theOther.myU1;
 395     }
 396   };
 397
 398   WorkWithBoundaries(const IntSurf_Quadric& theQuad1,
 399                      const IntSurf_Quadric& theQuad2,
 400                      const ComputationMethods::stCoeffsValue& theCoeffs,
 401                      const Bnd_Box2d& theUVSurf1,
 402                      const Bnd_Box2d& theUVSurf2,
 403                      const Standard_Integer theNbWLines,
 404                      const Standard_Real thePeriod,
 405                      const Standard_Real theTol3D,
 406                      const Standard_Real theTol2D,
 407                      const Standard_Boolean isTheReverse) :
 408       myQuad1(theQuad1), myQuad2(theQuad2), myCoeffs(theCoeffs),
 409       myUVSurf1(theUVSurf1), myUVSurf2(theUVSurf2), myNbWLines(theNbWLines),
 410       myPeriod(thePeriod), myTol3D(theTol3D), myTol2D(theTol2D),
 411       myIsReverse(isTheReverse)
 412   {
 413   };
 414
 415   // Returns parameters of system solved while finding
 416   // intersection line
 417   const ComputationMethods::stCoeffsValue &SICoeffs() const
 418   {
 419     return myCoeffs;
 420   }
 421
 422   // Returns quadric correspond to the index theIdx.
 423   const IntSurf_Quadric& GetQSurface(const Standard_Integer theIdx) const
 424   {
 425     if (theIdx <= 1)
 426       return myQuad1;
 427
 428     return myQuad2;
 429   }
 430
 431   // Returns TRUE in case of reverting surfaces
 432   Standard_Boolean IsReversed() const
 433   {
 434     return myIsReverse;
 435   }
 436
 437   // Returns 2D-tolerance
 438   Standard_Real Get2dTolerance() const
 439   {
 440     return myTol2D;
 441   }
 442
 443   // Returns 3D-tolerance
 444   Standard_Real Get3dTolerance() const
 445   {
 446     return myTol3D;
 447   }
 448
 449   // Returns UV-bounds of 1st surface
 450   const Bnd_Box2d& UVS1() const
 451   {
 452     return myUVSurf1;
 453   }
 454
 455   // Returns UV-bounds of 2nd surface
 456   const Bnd_Box2d& UVS2() const
 457   {
 458     return myUVSurf2;
 459   }
 460
 461   void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
 462                         const Standard_Real theU1,
 463                         const Standard_Real theU1Min,
 464                         const Standard_Real theU2,
 465                         const Standard_Real theV1,
 466                         const Standard_Real theV1Prev,
 467                         const Standard_Real theV2,
 468                         const Standard_Real theV2Prev,
 469                         const Standard_Integer theWLIndex,
 470                         const Standard_Boolean theFlForce,
 471                         Standard_Boolean& isTheFound1,
 472                         Standard_Boolean& isTheFound2) const;
 473
 474   static Standard_Boolean BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
 475                                               const Standard_Real thePeriod,
 476                                               Bnd_Range theURange[]);
 477
 478   void BoundaryEstimation(const gp_Cylinder& theCy1,
 479                           const gp_Cylinder& theCy2,
 480                           Bnd_Range& theOutBoxS1,
 481                           Bnd_Range& theOutBoxS2) const;
 482
 483 protected:
 484
 485   //Solves equation (2) (see declaration of ComputationMethods class) in case,
 486   //when V1 or V2 (is set by theSBType argument) is known (corresponds to the boundary
 487   //and equal to theVzad) but U1 is unknown. Computation is made by numeric methods and
 488   //requires initial values (theVInit, theInitU2 and theInitMainVar).
 489   Standard_Boolean
 490               SearchOnVBounds(const SearchBoundType theSBType,
 491                               const Standard_Real theVzad,
 492                               const Standard_Real theVInit,
 493                               const Standard_Real theInitU2,
 494                               const Standard_Real theInitMainVar,
 495                               Standard_Real& theMainVariableValue) const;
 496
 497   const WorkWithBoundaries& operator=(const WorkWithBoundaries&);
 498
 499 private:
 500   friend class ComputationMethods;
 501
 502   const IntSurf_Quadric& myQuad1;
 503   const IntSurf_Quadric& myQuad2;
 504   const ComputationMethods::stCoeffsValue& myCoeffs;
 505   const Bnd_Box2d& myUVSurf1;
 506   const Bnd_Box2d& myUVSurf2;
 507   const Standard_Integer myNbWLines;
 508   const Standard_Real myPeriod;
 509   const Standard_Real myTol3D;
 510   const Standard_Real myTol2D;
 511   const Standard_Boolean myIsReverse;
 512 };
 513
 514 static
 515   Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
 516                                 const gp_Cone& aCo,
 517                                 IntAna_Curve& aC,
 518                                 const Standard_Real aTol,
 519                                 IntAna_ListOfCurve& aLC);
 520
 521 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
 522                                   const IntSurf_Quadric& theQuad2,
 523                                   const Handle(IntSurf_LineOn2S)& theLine,
 524                                   const ComputationMethods::stCoeffsValue& theCoeffs,
 525                                   const Standard_Integer theWLIndex,
 526                                   const Standard_Integer theMinNbPoints,
 527                                   const Standard_Integer theStartPointOnLine,
 528                                   const Standard_Integer theEndPointOnLine,
 529                                   const Standard_Real theTol2D,
 530                                   const Standard_Real thePeriodOfSurf2,
 531                                   const Standard_Boolean isTheReverse);
 532
 533 //=======================================================================
 534 //function : MinMax
 535 //purpose  : Replaces theParMIN = MIN(theParMIN, theParMAX),
 536 //                    theParMAX = MAX(theParMIN, theParMAX).
 537 //=======================================================================
 538 static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
 539 {
 540   if(theParMIN > theParMAX)
 541   {
 542     const Standard_Real aux = theParMAX;
 543     theParMAX = theParMIN;
 544     theParMIN = aux;
 545   }
 546 }
 547
 548 //=======================================================================
 549 //function : ExtremaLineLine
 550 //purpose  : Computes extrema between the given lines. Returns parameters
 551 //          on correspond curve (see correspond method for Extrema_ExtElC class).
 552 //=======================================================================
 553 static inline void ExtremaLineLine(const gp_Ax1& theC1,
 554                                    const gp_Ax1& theC2,
 555                                    const Standard_Real theCosA,
 556                                    const Standard_Real theSqSinA,
 557                                    Standard_Real& thePar1,
 558                                    Standard_Real& thePar2)
 559 {
 560   const gp_Dir &aD1 = theC1.Direction(),
 561                &aD2 = theC2.Direction();
 562
 563   const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ();
 564   const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2),
 565                       aD2L = aD2.XYZ().Dot(aL1L2);
 566
 567   thePar1 = (aD1L - theCosA * aD2L) / theSqSinA;
 568   thePar2 = (theCosA * aD1L - aD2L) / theSqSinA;
 569 }
 570
 571 //=======================================================================
 572 //function : VBoundaryPrecise
 573 //purpose  : By default, we shall consider, that V1 and V2 will be increased
 574 //            if U1 is increased. But if it is not, new V1set and/or V2set
 575 //            must be computed as [V1current - DeltaV1] (analogically
 576 //            for V2). This function processes this case.
 577 //=======================================================================
 578 static void VBoundaryPrecise( const math_Matrix& theMatr,
 579                               const Standard_Real theV1AfterDecrByDelta,
 580                               const Standard_Real theV2AfterDecrByDelta,
 581                               Standard_Real& theV1Set,
 582                               Standard_Real& theV2Set)
 583 {
 584   //Now we are going to define if V1 (and V2) increases
 585   //(or decreases) when U1 will increase.
 586   const Standard_Integer aNbDim = 3;
 587   math_Matrix aSyst(1, aNbDim, 1, aNbDim);
 588
 589   aSyst.SetCol(1, theMatr.Col(1));
 590   aSyst.SetCol(2, theMatr.Col(2));
 591   aSyst.SetCol(3, theMatr.Col(4));
 592
 593   //We have the system (see comment to StepComputing(...) function)
 594   //    {a11*dV1 + a12*dV2 + a14*dU2 = -a13*dU1
 595   //    {a21*dV1 + a22*dV2 + a24*dU2 = -a23*dU1
 596   //    {a31*dV1 + a32*dV2 + a34*dU2 = -a33*dU1
 597
 598   const Standard_Real aDet = aSyst.Determinant();
 599
 600   aSyst.SetCol(1, theMatr.Col(3));
 601   const Standard_Real aDet1 = aSyst.Determinant();
 602
 603   aSyst.SetCol(1, theMatr.Col(1));
 604   aSyst.SetCol(2, theMatr.Col(3));
 605
 606   const Standard_Real aDet2 = aSyst.Determinant();
 607
 608   //Now,
 609   //    dV1 = -dU1*aDet1/aDet
 610   //    dV2 = -dU1*aDet2/aDet
 611
 612   //If U1 is increased then dU1 > 0.
 613   //If (aDet1/aDet > 0) then dV1 < 0 and
 614   //V1 will be decreased after increasing U1.
 615
 616   //We have analogical situation with V2-parameter.
 617
 618   if(aDet*aDet1 > 0.0)
 619   {
 620     theV1Set = theV1AfterDecrByDelta;
 621   }
 622
 623   if(aDet*aDet2 > 0.0)
 624   {
 625     theV2Set = theV2AfterDecrByDelta;
 626   }
 627 }
 628
 629 //=======================================================================
 630 //function : DeltaU1Computing
 631 //purpose  : Computes new step for U1 parameter.
 632 //=======================================================================
 633 static inline
 634         Standard_Boolean DeltaU1Computing(const math_Matrix& theSyst,
 635                                           const math_Vector& theFree,
 636                                           Standard_Real& theDeltaU1Found)
 637 {
 638   Standard_Real aDet = theSyst.Determinant();
 639
 640   if(Abs(aDet) > aNulValue)
 641   {
 642     math_Matrix aSyst1(theSyst);
 643     aSyst1.SetCol(2, theFree);
 644
 645     theDeltaU1Found = Abs(aSyst1.Determinant()/aDet);
 646     return Standard_True;
 647   }
 648
 649   return Standard_False;
 650 }
 651
 652 //=======================================================================
 653 //function : StepComputing
 654 //purpose  :
 655 //
 656 //Attention!!!:
 657 //            theMatr must have 3*5-dimension strictly.
 658 //            For system
 659 //                {a11*V1+a12*V2+a13*dU1+a14*dU2=b1;
 660 //                {a21*V1+a22*V2+a23*dU1+a24*dU2=b2;
 661 //                {a31*V1+a32*V2+a33*dU1+a34*dU2=b3;
 662 //            theMatr must be following:
 663 //                (a11 a12 a13 a14 b1)
 664 //                (a21 a22 a23 a24 b2)
 665 //                (a31 a32 a33 a34 b3)
 666 //=======================================================================
 667 static Standard_Boolean StepComputing(const math_Matrix& theMatr,
 668                                       const Standard_Real theV1Cur,
 669                                       const Standard_Real theV2Cur,
 670                                       const Standard_Real theDeltaV1,
 671                                       const Standard_Real theDeltaV2,
 672                                       Standard_Real& theDeltaU1Found/*,
 673                                       Standard_Real& theDeltaU2Found,
 674                                       Standard_Real& theV1Found,
 675                                       Standard_Real& theV2Found*/)
 676 {
 677 #ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
 678   bool flShow = false;
 679
 680   if(flShow)
 681   {
 682     printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
 683               theMatr(1,1), theMatr(1,2), theMatr(1,3), theMatr(1,4), theMatr(1,5));
 684     printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
 685               theMatr(2,1), theMatr(2,2), theMatr(2,3), theMatr(2,4), theMatr(2,5));
 686     printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
 687               theMatr(3,1), theMatr(3,2), theMatr(3,3), theMatr(3,4), theMatr(3,5));
 688   }
 689 #endif
 690
 691   Standard_Boolean isSuccess = Standard_False;
 692   theDeltaU1Found/* = theDeltaU2Found*/ = RealLast();
 693   //theV1Found = theV1set;
 694   //theV2Found = theV2Set;
 695   const Standard_Integer aNbDim = 3;
 696
 697   math_Matrix aSyst(1, aNbDim, 1, aNbDim);
 698   math_Vector aFree(1, aNbDim);
 699
 700   //By default, increasing V1(U1) and V2(U1) functions is
 701   //considered
 702   Standard_Real aV1Set = theV1Cur + theDeltaV1,
 703                 aV2Set = theV2Cur + theDeltaV2;
 704
 705   //However, what is indeed?
 706   VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1,
 707                     theV2Cur - theDeltaV2, aV1Set, aV2Set);
 708
 709   aSyst.SetCol(2, theMatr.Col(3));
 710   aSyst.SetCol(3, theMatr.Col(4));
 711
 712   for(Standard_Integer i = 0; i < 2; i++)
 713   {
 714     if(i == 0)
 715     {//V1 is known
 716       aSyst.SetCol(1, theMatr.Col(2));
 717       aFree.Set(1, aNbDim, theMatr.Col(5)-aV1Set*theMatr.Col(1));
 718     }
 719     else
 720     {//i==1 => V2 is known
 721       aSyst.SetCol(1, theMatr.Col(1));
 722       aFree.Set(1, aNbDim, theMatr.Col(5)-aV2Set*theMatr.Col(2));
 723     }
 724
 725     Standard_Real aNewDU = theDeltaU1Found;
 726     if(DeltaU1Computing(aSyst, aFree, aNewDU))
 727     {
 728       isSuccess = Standard_True;
 729       if(aNewDU < theDeltaU1Found)
 730       {
 731         theDeltaU1Found = aNewDU;
 732       }
 733     }
 734   }
 735
 736   if(!isSuccess)
 737   {
 738     aFree = theMatr.Col(5) - aV1Set*theMatr.Col(1) - aV2Set*theMatr.Col(2);
 739     math_Matrix aSyst1(1, aNbDim, 1, 2);
 740     aSyst1.SetCol(1, aSyst.Col(2));
 741     aSyst1.SetCol(2, aSyst.Col(3));
 742
 743     //Now we have overdetermined system.
 744
 745     const Standard_Real aDet1 = theMatr(1,3)*theMatr(2,4) - theMatr(2,3)*theMatr(1,4);
 746     const Standard_Real aDet2 = theMatr(1,3)*theMatr(3,4) - theMatr(3,3)*theMatr(1,4);
 747     const Standard_Real aDet3 = theMatr(2,3)*theMatr(3,4) - theMatr(3,3)*theMatr(2,4);
 748     const Standard_Real anAbsD1 = Abs(aDet1);
 749     const Standard_Real anAbsD2 = Abs(aDet2);
 750     const Standard_Real anAbsD3 = Abs(aDet3);
 751
 752     if(anAbsD1 >= anAbsD2)
 753     {
 754       if(anAbsD1 >= anAbsD3)
 755       {
 756         //Det1
 757         if(anAbsD1 <= aNulValue)
 758           return isSuccess;
 759
 760         theDeltaU1Found = Abs(aFree(1)*theMatr(2,4) - aFree(2)*theMatr(1,4))/anAbsD1;
 761         isSuccess = Standard_True;
 762       }
 763       else
 764       {
 765         //Det3
 766         if(anAbsD3 <= aNulValue)
 767           return isSuccess;
 768
 769         theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
 770         isSuccess = Standard_True;
 771       }
 772     }
 773     else
 774     {
 775       if(anAbsD2 >= anAbsD3)
 776       {
 777         //Det2
 778         if(anAbsD2 <= aNulValue)
 779           return isSuccess;
 780
 781         theDeltaU1Found = Abs(aFree(1)*theMatr(3,4) - aFree(3)*theMatr(1,4))/anAbsD2;
 782         isSuccess = Standard_True;
 783       }
 784       else
 785       {
 786         //Det3
 787         if(anAbsD3 <= aNulValue)
 788           return isSuccess;
 789
 790         theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
 791         isSuccess = Standard_True;
 792       }
 793     }
 794   }
 795
 796   return isSuccess;
 797 }
 798
 799 //=======================================================================
 800 //function : ProcessBounds
 801 //purpose  :
 802 //=======================================================================
 803 void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne courante
 804                    const IntPatch_SequenceOfLine& slin,
 805                    const IntSurf_Quadric& Quad1,
 806                    const IntSurf_Quadric& Quad2,
 807                    Standard_Boolean& procf,
 808                    const gp_Pnt& ptf,                     //-- Debut Ligne Courante
 809                    const Standard_Real first,             //-- Paramf
 810                    Standard_Boolean& procl,
 811                    const gp_Pnt& ptl,                     //-- Fin Ligne courante
 812                    const Standard_Real last,              //-- Paraml
 813                    Standard_Boolean& Multpoint,
 814                    const Standard_Real Tol)
 815 {
 816   Standard_Integer j,k;
 817   Standard_Real U1,V1,U2,V2;
 818   IntPatch_Point ptsol;
 819   Standard_Real d;
 820
 821   if (procf && procl) {
 822     j = slin.Length() + 1;
 823   }
 824   else {
 825     j = 1;
 826   }
 827
 828
 829   //-- On prend les lignes deja enregistrees
 830
 831   while (j <= slin.Length()) {
 832     if(slin.Value(j)->ArcType() == IntPatch_Analytic) {
 833       const Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(j));
 834       k = 1;
 835
 836       //-- On prend les vertex des lignes deja enregistrees
 837
 838       while (k <= aligold->NbVertex()) {
 839         ptsol = aligold->Vertex(k);
 840         if (!procf) {
 841           d=ptf.Distance(ptsol.Value());
 842           if (d <= Tol) {
 843             if (!ptsol.IsMultiple()) {
 844               //-- le point ptsol (de aligold) est declare multiple sur aligold
 845               Multpoint = Standard_True;
 846               ptsol.SetMultiple(Standard_True);
 847               aligold->Replace(k,ptsol);
 848             }
 849             ptsol.SetParameter(first);
 850             alig->AddVertex(ptsol);
 851             alig->SetFirstPoint(alig->NbVertex());
 852             procf = Standard_True;
 853
 854             //-- On restore le point avec son parametre sur aligold
 855             ptsol = aligold->Vertex(k);
 856
 857           }
 858         }
 859         if (!procl) {
 860           if (ptl.Distance(ptsol.Value()) <= Tol) {
 861             if (!ptsol.IsMultiple()) {
 862               Multpoint = Standard_True;
 863               ptsol.SetMultiple(Standard_True);
 864               aligold->Replace(k,ptsol);
 865             }
 866             ptsol.SetParameter(last);
 867             alig->AddVertex(ptsol);
 868             alig->SetLastPoint(alig->NbVertex());
 869             procl = Standard_True;
 870
 871             //-- On restore le point avec son parametre sur aligold
 872             ptsol = aligold->Vertex(k);
 873
 874           }
 875         }
 876         if (procf && procl) {
 877           k = aligold->NbVertex()+1;
 878         }
 879         else {
 880           k = k+1;
 881         }
 882       }
 883       if (procf && procl) {
 884         j = slin.Length()+1;
 885       }
 886       else {
 887         j = j+1;
 888       }
 889     }
 890   }
 891   if (!procf && !procl) {
 892     Quad1.Parameters(ptf,U1,V1);
 893     Quad2.Parameters(ptf,U2,V2);
 894     ptsol.SetValue(ptf,Tol,Standard_False);
 895     ptsol.SetParameters(U1,V1,U2,V2);
 896     ptsol.SetParameter(first);
 897     if (ptf.Distance(ptl) <= Tol) {
 898       ptsol.SetMultiple(Standard_True); // a voir
 899       Multpoint = Standard_True;        // a voir de meme
 900       alig->AddVertex(ptsol);
 901       alig->SetFirstPoint(alig->NbVertex());
 902
 903       ptsol.SetParameter(last);
 904       alig->AddVertex(ptsol);
 905       alig->SetLastPoint(alig->NbVertex());
 906     }
 907     else {
 908       alig->AddVertex(ptsol);
 909       alig->SetFirstPoint(alig->NbVertex());
 910       Quad1.Parameters(ptl,U1,V1);
 911       Quad2.Parameters(ptl,U2,V2);
 912       ptsol.SetValue(ptl,Tol,Standard_False);
 913       ptsol.SetParameters(U1,V1,U2,V2);
 914       ptsol.SetParameter(last);
 915       alig->AddVertex(ptsol);
 916       alig->SetLastPoint(alig->NbVertex());
 917     }
 918   }
 919   else if (!procf) {
 920     Quad1.Parameters(ptf,U1,V1);
 921     Quad2.Parameters(ptf,U2,V2);
 922     ptsol.SetValue(ptf,Tol,Standard_False);
 923     ptsol.SetParameters(U1,V1,U2,V2);
 924     ptsol.SetParameter(first);
 925     alig->AddVertex(ptsol);
 926     alig->SetFirstPoint(alig->NbVertex());
 927   }
 928   else if (!procl) {
 929     Quad1.Parameters(ptl,U1,V1);
 930     Quad2.Parameters(ptl,U2,V2);
 931     ptsol.SetValue(ptl,Tol,Standard_False);
 932     ptsol.SetParameters(U1,V1,U2,V2);
 933     ptsol.SetParameter(last);
 934     alig->AddVertex(ptsol);
 935     alig->SetLastPoint(alig->NbVertex());
 936   }
 937 }
 938
 939 //=======================================================================
 940 //function : CyCyAnalyticalIntersect
 941 //purpose  : Checks if intersection curve is analytical (line, circle, ellipse)
 942 //            and returns these curves.
 943 //=======================================================================
 944 Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1,
 945                                           const IntSurf_Quadric& Quad2,
 946                                           const IntAna_QuadQuadGeo& theInter,
 947                                           const Standard_Real Tol,
 948                                           Standard_Boolean& Empty,
 949                                           Standard_Boolean& Same,
 950                                           Standard_Boolean& Multpoint,
 951                                           IntPatch_SequenceOfLine& slin,
 952                                           IntPatch_SequenceOfPoint& spnt)
 953 {
 954   IntPatch_Point ptsol;
 955
 956   Standard_Integer i;
 957
 958   IntSurf_TypeTrans trans1,trans2;
 959   IntAna_ResultType typint;
 960
 961   gp_Elips elipsol;
 962   gp_Lin linsol;
 963
 964   gp_Cylinder Cy1(Quad1.Cylinder());
 965   gp_Cylinder Cy2(Quad2.Cylinder());
 966
 967   typint = theInter.TypeInter();
 968   Standard_Integer NbSol = theInter.NbSolutions();
 969   Empty = Standard_False;
 970   Same  = Standard_False;
 971
 972   switch (typint)
 973       {
 974   case IntAna_Empty:
 975     {
 976       Empty = Standard_True;
 977       }
 978     break;
 979
 980   case IntAna_Same:
 981       {
 982       Same  = Standard_True;
 983       }
 984     break;
 985
 986   case IntAna_Point:
 987     {
 988       gp_Pnt psol(theInter.Point(1));
 989       ptsol.SetValue(psol,Tol,Standard_True);
 990
 991       Standard_Real U1,V1,U2,V2;
 992       Quad1.Parameters(psol, U1, V1);
 993       Quad2.Parameters(psol, U2, V2);
 994
 995       ptsol.SetParameters(U1,V1,U2,V2);
 996       spnt.Append(ptsol);
 997     }
 998     break;
 999
1000   case IntAna_Line:
1001     {
1002       gp_Pnt ptref;
1003       if (NbSol == 1)
1004       { // Cylinders are tangent to each other by line
1005         linsol = theInter.Line(1);
1006         ptref = linsol.Location();
1007
1008         //Radius-vectors
1009         gp_Dir crb1(gp_Vec(ptref,Cy1.Location()));
1010         gp_Dir crb2(gp_Vec(ptref,Cy2.Location()));
1011
1012         //outer normal lines
1013         gp_Vec norm1(Quad1.Normale(ptref));
1014         gp_Vec norm2(Quad2.Normale(ptref));
1015         IntSurf_Situation situcyl1;
1016         IntSurf_Situation situcyl2;
1017
1018         if (crb1.Dot(crb2) < 0.)
1019         { // centre de courbures "opposes"
1020             //ATTENTION!!!
1021             //        Normal and Radius-vector of the 1st(!) cylinder
1022             //        is used for judging what the situation of the 2nd(!)
1023             //        cylinder is.
1024
1025           if (norm1.Dot(crb1) > 0.)
1026           {
1027             situcyl2 = IntSurf_Inside;
1028           }
1029           else
1030           {
1031             situcyl2 = IntSurf_Outside;
1032           }
1033
1034           if (norm2.Dot(crb2) > 0.)
1035           {
1036             situcyl1 = IntSurf_Inside;
1037           }
1038           else
1039           {
1040             situcyl1 = IntSurf_Outside;
1041           }
1042         }
1043         else
1044           {
1045           if (Cy1.Radius() < Cy2.Radius())
1046           {
1047             if (norm1.Dot(crb1) > 0.)
1048             {
1049               situcyl2 = IntSurf_Inside;
1050             }
1051             else
1052             {
1053               situcyl2 = IntSurf_Outside;
1054             }
1055
1056             if (norm2.Dot(crb2) > 0.)
1057             {
1058               situcyl1 = IntSurf_Outside;
1059             }
1060             else
1061             {
1062               situcyl1 = IntSurf_Inside;
1063             }
1064           }
1065           else
1066           {
1067             if (norm1.Dot(crb1) > 0.)
1068             {
1069               situcyl2 = IntSurf_Outside;
1070             }
1071             else
1072             {
1073               situcyl2 = IntSurf_Inside;
1074             }
1075
1076             if (norm2.Dot(crb2) > 0.)
1077             {
1078               situcyl1 = IntSurf_Inside;
1079             }
1080             else
1081             {
1082               situcyl1 = IntSurf_Outside;
1083             }
1084           }
1085         }
1086
1087         Handle(IntPatch_GLine) glig =  new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2);
1088         slin.Append(glig);
1089       }
1090       else
1091       {
1092         for (i=1; i <= NbSol; i++)
1093         {
1094           linsol = theInter.Line(i);
1095           ptref = linsol.Location();
1096           gp_Vec lsd = linsol.Direction();
1097
1098           //Theoretically, qwe = +/- 1.0.
1099           Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
1100           if (qwe >0.00000001)
1101           {
1102             trans1 = IntSurf_Out;
1103             trans2 = IntSurf_In;
1104           }
1105           else if (qwe <-0.00000001)
1106           {
1107             trans1 = IntSurf_In;
1108             trans2 = IntSurf_Out;
1109           }
1110           else
1111           {
1112             trans1=trans2=IntSurf_Undecided;
1113           }
1114
1115           Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2);
1116           slin.Append(glig);
1117         }
1118       }
1119     }
1120     break;
1121
1122   case IntAna_Ellipse:
1123     {
1124       gp_Vec Tgt;
1125       gp_Pnt ptref;
1126       IntPatch_Point pmult1, pmult2;
1127
1128       elipsol = theInter.Ellipse(1);
1129
1130       gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
1131       gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));
1132
1133       Multpoint = Standard_True;
1134       pmult1.SetValue(pttang1,Tol,Standard_True);
1135       pmult2.SetValue(pttang2,Tol,Standard_True);
1136       pmult1.SetMultiple(Standard_True);
1137       pmult2.SetMultiple(Standard_True);
1138
1139       Standard_Real oU1,oV1,oU2,oV2;
1140       Quad1.Parameters(pttang1, oU1, oV1);
1141       Quad2.Parameters(pttang1, oU2, oV2);
1142
1143       pmult1.SetParameters(oU1,oV1,oU2,oV2);
1144       Quad1.Parameters(pttang2,oU1,oV1);
1145       Quad2.Parameters(pttang2,oU2,oV2);
1146
1147       pmult2.SetParameters(oU1,oV1,oU2,oV2);
1148
1149       // on traite la premiere ellipse
1150
1151       //-- Calcul de la Transition de la ligne
1152       ElCLib::D1(0.,elipsol,ptref,Tgt);
1153
1154       //Theoretically, qwe = +/- |Tgt|.
1155       Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
1156       if (qwe>0.00000001)
1157       {
1158         trans1 = IntSurf_Out;
1159         trans2 = IntSurf_In;
1160       }
1161       else if (qwe<-0.00000001)
1162       {
1163         trans1 = IntSurf_In;
1164         trans2 = IntSurf_Out;
1165       }
1166       else
1167       {
1168         trans1=trans2=IntSurf_Undecided;
1169       }
1170
1171       //-- Transition calculee au point 0 -> Trans2 , Trans1
1172       //-- car ici, on devarit calculer en PI
1173       Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans2,trans1);
1174       //
1175       {
1176         Standard_Real aU1, aV1, aU2, aV2;
1177         IntPatch_Point aIP;
1178         gp_Pnt aP (ElCLib::Value(0., elipsol));
1179         //
1180         aIP.SetValue(aP,Tol,Standard_False);
1181         aIP.SetMultiple(Standard_False);
1182         //
1183         Quad1.Parameters(aP, aU1, aV1);
1184         Quad2.Parameters(aP, aU2, aV2);
1185
1186         aIP.SetParameters(aU1, aV1, aU2, aV2);
1187         //
1188         aIP.SetParameter(0.);
1189         glig->AddVertex(aIP);
1190         glig->SetFirstPoint(1);
1191         //
1192         aIP.SetParameter(2.*M_PI);
1193         glig->AddVertex(aIP);
1194         glig->SetLastPoint(2);
1195       }
1196       //
1197       pmult1.SetParameter(0.5*M_PI);
1198       glig->AddVertex(pmult1);
1199       //
1200       pmult2.SetParameter(1.5*M_PI);
1201       glig->AddVertex(pmult2);
1202
1203       //
1204       slin.Append(glig);
1205
1206       //-- Transitions calculee au point 0    OK
1207       //
1208       // on traite la deuxieme ellipse
1209       elipsol = theInter.Ellipse(2);
1210
1211       Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
1212       Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
1213       Standard_Real parampourtransition = 0.0;
1214       if (param1 < param2)
1215       {
1216         pmult1.SetParameter(0.5*M_PI);
1217         pmult2.SetParameter(1.5*M_PI);
1218         parampourtransition = M_PI;
1219       }
1220       else {
1221         pmult1.SetParameter(1.5*M_PI);
1222         pmult2.SetParameter(0.5*M_PI);
1223         parampourtransition = 0.0;
1224       }
1225
1226       //-- Calcul des transitions de ligne pour la premiere ligne
1227       ElCLib::D1(parampourtransition,elipsol,ptref,Tgt);
1228
1229       //Theoretically, qwe = +/- |Tgt|.
1230       qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
1231       if(qwe> 0.00000001)
1232       {
1233         trans1 = IntSurf_Out;
1234         trans2 = IntSurf_In;
1235       }
1236       else if(qwe< -0.00000001)
1237       {
1238         trans1 = IntSurf_In;
1239         trans2 = IntSurf_Out;
1240       }
1241       else
1242       {
1243         trans1=trans2=IntSurf_Undecided;
1244       }
1245
1246       //-- La transition a ete calculee sur un point de cette ligne
1247       glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2);
1248       //
1249       {
1250         Standard_Real aU1, aV1, aU2, aV2;
1251         IntPatch_Point aIP;
1252         gp_Pnt aP (ElCLib::Value(0., elipsol));
1253         //
1254         aIP.SetValue(aP,Tol,Standard_False);
1255         aIP.SetMultiple(Standard_False);
1256         //
1257
1258         Quad1.Parameters(aP, aU1, aV1);
1259         Quad2.Parameters(aP, aU2, aV2);
1260
1261         aIP.SetParameters(aU1, aV1, aU2, aV2);
1262         //
1263         aIP.SetParameter(0.);
1264         glig->AddVertex(aIP);
1265         glig->SetFirstPoint(1);
1266         //
1267         aIP.SetParameter(2.*M_PI);
1268         glig->AddVertex(aIP);
1269         glig->SetLastPoint(2);
1270       }
1271       //
1272       glig->AddVertex(pmult1);
1273       glig->AddVertex(pmult2);
1274       //
1275       slin.Append(glig);
1276     }
1277     break;
1278
1279   case IntAna_Parabola:
1280   case IntAna_Hyperbola:
1281     throw Standard_Failure("IntCyCy(): Wrong intersection type!");
1282
1283   case IntAna_Circle:
1284     // Circle is useful when we will work with trimmed surfaces
1285     // (two cylinders can be tangent by their basises, e.g. circle)
1286   case IntAna_NoGeometricSolution:
1287   default:
1288     return Standard_False;
1289   }
1290
1291   return Standard_True;
1292 }
1293
1294 //=======================================================================
1295 //function : ShortCosForm
1296 //purpose  : Represents theCosFactor*cosA+theSinFactor*sinA as
1297 //            theCoeff*cos(A-theAngle) if it is possibly (all angles are
1298 //            in radians).
1299 //=======================================================================
1300 static void ShortCosForm( const Standard_Real theCosFactor,
1301                           const Standard_Real theSinFactor,
1302                           Standard_Real& theCoeff,
1303                           Standard_Real& theAngle)
1304 {
1305   theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
1306   theAngle = 0.0;
1307   if(IsEqual(theCoeff, 0.0))
1308   {
1309     theAngle = 0.0;
1310     return;
1311   }
1312
1313   theAngle = acos(Abs(theCosFactor/theCoeff));
1314
1315   if(theSinFactor > 0.0)
1316   {
1317     if(IsEqual(theCosFactor, 0.0))
1318     {
1319       theAngle = M_PI/2.0;
1320     }
1321     else if(theCosFactor < 0.0)
1322     {
1323       theAngle = M_PI-theAngle;
1324     }
1325   }
1326   else if(IsEqual(theSinFactor, 0.0))
1327   {
1328     if(theCosFactor < 0.0)
1329     {
1330       theAngle = M_PI;
1331     }
1332   }
1333   if(theSinFactor < 0.0)
1334   {
1335     if(theCosFactor > 0.0)
1336     {
1337       theAngle = 2.0*M_PI-theAngle;
1338     }
1339     else if(IsEqual(theCosFactor, 0.0))
1340     {
1341       theAngle = 3.0*M_PI/2.0;
1342     }
1343     else if(theCosFactor < 0.0)
1344     {
1345       theAngle = M_PI+theAngle;
1346     }
1347   }
1348 }
1349
1350 //=======================================================================
1351 //function : CylCylMonotonicity
1352 //purpose  : Determines, if U2(U1) function is increasing.
1353 //=======================================================================
1354 Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par,
1355                                                         const Standard_Integer theWLIndex,
1356                                                         const stCoeffsValue& theCoeffs,
1357                                                         const Standard_Real thePeriod,
1358                                                         Standard_Boolean& theIsIncreasing)
1359 {
1360   // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
1361   //Accordingly,
1362   //Func. U2(X1) = FI2 + X1;
1363   //Func. X1(X2) = anArccosFactor*X2;
1364   //Func. X2(X3) = acos(X3);
1365   //Func. X3(X4) = B*X4 + C;
1366   //Func. X4(U1) = cos(U1 - FI1).
1367   //
1368   //Consequently,
1369   //U2(X1) is always increasing.
1370   //X1(X2) is increasing if anArccosFactor > 0.0 and
1371   //is decreasing otherwise.
1372   //X2(X3) is always decreasing.
1373   //Therefore, U2(X3) is decreasing if anArccosFactor > 0.0 and
1374   //is increasing otherwise.
1375   //X3(X4) is increasing if B > 0 and is decreasing otherwise.
1376   //X4(U1) is increasing if U1 - FI1 in [PI, 2*PI) and
1377   //is decreasing U1 - FI1 in [0, PI) or if (U1 - FI1 == 2PI).
1378   //After that, we can predict behaviour of U2(U1) function:
1379   //if it is increasing or decreasing.
1380
1381   //For "+/-" sign. If isPlus == TRUE, "+" is chosen, otherwise, we choose "-".
1382   Standard_Boolean isPlus = Standard_False;
1383
1384   switch(theWLIndex)
1385   {
1386   case 0:
1387     isPlus = Standard_True;
1388     break;
1389   case 1:
1390     isPlus = Standard_False;
1391     break;
1392   default:
1393     //throw Standard_Failure("Error. Range Error!!!!");
1394     return Standard_False;
1395   }
1396
1397   Standard_Real aU1Temp = theU1par - theCoeffs.mFI1;
1398   InscribePoint(0, thePeriod, aU1Temp, 0.0, thePeriod, Standard_False);
1399
1400   theIsIncreasing = Standard_True;
1401
1402   if(((M_PI - aU1Temp) < RealSmall()) && (aU1Temp < thePeriod))
1403   {
1404     theIsIncreasing = Standard_False;
1405   }
1406
1407   if(theCoeffs.mB < 0.0)
1408   {
1409     theIsIncreasing = !theIsIncreasing;
1410   }
1411
1412   if(!isPlus)
1413   {
1414     theIsIncreasing = !theIsIncreasing;
1415   }
1416
1417   return Standard_True;
1418 }
1419
1420 //=======================================================================
1421 //function : CylCylComputeParameters
1422 //purpose  : Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
1423 //            estimates the tolerance of U2-computing (estimation result is
1424 //            assigned to *theDelta value).
1425 //=======================================================================
1426 Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
1427                                                 const Standard_Integer theWLIndex,
1428                                                 const stCoeffsValue& theCoeffs,
1429                                                 Standard_Real& theU2,
1430                                                 Standard_Real* const theDelta)
1431 {
1432   //This formula is got from some experience and can be changed.
1433   const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue);
1434   const Standard_Real aTol = 1.0 - aTol0;
1435
1436   if(theWLIndex < 0 || theWLIndex > 1)
1437     return Standard_False;
1438
1439   const Standard_Real aSign = theWLIndex ? -1.0 : 1.0;
1440
1441   Standard_Real anArg = cos(theU1par - theCoeffs.mFI1);
1442   anArg = theCoeffs.mB*anArg + theCoeffs.mC;
1443
1444   if(anArg >= aTol)
1445   {
1446     if(theDelta)
1447       *theDelta = 0.0;
1448
1449     anArg = 1.0;
1450   }
1451   else if(anArg <= -aTol)
1452   {
1453     if(theDelta)
1454       *theDelta = 0.0;
1455
1456     anArg = -1.0;
1457   }
1458   else if(theDelta)
1459   {
1460     //There is a case, when
1461     //  const double aPar = 0.99999999999721167;
1462     //  const double aFI2 = 2.3593296083566181e-006;
1463
1464     //Then
1465     //  aPar - cos(aFI2) == -5.10703e-015 ==> cos(aFI2) == aPar.
1466     //Theoreticaly, in this case
1467     //  aFI2 == +/- acos(aPar).
1468     //However,
1469     //  acos(aPar) - aFI2 == 2.16362e-009.
1470     //Error is quite big.
1471
1472     //This error should be estimated. Let use following way, which is based
1473     //on expanding to Taylor series.
1474
1475     //  acos(p)-acos(p+x) = x/sqrt(1-p*p).
1476
1477     //If p == (1-d) (when p > 0) or p == (-1+d) (when p < 0) then
1478     //  acos(p)-acos(p+x) = x/sqrt(d*(2-d)).
1479
1480     //Here always aTol0 <= d <= 1. Max(x) is considered (!) to be equal to aTol0.
1481     //In this case
1482     //  8*aTol0 <= acos(p)-acos(p+x) <= sqrt(2/(2-aTol0)-1),
1483     //                                              because 0 < aTol0 < 1.
1484     //E.g. when aTol0 = 1.0e-11,
1485     //   8.0e-11 <= acos(p)-acos(p+x) < 2.24e-6.
1486
1487     const Standard_Real aDelta = Min(1.0-anArg, 1.0+anArg);
1488     Standard_DivideByZero_Raise_if((aDelta*aDelta < RealSmall()) || (aDelta >= 2.0),
1489           "IntPatch_ImpImpIntersection_4.gxx, CylCylComputeParameters()");
1490     *theDelta = aTol0/sqrt(aDelta*(2.0-aDelta));
1491   }
1492
1493   theU2 = acos(anArg);
1494   theU2 = theCoeffs.mFI2 + aSign*theU2;
1495
1496   return Standard_True;
1497 }
1498
1499 //=======================================================================
1500 //function : CylCylComputeParameters
1501 //purpose  : Computes V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
1502 //=======================================================================
1503 Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1,
1504                                                 const Standard_Real theU2,
1505                                                 const stCoeffsValue& theCoeffs,
1506                                                 Standard_Real& theV1,
1507                                                 Standard_Real& theV2)
1508 {
1509   theV1 = theCoeffs.mK21 * sin(theU2) +
1510           theCoeffs.mK11 * sin(theU1) +
1511           theCoeffs.mL21 * cos(theU2) +
1512           theCoeffs.mL11 * cos(theU1) + theCoeffs.mM1;
1513
1514   theV2 = theCoeffs.mK22 * sin(theU2) +
1515           theCoeffs.mK12 * sin(theU1) +
1516           theCoeffs.mL22 * cos(theU2) +
1517           theCoeffs.mL12 * cos(theU1) + theCoeffs.mM2;
1518
1519   return Standard_True;
1520 }
1521
1522 //=======================================================================
1523 //function : CylCylComputeParameters
1524 //purpose  : Computes U2 (U-parameter of the 2nd cylinder),
1525 //            V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
1526 //=======================================================================
1527 Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
1528                                                 const Standard_Integer theWLIndex,
1529                                                 const stCoeffsValue& theCoeffs,
1530                                                 Standard_Real& theU2,
1531                                                 Standard_Real& theV1,
1532                                                 Standard_Real& theV2)
1533 {
1534   if(!CylCylComputeParameters(theU1par, theWLIndex, theCoeffs, theU2))
1535     return Standard_False;
1536
1537   if(!CylCylComputeParameters(theU1par, theU2, theCoeffs, theV1, theV2))
1538     return Standard_False;
1539
1540   return Standard_True;
1541 }
1542
1543 //=======================================================================
1544 //function : SearchOnVBounds
1545 //purpose  :
1546 //=======================================================================
1547 Standard_Boolean WorkWithBoundaries::
1548                         SearchOnVBounds(const SearchBoundType theSBType,
1549                                         const Standard_Real theVzad,
1550                                         const Standard_Real theVInit,
1551                                         const Standard_Real theInitU2,
1552                                         const Standard_Real theInitMainVar,
1553                                         Standard_Real& theMainVariableValue) const
1554 {
1555   const Standard_Integer aNbDim = 3;
1556   const Standard_Real aMaxError = 4.0*M_PI; // two periods
1557
1558   theMainVariableValue = theInitMainVar;
1559   const Standard_Real aTol2 = 1.0e-18;
1560   Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit;
1561
1562   //Structure of aMatr:
1563   //  C_{1}*U_{1} & C_{2}*U_{2} & C_{3}*V_{*},
1564   //where C_{1}, C_{2} and C_{3} are math_Vector.
1565   math_Matrix aMatr(1, aNbDim, 1, aNbDim);
1566
1567   Standard_Real anError = RealLast();
1568   Standard_Real anErrorPrev = anError;
1569   Standard_Integer aNbIter = 0;
1570   do
1571   {
1572     if(++aNbIter > 1000)
1573       return Standard_False;
1574
1575     const Standard_Real aSinU1 = sin(aMainVarPrev),
1576                         aCosU1 = cos(aMainVarPrev),
1577                         aSinU2 = sin(aU2Prev),
1578                         aCosU2 = cos(aU2Prev);
1579
1580     math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev +
1581                                               myCoeffs.mVecB2) * aSinU2 -
1582                               (myCoeffs.mVecB2 * aU2Prev -
1583                                               myCoeffs.mVecA2) * aCosU2 +
1584                               (myCoeffs.mVecA1 * aMainVarPrev +
1585                                               myCoeffs.mVecB1) * aSinU1 -
1586                               (myCoeffs.mVecB1 * aMainVarPrev -
1587                                               myCoeffs.mVecA1) * aCosU1 +
1588                                                             myCoeffs.mVecD;
1589
1590     math_Vector aMSum(1, 3);
1591
1592     switch(theSBType)
1593     {
1594     case SearchV1:
1595       aMatr.SetCol(3, myCoeffs.mVecC2);
1596       aMSum = myCoeffs.mVecC1 * theVzad;
1597       aVecFreeMem -= aMSum;
1598       aMSum += myCoeffs.mVecC2*anOtherVar;
1599       break;
1600
1601     case SearchV2:
1602       aMatr.SetCol(3, myCoeffs.mVecC1);
1603       aMSum = myCoeffs.mVecC2 * theVzad;
1604       aVecFreeMem -= aMSum;
1605       aMSum += myCoeffs.mVecC1*anOtherVar;
1606       break;
1607
1608     default:
1609       return Standard_False;
1610     }
1611
1612     aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1);
1613     aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2);
1614
1615     Standard_Real aDetMainSyst = aMatr.Determinant();
1616
1617     if(Abs(aDetMainSyst) < aNulValue)
1618     {
1619       return Standard_False;
1620     }
1621
1622     math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr);
1623     aM1.SetCol(1, aVecFreeMem);
1624     aM2.SetCol(2, aVecFreeMem);
1625     aM3.SetCol(3, aVecFreeMem);
1626
1627     const Standard_Real aDetMainVar = aM1.Determinant();
1628     const Standard_Real aDetVar1    = aM2.Determinant();
1629     const Standard_Real aDetVar2    = aM3.Determinant();
1630
1631     Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
1632
1633     if(Abs(aDelta) > aMaxError)
1634       return Standard_False;
1635
1636     anError = aDelta*aDelta;
1637     aMainVarPrev += aDelta;
1638
1639     ///
1640     aDelta = aDetVar1/aDetMainSyst-aU2Prev;
1641
1642     if(Abs(aDelta) > aMaxError)
1643       return Standard_False;
1644
1645     anError += aDelta*aDelta;
1646     aU2Prev += aDelta;
1647
1648     ///
1649     aDelta = aDetVar2/aDetMainSyst-anOtherVar;
1650     anError += aDelta*aDelta;
1651     anOtherVar += aDelta;
1652
1653     if(anError > anErrorPrev)
1654     {//Method diverges. Keep the best result
1655       const Standard_Real aSinU1Last = sin(aMainVarPrev),
1656                           aCosU1Last = cos(aMainVarPrev),
1657                           aSinU2Last = sin(aU2Prev),
1658                           aCosU2Last = cos(aU2Prev);
1659       aMSum -= (myCoeffs.mVecA1*aCosU1Last +
1660                 myCoeffs.mVecB1*aSinU1Last +
1661                 myCoeffs.mVecA2*aCosU2Last +
1662                 myCoeffs.mVecB2*aSinU2Last +
1663                 myCoeffs.mVecD);
1664       const Standard_Real aSQNorm = aMSum.Norm2();
1665       return (aSQNorm < aTol2);
1666     }
1667     else
1668     {
1669       theMainVariableValue = aMainVarPrev;
1670     }
1671
1672     anErrorPrev = anError;
1673   }
1674   while(anError > aTol2);
1675
1676   theMainVariableValue = aMainVarPrev;
1677
1678   return Standard_True;
1679 }
1680
1681 //=======================================================================
1682 //function : InscribePoint
1683 //purpose  : If theFlForce==TRUE theUGiven will be changed forcefully
1684 //            even if theUGiven is already inscibed in the boundary
1685 //            (if it is possible; i.e. if new theUGiven is inscribed
1686 //            in the boundary, too).
1687 //=======================================================================
1688 Standard_Boolean InscribePoint( const Standard_Real theUfTarget,
1689                                 const Standard_Real theUlTarget,
1690                                 Standard_Real& theUGiven,
1691                                 const Standard_Real theTol2D,
1692                                 const Standard_Real thePeriod,
1693                                 const Standard_Boolean theFlForce)
1694 {
1695   if(Precision::IsInfinite(theUGiven))
1696   {
1697     return Standard_False;
1698   }
1699
1700   if((theUfTarget - theUGiven <= theTol2D) &&
1701               (theUGiven - theUlTarget <= theTol2D))
1702   {//it has already inscribed
1703
1704     /*
1705              Utf    U      Utl
1706               +     *       +
1707     */
1708
1709     if(theFlForce)
1710     {
1711       Standard_Real anUtemp = theUGiven + thePeriod;
1712       if((theUfTarget - anUtemp <= theTol2D) &&
1713                 (anUtemp - theUlTarget <= theTol2D))
1714       {
1715         theUGiven = anUtemp;
1716         return Standard_True;
1717       }
1718
1719       anUtemp = theUGiven - thePeriod;
1720       if((theUfTarget - anUtemp <= theTol2D) &&
1721                 (anUtemp - theUlTarget <= theTol2D))
1722       {
1723         theUGiven = anUtemp;
1724       }
1725     }
1726
1727     return Standard_True;
1728   }
1729
1730   const Standard_Real aUf = theUfTarget - theTol2D;
1731   const Standard_Real aUl = aUf + thePeriod;
1732
1733   theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl);
1734
1735   return ((theUfTarget - theUGiven <= theTol2D) &&
1736           (theUGiven - theUlTarget <= theTol2D));
1737 }
1738
1739 //=======================================================================
1740 //function : InscribeInterval
1741 //purpose  : Shifts theRange in order to make at least one of its
1742 //            boundary in the range [theUfTarget, theUlTarget]
1743 //=======================================================================
1744 static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
1745                                          const Standard_Real theUlTarget,
1746                                          Bnd_Range &theRange,
1747                                          const Standard_Real theTol2D,
1748                                          const Standard_Real thePeriod)
1749 {
1750   Standard_Real anUpar = 0.0;
1751   if (!theRange.GetMin(anUpar))
1752   {
1753     return Standard_False;
1754   }
1755
1756   const Standard_Real aDelta = theRange.Delta();
1757   if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1758           theTol2D, thePeriod, (Abs(theUlTarget-anUpar) < theTol2D)))
1759   {
1760     theRange.SetVoid();
1761     theRange.Add(anUpar);
1762     theRange.Add(anUpar + aDelta);
1763   }
1764   else
1765   {
1766     if (!theRange.GetMAX(anUpar))
1767     {
1768       return Standard_False;
1769     }
1770
1771     if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1772           theTol2D, thePeriod, (Abs(theUfTarget-anUpar) < theTol2D)))
1773     {
1774       theRange.SetVoid();
1775       theRange.Add(anUpar);
1776       theRange.Add(anUpar - aDelta);
1777     }
1778     else
1779     {
1780       return Standard_False;
1781     }
1782   }
1783
1784   return Standard_True;
1785 }
1786
1787 //=======================================================================
1788 //function : ExcludeNearElements
1789 //purpose  : Checks if theArr contains two almost equal elements.
1790 //            If it is true then one of equal elements will be excluded
1791 //            (made infinite).
1792 //           Returns TRUE if at least one element of theArr has been changed.
1793 //ATTENTION!!!
1794 //           1. Every not infinite element of theArr is considered to be
1795 //            in [0, T] interval (where T is the period);
1796 //           2. theArr must be sorted in ascending order.
1797 //=======================================================================
1798 static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
1799                                             const Standard_Integer theNOfMembers,
1800                                             const Standard_Real theUSurf1f,
1801                                             const Standard_Real theUSurf1l,
1802                                             const Standard_Real theTol)
1803 {
1804   Standard_Boolean aRetVal = Standard_False;
1805   for(Standard_Integer i = 1; i < theNOfMembers; i++)
1806   {
1807     Standard_Real &anA = theArr[i],
1808                   &anB = theArr[i-1];
1809
1810     //Here, anA >= anB
1811
1812     if(Precision::IsInfinite(anA))
1813       break;
1814
1815     if((anA-anB) < theTol)
1816     {
1817       if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l))
1818       anA = (anA + anB)/2.0;
1819       else
1820         anA = anB;
1821
1822       //Make this element infinite an forget it
1823       //(we will not use it in next iterations).
1824       anB = Precision::Infinite();
1825       aRetVal = Standard_True;
1826     }
1827   }
1828
1829   return aRetVal;
1830 }
1831
1832 //=======================================================================
1833 //function : AddPointIntoWL
1834 //purpose  : Surf1 is a surface, whose U-par is variable.
1835 //           If theFlBefore == TRUE then we enable the U1-parameter
1836 //            of the added point to be less than U1-parameter of
1837 //           previously added point (in general U1-parameter is
1838 //           always increased; therefore, this situation is abnormal).
1839 //           If theOnlyCheck==TRUE then no point will be added to theLine.
1840 //=======================================================================
1841 static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
1842                                         const IntSurf_Quadric& theQuad2,
1843                                         const ComputationMethods::stCoeffsValue& theCoeffs,
1844                                         const Standard_Boolean isTheReverse,
1845                                         const Standard_Boolean isThePrecise,
1846                                         const gp_Pnt2d& thePntOnSurf1,
1847                                         const gp_Pnt2d& thePntOnSurf2,
1848                                         const Standard_Real theUfSurf1,
1849                                         const Standard_Real theUlSurf1,
1850                                         const Standard_Real theUfSurf2,
1851                                         const Standard_Real theUlSurf2,
1852                                         const Standard_Real theVfSurf1,
1853                                         const Standard_Real theVlSurf1,
1854                                         const Standard_Real theVfSurf2,
1855                                         const Standard_Real theVlSurf2,
1856                                         const Standard_Real thePeriodOfSurf1,
1857                                         const Handle(IntSurf_LineOn2S)& theLine,
1858                                         const Standard_Integer theWLIndex,
1859                                         const Standard_Real theTol3D,
1860                                         const Standard_Real theTol2D,
1861                                         const Standard_Boolean theFlBefore = Standard_False,
1862                                         const Standard_Boolean theOnlyCheck = Standard_False)
1863 {
1864   //Check if the point is in the domain or can be inscribed in the domain after adjusting.
1865
1866   gp_Pnt  aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
1867           aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
1868
1869   Standard_Real aU1par = thePntOnSurf1.X();
1870
1871   // aU1par always increases. Therefore, we must reduce its
1872   // value in order to continue creation of WLine.
1873   if(!InscribePoint(theUfSurf1, theUlSurf1, aU1par, theTol2D,
1874                   thePeriodOfSurf1, aU1par > 0.5*(theUfSurf1+theUlSurf1)))
1875     return Standard_False;
1876
1877   if ((theLine->NbPoints() > 0) &&
1878       ((theUlSurf1 - theUfSurf1) >= thePeriodOfSurf1) &&
1879       (((aU1par + thePeriodOfSurf1 - theUlSurf1) <= theTol2D) ||
1880        ((aU1par - thePeriodOfSurf1 - theUfSurf1) >= theTol2D)))
1881   {
1882     // aU1par can be adjusted to both theUlSurf1 and theUfSurf1
1883     // with equal possibilities. This code fragment allows choosing
1884     // correct parameter from these two variants.
1885
1886     Standard_Real aU1 = 0.0, aV1 = 0.0;
1887     if (isTheReverse)
1888     {
1889       theLine->Value(theLine->NbPoints()).ParametersOnS2(aU1, aV1);
1890     }
1891     else
1892     {
1893       theLine->Value(theLine->NbPoints()).ParametersOnS1(aU1, aV1);
1894     }
1895
1896     const Standard_Real aDelta = aU1 - aU1par;
1897     if (2.0*Abs(aDelta) > thePeriodOfSurf1)
1898     {
1899       aU1par += Sign(thePeriodOfSurf1, aDelta);
1900     }
1901   }
1902
1903   Standard_Real aU2par = thePntOnSurf2.X();
1904   if(!InscribePoint(theUfSurf2, theUlSurf2, aU2par, theTol2D,
1905                                     thePeriodOfSurf1, Standard_False))
1906     return Standard_False;
1907
1908   Standard_Real aV1par = thePntOnSurf1.Y();
1909   if((aV1par - theVlSurf1 > theTol2D) || (theVfSurf1 - aV1par > theTol2D))
1910     return Standard_False;
1911
1912   Standard_Real aV2par = thePntOnSurf2.Y();
1913   if((aV2par -  theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
1914     return Standard_False;
1915
1916   //Get intersection point and add it in the WL
1917   IntSurf_PntOn2S aPnt;
1918
1919   if(isTheReverse)
1920   {
1921     aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
1922                         aU2par, aV2par,
1923                         aU1par, aV1par);
1924   }
1925   else
1926   {
1927     aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
1928                         aU1par, aV1par,
1929                         aU2par, aV2par);
1930   }
1931
1932   Standard_Integer aNbPnts = theLine->NbPoints();
1933   if(aNbPnts > 0)
1934   {
1935     Standard_Real aUl = 0.0, aVl = 0.0;
1936     const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts);
1937     if(isTheReverse)
1938       aPlast.ParametersOnS2(aUl, aVl);
1939     else
1940       aPlast.ParametersOnS1(aUl, aVl);
1941
1942     if(!theFlBefore && (aU1par <= aUl))
1943     {
1944       //Parameter value must be increased if theFlBefore == FALSE.
1945
1946       aU1par += thePeriodOfSurf1;
1947
1948       //The condition is as same as in
1949       //InscribePoint(...) function
1950       if((theUfSurf1 - aU1par > theTol2D) ||
1951          (aU1par - theUlSurf1 > theTol2D))
1952       {
1953         //New aU1par is out of target interval.
1954         //Go back to old value.
1955
1956         return Standard_False;
1957       }
1958     }
1959
1960     if (theOnlyCheck)
1961       return Standard_True;
1962
1963     //theTol2D is minimal step along parameter changed.
1964     //Therefore, if we apply this minimal step two
1965     //neighbour points will be always "same". Consequently,
1966     //we should reduce tolerance for IsSame checking.
1967     const Standard_Real aDTol = 1.0-Epsilon(1.0);
1968     if(aPnt.IsSame(aPlast, theTol3D*aDTol, theTol2D*aDTol))
1969     {
1970       theLine->RemovePoint(aNbPnts);
1971     }
1972   }
1973
1974   if (theOnlyCheck)
1975     return Standard_True;
1976
1977   theLine->Add(aPnt);
1978
1979   if(!isThePrecise)
1980     return Standard_True;
1981
1982   //Try to precise existing WLine
1983   aNbPnts = theLine->NbPoints();
1984   if(aNbPnts >= 3)
1985   {
1986     Standard_Real aU1 = 0.0, aU2 = 0.0, aU3 = 0.0, aV = 0.0;
1987     if(isTheReverse)
1988     {
1989       theLine->Value(aNbPnts).ParametersOnS2(aU3, aV);
1990       theLine->Value(aNbPnts-1).ParametersOnS2(aU2, aV);
1991       theLine->Value(aNbPnts-2).ParametersOnS2(aU1, aV);
1992     }
1993     else
1994     {
1995       theLine->Value(aNbPnts).ParametersOnS1(aU3, aV);
1996       theLine->Value(aNbPnts-1).ParametersOnS1(aU2, aV);
1997       theLine->Value(aNbPnts-2).ParametersOnS1(aU1, aV);
1998     }
1999
2000     const Standard_Real aStepPrev = aU2-aU1;
2001     const Standard_Real aStep = aU3-aU2;
2002
2003     const Standard_Integer aDeltaStep = RealToInt(aStepPrev/aStep);
2004
2005     if((1 < aDeltaStep) && (aDeltaStep < 2000))
2006     {
2007       //Add new points in case of non-uniform distribution of existing points
2008       SeekAdditionalPoints( theQuad1, theQuad2, theLine,
2009                             theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
2010                             aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse);
2011     }
2012   }
2013
2014   return Standard_True;
2015 }
2016
2017 //=======================================================================
2018 //function : AddBoundaryPoint
2019 //purpose  : Find intersection point on V-boundary.
2020 //=======================================================================
2021 void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
2022                                           const Standard_Real theU1,
2023                                           const Standard_Real theU1Min,
2024                                           const Standard_Real theU2,
2025                                           const Standard_Real theV1,
2026                                           const Standard_Real theV1Prev,
2027                                           const Standard_Real theV2,
2028                                           const Standard_Real theV2Prev,
2029                                           const Standard_Integer theWLIndex,
2030                                           const Standard_Boolean theFlForce,
2031                                           Standard_Boolean& isTheFound1,
2032                                           Standard_Boolean& isTheFound2) const
2033 {
2034   Standard_Real aUSurf1f = 0.0, //const
2035                 aUSurf1l = 0.0,
2036                 aVSurf1f = 0.0,
2037                 aVSurf1l = 0.0;
2038   Standard_Real aUSurf2f = 0.0, //const
2039                 aUSurf2l = 0.0,
2040                 aVSurf2f = 0.0,
2041                 aVSurf2l = 0.0;
2042
2043   myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
2044   myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
2045
2046   const Standard_Integer aSize = 4;
2047   const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l};
2048
2049   StPInfo aUVPoint[aSize];
2050
2051   for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2)
2052   {
2053     const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2,
2054                         aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev;
2055
2056     const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2;
2057
2058     for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++)
2059     {
2060       const Standard_Integer anIndex = anIDSurf+anIDBound;
2061
2062       aUVPoint[anIndex].mySurfID = anIDSurf;
2063
2064       if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) &&
2065           ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0))
2066       {
2067         continue;
2068       }
2069
2070       //Segment [aVf, aVl] intersects at least one V-boundary (first or last)
2071       // (in general, case is possible, when aVf > aVl).
2072
2073       // Precise intersection point
2074       const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex],
2075                                                     (anIDSurf == 0) ? theV2 : theV1,
2076                                                     theU2, theU1,
2077                                                     aUVPoint[anIndex].myU1);
2078
2079       // aUVPoint[anIndex].myU1 is considered to be nearer to theU1 than
2080       // to theU1+/-Period
2081       if (!aRes || (aUVPoint[anIndex].myU1 >= theU1) ||
2082                               (aUVPoint[anIndex].myU1 < theU1Min))
2083       {
2084         //Intersection point is not found or out of the domain
2085         aUVPoint[anIndex].myU1 = RealLast();
2086         continue;
2087       }
2088       else
2089       {
2090         //intersection point is found
2091
2092         Standard_Real &aU1 = aUVPoint[anIndex].myU1,
2093                       &aU2 = aUVPoint[anIndex].myU2,
2094                       &aV1 = aUVPoint[anIndex].myV1,
2095                       &aV2 = aUVPoint[anIndex].myV2;
2096         aU2 = theU2;
2097         aV1 = theV1;
2098         aV2 = theV2;
2099
2100         if(!ComputationMethods::
2101                   CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2))
2102         {
2103           // Found point is wrong
2104           aU1 = RealLast();
2105           continue;
2106         }
2107
2108         //Point on true V-boundary.
2109         if(aTS == SearchV1)
2110           aV1 = anArrVzad[anIndex];
2111         else //if(aTS[anIndex] == SearchV2)
2112           aV2 = anArrVzad[anIndex];
2113       }
2114     }
2115   }
2116
2117   // Sort with acceding U1-parameter.
2118   std::sort(aUVPoint, aUVPoint+aSize);
2119
2120   isTheFound1 = isTheFound2 = Standard_False;
2121
2122   //Adding found points on boundary in the WLine.
2123   for(Standard_Integer i = 0; i < aSize; i++)
2124   {
2125     if(aUVPoint[i].myU1 == RealLast())
2126       break;
2127
2128     if(!AddPointIntoWL(myQuad1, myQuad2, myCoeffs, myIsReverse, Standard_False,
2129                         gp_Pnt2d(aUVPoint[i].myU1, aUVPoint[i].myV1),
2130                         gp_Pnt2d(aUVPoint[i].myU2, aUVPoint[i].myV2),
2131                         aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2132                         aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, myPeriod,
2133                         theWL->Curve(), theWLIndex, myTol3D, myTol2D, theFlForce))
2134     {
2135       continue;
2136     }
2137
2138     if(aUVPoint[i].mySurfID == 0)
2139     {
2140       isTheFound1 = Standard_True;
2141     }
2142     else
2143     {
2144       isTheFound2 = Standard_True;
2145     }
2146   }
2147 }
2148
2149 //=======================================================================
2150 //function : SeekAdditionalPoints
2151 //purpose  : Inserts additional intersection points between neighbor points.
2152 //            This action is bone with several iterations. During every iteration,
2153 //          new point is inserted in middle of every interval.
2154 //            The process will be finished if NbPoints >= theMinNbPoints.
2155 //=======================================================================
2156 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
2157                                   const IntSurf_Quadric& theQuad2,
2158                                   const Handle(IntSurf_LineOn2S)& theLine,
2159                                   const ComputationMethods::stCoeffsValue& theCoeffs,
2160                                   const Standard_Integer theWLIndex,
2161                                   const Standard_Integer theMinNbPoints,
2162                                   const Standard_Integer theStartPointOnLine,
2163                                   const Standard_Integer theEndPointOnLine,
2164                                   const Standard_Real theTol2D,
2165                                   const Standard_Real thePeriodOfSurf2,
2166                                   const Standard_Boolean isTheReverse)
2167 {
2168   if(theLine.IsNull())
2169     return;
2170
2171   Standard_Integer aNbPoints = theEndPointOnLine - theStartPointOnLine + 1;
2172   if(aNbPoints >= theMinNbPoints)
2173   {
2174     return;
2175   }
2176
2177   Standard_Real aMinDeltaParam = theTol2D;
2178
2179   {
2180     Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
2181
2182     if(isTheReverse)
2183     {
2184       theLine->Value(theStartPointOnLine).ParametersOnS2(u1, v1);
2185       theLine->Value(theEndPointOnLine).ParametersOnS2(u2, v2);
2186     }
2187     else
2188     {
2189       theLine->Value(theStartPointOnLine).ParametersOnS1(u1, v1);
2190       theLine->Value(theEndPointOnLine).ParametersOnS1(u2, v2);
2191     }
2192
2193     aMinDeltaParam = Max(Abs(u2 - u1)/IntToReal(theMinNbPoints), aMinDeltaParam);
2194   }
2195
2196   Standard_Integer aLastPointIndex = theEndPointOnLine;
2197   Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
2198
2199   Standard_Integer aNbPointsPrev = 0;
2200   while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
2201   {
2202     aNbPointsPrev = aNbPoints;
2203     for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
2204     {
2205       Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
2206       Standard_Real U1l = 0.0, V1l = 0.0; //last  point in 1st suraface
2207
2208       Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface
2209       Standard_Real U2l = 0.0, V2l = 0.0; //last  point in 2nd suraface
2210
2211       lp = fp+1;
2212
2213       if(isTheReverse)
2214       {
2215         theLine->Value(fp).ParametersOnS2(U1f, V1f);
2216         theLine->Value(lp).ParametersOnS2(U1l, V1l);
2217
2218         theLine->Value(fp).ParametersOnS1(U2f, V2f);
2219         theLine->Value(lp).ParametersOnS1(U2l, V2l);
2220       }
2221       else
2222       {
2223         theLine->Value(fp).ParametersOnS1(U1f, V1f);
2224         theLine->Value(lp).ParametersOnS1(U1l, V1l);
2225
2226         theLine->Value(fp).ParametersOnS2(U2f, V2f);
2227         theLine->Value(lp).ParametersOnS2(U2l, V2l);
2228       }
2229
2230       if(Abs(U1l - U1f) <= aMinDeltaParam)
2231       {
2232         //Step is minimal. It is not necessary to divide it.
2233         continue;
2234       }
2235
2236       U1prec = 0.5*(U1f+U1l);
2237
2238       if(!ComputationMethods::
2239             CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
2240       {
2241         continue;
2242       }
2243
2244       MinMax(U2f, U2l);
2245       if(!InscribePoint(U2f, U2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False))
2246       {
2247         continue;
2248       }
2249
2250       const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec));
2251       const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
2252
2253 #ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
2254       std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl;
2255 #endif
2256
2257       IntSurf_PntOn2S anIP;
2258       if(isTheReverse)
2259       {
2260         anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
2261       }
2262       else
2263       {
2264         anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
2265       }
2266
2267       theLine->InsertBefore(lp, anIP);
2268
2269       aNbPoints++;
2270       aLastPointIndex++;
2271     }
2272
2273     if(aNbPoints >= theMinNbPoints)
2274     {
2275       return;
2276     }
2277   }
2278 }
2279
2280 //=======================================================================
2281 //function : BoundariesComputing
2282 //purpose  : Computes true domain of future intersection curve (allows
2283 //            avoiding excess iterations)
2284 //=======================================================================
2285 Standard_Boolean WorkWithBoundaries::
2286             BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
2287                                 const Standard_Real thePeriod,
2288                                 Bnd_Range theURange[])
2289 {
2290   //All comments to this method is related to the comment
2291   //to ComputationMethods class
2292
2293   //So, we have the equation
2294   //    cos(U2-FI2)=B*cos(U1-FI1)+C
2295   //Evidently,
2296   //    -1 <= B*cos(U1-FI1)+C <= 1
2297
2298   if (theCoeffs.mB > 0.0)
2299   {
2300     // -(1+C)/B <= cos(U1-FI1) <= (1-C)/B
2301
2302     if (theCoeffs.mB + Abs(theCoeffs.mC) < -1.0)
2303     {
2304       //(1-C)/B < -1 or -(1+C)/B > 1  ==> No solution
2305
2306       return Standard_False;
2307     }
2308     else if (theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
2309     {
2310       //(1-C)/B >= 1 and -(1+C)/B <= -1 ==> U=[0;2*PI]+aFI1
2311       theURange[0].Add(theCoeffs.mFI1);
2312       theURange[0].Add(thePeriod + theCoeffs.mFI1);
2313     }
2314     else if ((1 + theCoeffs.mC <= theCoeffs.mB) &&
2315              (theCoeffs.mB <= 1 - theCoeffs.mC))
2316     {
2317       //(1-C)/B >= 1 and -(1+C)/B >= -1 ==>
2318       //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
2319       //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB)
2320
2321       Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
2322       if(anArg > 1.0)
2323         anArg = 1.0;
2324       if(anArg < -1.0)
2325         anArg = -1.0;
2326
2327       const Standard_Real aDAngle = acos(anArg);
2328       theURange[0].Add(theCoeffs.mFI1);
2329       theURange[0].Add(aDAngle + theCoeffs.mFI1);
2330       theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
2331       theURange[1].Add(thePeriod + theCoeffs.mFI1);
2332     }
2333     else if ((1 - theCoeffs.mC <= theCoeffs.mB) &&
2334              (theCoeffs.mB <= 1 + theCoeffs.mC))
2335     {
2336       //(1-C)/B <= 1 and -(1+C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1
2337       //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)
2338
2339       Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
2340       if(anArg > 1.0)
2341         anArg = 1.0;
2342       if(anArg < -1.0)
2343         anArg = -1.0;
2344
2345       const Standard_Real aDAngle = acos(anArg);
2346       theURange[0].Add(aDAngle + theCoeffs.mFI1);
2347       theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
2348     }
2349     else if (theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
2350     {
2351       //(1-C)/B <= 1 and -(1+C)/B >= -1 ==>
2352       //(U=[aDAngle1;aDAngle2]+aFI1) ||
2353       //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2354       //where aDAngle1 = acos((1 - myCoeffs.mC) / myCoeffs.mB),
2355       //      aDAngle2 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).
2356
2357       Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB,
2358                     anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB;
2359       if(anArg1 > 1.0)
2360         anArg1 = 1.0;
2361       if(anArg1 < -1.0)
2362         anArg1 = -1.0;
2363
2364       if(anArg2 > 1.0)
2365         anArg2 = 1.0;
2366       if(anArg2 < -1.0)
2367         anArg2 = -1.0;
2368
2369       const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2370       //(U=[aDAngle1;aDAngle2]+aFI1) ||
2371       //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
2372       theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
2373       theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
2374       theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
2375       theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
2376     }
2377     else
2378     {
2379       return Standard_False;
2380     }
2381   }
2382   else if (theCoeffs.mB < 0.0)
2383   {
2384     // (1-C)/B <= cos(U1-FI1) <= -(1+C)/B
2385
2386     if (theCoeffs.mB + Abs(theCoeffs.mC) > 1.0)
2387     {
2388       // -(1+C)/B < -1 or (1-C)/B > 1 ==> No solutions
2389       return Standard_False;
2390     }
2391     else if (-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
2392     {
2393       //  -(1+C)/B >= 1 and (1-C)/B <= -1 ==> U=[0;2*PI]+aFI1
2394       theURange[0].Add(theCoeffs.mFI1);
2395       theURange[0].Add(thePeriod + theCoeffs.mFI1);
2396     }
2397     else if ((-theCoeffs.mC - 1 <= theCoeffs.mB) &&
2398              (theCoeffs.mB <= theCoeffs.mC - 1))
2399     {
2400       //  -(1+C)/B >= 1 and (1-C)/B >= -1 ==>
2401       //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
2402       //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)
2403
2404       Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
2405       if(anArg > 1.0)
2406         anArg = 1.0;
2407       if(anArg < -1.0)
2408         anArg = -1.0;
2409
2410       const Standard_Real aDAngle = acos(anArg);
2411       theURange[0].Add(theCoeffs.mFI1);
2412       theURange[0].Add(aDAngle + theCoeffs.mFI1);
2413       theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
2414       theURange[1].Add(thePeriod + theCoeffs.mFI1);
2415     }
2416     else if ((theCoeffs.mC - 1 <= theCoeffs.mB) &&
2417              (theCoeffs.mB <= -theCoeffs.mB - 1))
2418     {
2419       //  -(1+C)/B <= 1 and (1-C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1,
2420       //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).
2421
2422       Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
2423       if(anArg > 1.0)
2424         anArg = 1.0;
2425       if(anArg < -1.0)
2426         anArg = -1.0;
2427
2428       const Standard_Real aDAngle = acos(anArg);
2429       theURange[0].Add(aDAngle + theCoeffs.mFI1);
2430       theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
2431     }
2432     else if (-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
2433     {
2434       //  -(1+C)/B <= 1 and (1-C)/B >= -1 ==>
2435       //(U=[aDAngle1;aDAngle2]+aFI1) || (U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1),
2436       //where aDAngle1 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB),
2437       //      aDAngle2 = acos((1 - myCoeffs.mC) / myCoeffs.mB)
2438
2439       Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB,
2440                     anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB;
2441       if(anArg1 > 1.0)
2442         anArg1 = 1.0;
2443       if(anArg1 < -1.0)
2444         anArg1 = -1.0;
2445
2446       if(anArg2 > 1.0)
2447         anArg2 = 1.0;
2448       if(anArg2 < -1.0)
2449         anArg2 = -1.0;
2450
2451       const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
2452       theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
2453       theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
2454       theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
2455       theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
2456     }
2457     else
2458     {
2459       return Standard_False;
2460     }
2461   }
2462   else
2463   {
2464     return Standard_False;
2465   }
2466
2467   return Standard_True;
2468 }
2469
2470 //=======================================================================
2471 //function : CriticalPointsComputing
2472 //purpose  : theNbCritPointsMax contains true number of critical points.
2473 //            It must be initialized correctly before function calling
2474 //=======================================================================
2475 static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs,
2476                                     const Standard_Real theUSurf1f,
2477                                     const Standard_Real theUSurf1l,
2478                                     const Standard_Real theUSurf2f,
2479                                     const Standard_Real theUSurf2l,
2480                                     const Standard_Real thePeriod,
2481                                     const Standard_Real theTol2D,
2482                                     Standard_Integer& theNbCritPointsMax,
2483                                     Standard_Real theU1crit[])
2484 {
2485   //[0...1] - in these points parameter U1 goes through
2486   //          the seam-edge of the first cylinder.
2487   //[2...3] - First and last U1 parameter.
2488   //[4...5] - in these points parameter U2 goes through
2489   //          the seam-edge of the second cylinder.
2490   //[6...9] - in these points an intersection line goes through
2491   //          U-boundaries of the second surface.
2492   //[10...11] - Boundary of monotonicity interval of U2(U1) function
2493   //            (see CylCylMonotonicity() function)
2494
2495   theU1crit[0] = 0.0;
2496   theU1crit[1] = thePeriod;
2497   theU1crit[2] = theUSurf1f;
2498   theU1crit[3] = theUSurf1l;
2499
2500   const Standard_Real aCOS = cos(theCoeffs.mFI2);
2501   const Standard_Real aBSB = Abs(theCoeffs.mB);
2502   if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
2503   {
2504     Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
2505     if(anArg > 1.0)
2506       anArg = 1.0;
2507     if(anArg < -1.0)
2508       anArg = -1.0;
2509
2510     theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
2511     theU1crit[5] =  acos(anArg) + theCoeffs.mFI1;
2512   }
2513
2514   Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
2515   Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
2516   MinMax(aSf, aSl);
2517
2518   //In accorance with pure mathematic, theU1crit[6] and [8]
2519   //must be -Precision::Infinite() instead of used +Precision::Infinite()
2520   theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2521                   -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2522                   Precision::Infinite();
2523   theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2524                   -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2525                   Precision::Infinite();
2526   theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2527                   acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2528                   Precision::Infinite();
2529   theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
2530                   acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
2531                   Precision::Infinite();
2532
2533   theU1crit[10] = theCoeffs.mFI1;
2534   theU1crit[11] = M_PI+theCoeffs.mFI1;
2535
2536   //preparative treatment of array. This array must have faled to contain negative
2537   //infinity number
2538
2539   for(Standard_Integer i = 0; i < theNbCritPointsMax; i++)
2540   {
2541     if(Precision::IsInfinite(theU1crit[i]))
2542     {
2543       continue;
2544     }
2545
2546     theU1crit[i] = fmod(theU1crit[i], thePeriod);
2547     if(theU1crit[i] < 0.0)
2548       theU1crit[i] += thePeriod;
2549   }
2550
2551   //Here all not infinite elements of theU1crit are in [0, thePeriod) range
2552
2553   do
2554   {
2555     std::sort(theU1crit, theU1crit + theNbCritPointsMax);
2556   }
2557   while(ExcludeNearElements(theU1crit, theNbCritPointsMax,
2558                             theUSurf1f, theUSurf1l, theTol2D));
2559
2560   //Here all not infinite elements in theU1crit are different and sorted.
2561
2562   while(theNbCritPointsMax > 0)
2563   {
2564     Standard_Real &anB = theU1crit[theNbCritPointsMax-1];
2565     if(Precision::IsInfinite(anB))
2566     {
2567       theNbCritPointsMax--;
2568       continue;
2569     }
2570
2571     //1st not infinte element is found
2572
2573     if(theNbCritPointsMax == 1)
2574       break;
2575
2576     //Here theNbCritPointsMax > 1
2577
2578     Standard_Real &anA = theU1crit[0];
2579
2580     //Compare 1st and last significant elements of theU1crit
2581     //They may still differs by period.
2582
2583     if (Abs(anB - anA - thePeriod) < theTol2D)
2584     {//E.g. anA == 2.0e-17, anB == (thePeriod-1.0e-18)
2585       anA = (anA + anB - thePeriod)/2.0;
2586       anB = Precision::Infinite();
2587       theNbCritPointsMax--;
2588     }
2589
2590     //Out of "while(theNbCritPointsMax > 0)" cycle.
2591     break;
2592   }
2593
2594   //Attention! Here theU1crit may be unsorted.
2595 }
2596
2597 //=======================================================================
2598 //function : BoundaryEstimation
2599 //purpose  : Rough estimation of the parameter range.
2600 //=======================================================================
2601 void WorkWithBoundaries::BoundaryEstimation(const gp_Cylinder& theCy1,
2602                                             const gp_Cylinder& theCy2,
2603                                             Bnd_Range& theOutBoxS1,
2604                                             Bnd_Range& theOutBoxS2) const
2605 {
2606   const gp_Dir &aD1 = theCy1.Axis().Direction(),
2607                &aD2 = theCy2.Axis().Direction();
2608   const Standard_Real aR1 = theCy1.Radius(),
2609                       aR2 = theCy2.Radius();
2610
2611   //Let consider a parallelogram. Its edges are parallel to aD1 and aD2.
2612   //Its altitudes are equal to 2*aR1 and 2*aR2 (diameters of the cylinders).
2613   //In fact, this parallelogram is a projection of the cylinders to the plane
2614   //created by the intersected axes aD1 and aD2 (if the axes are skewed then
2615   //one axis can be translated by parallel shifting till intersection).
2616
2617   const Standard_Real aCosA = aD1.Dot(aD2);
2618   const Standard_Real aSqSinA = aD1.XYZ().CrossSquareMagnitude(aD2.XYZ());
2619
2620   //If sine is small then it can be compared with angle.
2621   if (aSqSinA < Precision::Angular()*Precision::Angular())
2622     return;
2623
2624   //Half of delta V. Delta V is a distance between
2625   //projections of two opposite parallelogram vertices
2626   //(joined by the maximal diagonal) to the cylinder axis.
2627   const Standard_Real aSinA = sqrt(aSqSinA);
2628   const Standard_Real anAbsCosA = Abs(aCosA);
2629   const Standard_Real aHDV1 = (aR1 * anAbsCosA + aR2) / aSinA,
2630                       aHDV2 = (aR2 * anAbsCosA + aR1) / aSinA;
2631
2632 #ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
2633   //The code in this block is created for test only.It is stupidly to create
2634   //OCCT-test for the method, which will be changed possibly never.
2635   std::cout << "Reference values: aHDV1 = " << aHDV1 << "; aHDV2 = " << aHDV2 << std::endl;
2636 #endif
2637
2638   //V-parameters of intersection point of the axes (in case of skewed axes,
2639   //see comment above).
2640   Standard_Real aV01 = 0.0, aV02 = 0.0;
2641   ExtremaLineLine(theCy1.Axis(), theCy2.Axis(), aCosA, aSqSinA, aV01, aV02);
2642
2643   theOutBoxS1.Add(aV01 - aHDV1);
2644   theOutBoxS1.Add(aV01 + aHDV1);
2645
2646   theOutBoxS2.Add(aV02 - aHDV2);
2647   theOutBoxS2.Add(aV02 + aHDV2);
2648
2649   theOutBoxS1.Enlarge(Precision::Confusion());
2650   theOutBoxS2.Enlarge(Precision::Confusion());
2651
2652   Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
2653   myUVSurf1.Get(aU1, aV1, aU2, aV2);
2654   theOutBoxS1.Common(Bnd_Range(aV1, aV2));
2655
2656   myUVSurf2.Get(aU1, aV1, aU2, aV2);
2657   theOutBoxS2.Common(Bnd_Range(aV1, aV2));
2658 }
2659
2660 //=======================================================================
2661 //function : CyCyNoGeometric
2662 //purpose  :
2663 //=======================================================================
2664 static IntPatch_ImpImpIntersection::IntStatus
2665                     CyCyNoGeometric(const gp_Cylinder &theCyl1,
2666                                     const gp_Cylinder &theCyl2,
2667                                     const WorkWithBoundaries &theBW,
2668                                     Bnd_Range theRange[],
2669                                     const Standard_Integer theNbOfRanges /*=2*/,
2670                                     Standard_Boolean& isTheEmpty,
2671                                     IntPatch_SequenceOfLine& theSlin,
2672                                     IntPatch_SequenceOfPoint& theSPnt)
2673 {
2674   Standard_Real aUSurf1f = 0.0, aUSurf1l = 0.0,
2675                 aUSurf2f = 0.0, aUSurf2l = 0.0,
2676                 aVSurf1f = 0.0, aVSurf1l = 0.0,
2677                 aVSurf2f = 0.0, aVSurf2l = 0.0;
2678
2679   theBW.UVS1().Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
2680   theBW.UVS2().Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
2681
2682   Bnd_Range aRangeS1, aRangeS2;
2683   theBW.BoundaryEstimation(theCyl1, theCyl2, aRangeS1, aRangeS2);
2684   if (aRangeS1.IsVoid() || aRangeS2.IsVoid())
2685     return IntPatch_ImpImpIntersection::IntStatus_OK;
2686
2687   {
2688     //Quotation of the message from issue #26894 (author MSV):
2689     //"We should return fail status from intersector if the result should be an
2690     //infinite curve of non-analytical type... I propose to define the limit for the
2691     //extent as the radius divided by 1e+2 and multiplied by 1e+7.
2692     //Thus, taking into account the number of valuable digits (15), we provide reliable
2693     //computations with an error not exceeding R/100."
2694     const Standard_Real aF = 1.0e+5;
2695     const Standard_Real aMaxV1Range = aF*theCyl1.Radius(), aMaxV2Range = aF*theCyl2.Radius();
2696     if ((aRangeS1.Delta() > aMaxV1Range) || (aRangeS2.Delta() > aMaxV2Range))
2697       return IntPatch_ImpImpIntersection::IntStatus_InfiniteSectionCurve;
2698   }
2699
2700   const ComputationMethods::stCoeffsValue &anEquationCoeffs = theBW.SICoeffs();
2701   const IntSurf_Quadric& aQuad1 = theBW.GetQSurface(1);
2702   const IntSurf_Quadric& aQuad2 = theBW.GetQSurface(2);
2703   const Standard_Boolean isReversed = theBW.IsReversed();
2704   const Standard_Real aTol2D = theBW.Get2dTolerance();
2705   const Standard_Real aTol3D = theBW.Get3dTolerance();
2706   const Standard_Real aPeriod = 2.0*M_PI;
2707   const Standard_Integer aNbMaxPoints = 2000;
2708   const Standard_Integer aNbMinPoints = 200;
2709   const Standard_Integer aNbPoints = Min(Max(aNbMinPoints,
2710                                   RealToInt(20.0*theCyl1.Radius())), aNbMaxPoints);
2711   const Standard_Real aStepMin = aTol2D,
2712                                   aStepMax = (aUSurf1l - aUSurf1f > M_PI / 100.0) ?
2713                                   (aUSurf1l - aUSurf1f) / IntToReal(aNbPoints) :
2714                                   aUSurf1l - aUSurf1f;
2715
2716   //The main idea of the algorithm is to change U1-parameter
2717   //(U-parameter of theCyl1) from aU1f to aU1l with some step
2718   //(step is adaptive) and to obtain set of intersection points.
2719
2720   for (Standard_Integer i = 0; i < theNbOfRanges; i++)
2721   {
2722     if (theRange[i].IsVoid())
2723       continue;
2724
2725     InscribeInterval(aUSurf1f, aUSurf1l, theRange[i], aTol2D, aPeriod);
2726   }
2727
2728   if (theRange[0].Union(theRange[1]))
2729   {
2730     // Works only if (theNbOfRanges == 2).
2731     theRange[1].SetVoid();
2732   }
2733
2734   //Critical points are the value of U1-parameter in the points
2735   //where WL must be decomposed.
2736
2737   //When U1 goes through critical points its value is set up to this
2738   //parameter forcefully and the intersection point is added in the line.
2739   //After that, the WL is broken (next U1 value will be correspond to the new WL).
2740
2741   //See CriticalPointsComputing(...) function to get detail information about this array.
2742   const Standard_Integer aNbCritPointsMax = 12;
2743   Standard_Real anU1crit[aNbCritPointsMax] = { Precision::Infinite(),
2744                                                Precision::Infinite(),
2745                                                Precision::Infinite(),
2746                                                Precision::Infinite(),
2747                                                Precision::Infinite(),
2748                                                Precision::Infinite(),
2749                                                Precision::Infinite(),
2750                                                Precision::Infinite(),
2751                                                Precision::Infinite(),
2752                                                Precision::Infinite(),
2753                                                Precision::Infinite(),
2754                                                Precision::Infinite() };
2755
2756   //This list of critical points is not full because it does not contain any points
2757   //which intersection line goes through V-bounds of cylinders in.
2758   //They are computed by numerical methods on - line (during algorithm working).
2759   //The moment is caught, when intersection line goes through V-bounds of any cylinder.
2760
2761   Standard_Integer aNbCritPoints = aNbCritPointsMax;
2762   CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2763                           aPeriod, aTol2D, aNbCritPoints, anU1crit);
2764
2765   //Getting Walking-line
2766
2767   enum WLFStatus
2768   {
2769     // No points have been added in WL
2770     WLFStatus_Absent = 0,
2771     // WL contains at least one point
2772     WLFStatus_Exist = 1,
2773     // WL has been finished in some critical point
2774     // We should start new line
2775     WLFStatus_Broken = 2
2776   };
2777
2778   const Standard_Integer aNbWLines = 2;
2779   for (Standard_Integer aCurInterval = 0; aCurInterval < theNbOfRanges; aCurInterval++)
2780   {
2781     //Process every continuous region
2782     Standard_Boolean isAddedIntoWL[aNbWLines];
2783     for (Standard_Integer i = 0; i < aNbWLines; i++)
2784       isAddedIntoWL[i] = Standard_False;
2785
2786     Standard_Real anUf = 1.0, anUl = 0.0;
2787     if (!theRange[aCurInterval].GetBounds(anUf, anUl))
2788       continue;
2789
2790     const Standard_Boolean isDeltaPeriod = IsEqual(anUl - anUf, aPeriod);
2791
2792     //Inscribe and sort critical points
2793     for (Standard_Integer i = 0; i < aNbCritPoints; i++)
2794     {
2795       InscribePoint(anUf, anUl, anU1crit[i], 0.0, aPeriod, Standard_False);
2796     }
2797
2798     std::sort(anU1crit, anU1crit + aNbCritPoints);
2799
2800     while (anUf < anUl)
2801     {
2802       //Change value of U-parameter on the 1st surface from anUf to anUl
2803       //(anUf will be modified in the cycle body).
2804       //Step is computed adaptively (see comments below).
2805
2806       Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
2807       WLFStatus aWLFindStatus[aNbWLines];
2808       Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
2809       Standard_Real anUexpect[aNbWLines];
2810       Standard_Boolean isAddingWLEnabled[aNbWLines];
2811
2812       Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
2813       Handle(IntPatch_WLine) aWLine[aNbWLines];
2814       for (Standard_Integer i = 0; i < aNbWLines; i++)
2815       {
2816         aL2S[i] = new IntSurf_LineOn2S();
2817         aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
2818         aWLFindStatus[i] = WLFStatus_Absent;
2819         isAddingWLEnabled[i] = Standard_True;
2820         aU2[i] = aV1[i] = aV2[i] = 0.0;
2821         aV1Prev[i] = aV2Prev[i] = 0.0;
2822         anUexpect[i] = anUf;
2823       }
2824
2825       Standard_Real aCriticalDelta[aNbCritPointsMax] = { 0 };
2826       for (Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++)
2827       { //We are not interested in elements of aCriticalDelta array
2828         //if their index is greater than or equal to aNbCritPoints
2829
2830         aCriticalDelta[aCritPID] = anUf - anU1crit[aCritPID];
2831       }
2832
2833       Standard_Real anU1 = anUf, aMinCriticalParam = anUf;
2834       Standard_Boolean isFirst = Standard_True;
2835
2836       while (anU1 <= anUl)
2837       {
2838         //Change value of U-parameter on the 1st surface from anUf to anUl
2839         //(anUf will be modified in the cycle body). However, this cycle
2840         //can be broken if WL goes though some critical point.
2841         //Step is computed adaptively (see comments below).
2842
2843         for (Standard_Integer i = 0; i < aNbCritPoints; i++)
2844         {
2845           if ((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
2846           {
2847             //WL has gone through i-th critical point
2848             anU1 = anU1crit[i];
2849
2850             for (Standard_Integer j = 0; j < aNbWLines; j++)
2851             {
2852               aWLFindStatus[j] = WLFStatus_Broken;
2853               anUexpect[j] = anU1;
2854             }
2855
2856             break;
2857           }
2858         }
2859
2860         if (IsEqual(anU1, anUl))
2861         {
2862           for (Standard_Integer i = 0; i < aNbWLines; i++)
2863           {
2864             aWLFindStatus[i] = WLFStatus_Broken;
2865             anUexpect[i] = anU1;
2866
2867             if (isDeltaPeriod)
2868             {
2869               //if isAddedIntoWL[i] == TRUE WLine contains only one point
2870               //(which was end point of previous WLine). If we will
2871               //add point found on the current step WLine will contain only
2872               //two points. At that both these points will be equal to the
2873               //points found earlier. Therefore, new WLine will repeat
2874               //already existing WLine. Consequently, it is necessary
2875               //to forbid building new line in this case.
2876
2877               isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
2878             }
2879             else
2880             {
2881               isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
2882                                       (aWLFindStatus[i] == WLFStatus_Absent));
2883             }
2884           }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2885         }
2886         else
2887         {
2888           for (Standard_Integer i = 0; i < aNbWLines; i++)
2889           {
2890             isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
2891                                     (aWLFindStatus[i] == WLFStatus_Absent));
2892           }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2893         }
2894
2895         for (Standard_Integer i = 0; i < aNbWLines; i++)
2896         {
2897           const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
2898             aWLine[i]->Curve()->NbPoints();
2899
2900           if ((aWLFindStatus[i] == WLFStatus_Broken) ||
2901             (aWLFindStatus[i] == WLFStatus_Absent))
2902           {//Begin and end of WLine must be on boundary point
2903            //or on seam-edge strictly (if it is possible).
2904
2905             Standard_Real aTol = aTol2D;
2906             ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs,
2907                                                         aU2[i], &aTol);
2908             InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
2909
2910             aTol = Max(aTol, aTol2D);
2911
2912             if (Abs(aU2[i]) <= aTol)
2913               aU2[i] = 0.0;
2914             else if (Abs(aU2[i] - aPeriod) <= aTol)
2915               aU2[i] = aPeriod;
2916             else if (Abs(aU2[i] - aUSurf2f) <= aTol)
2917               aU2[i] = aUSurf2f;
2918             else if (Abs(aU2[i] - aUSurf2l) <= aTol)
2919               aU2[i] = aUSurf2l;
2920           }
2921           else
2922           {
2923             ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
2924             InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
2925           }
2926
2927           if (aNbPntsWL == 0)
2928           {//the line has not contained any points yet
2929             if (((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*aTol2D) &&
2930                 ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
2931                   (Abs(aU2[i] - aUSurf2l) < aTol2D)))
2932             {
2933               //In this case aU2[i] can have two values: current aU2[i] or
2934               //aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
2935               //correct value.
2936
2937               Standard_Boolean isIncreasing = Standard_True;
2938               ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs,
2939                                                       aPeriod, isIncreasing);
2940
2941               //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
2942               //then after the next step (when U1 will be increased) U2 will be
2943               //increased too. And we will go out of surface boundary.
2944               //Therefore, If U2(U1) is increasing then U2 must be equal aUSurf2f.
2945               //Analogically, if U2(U1) is decreasing.
2946               if (isIncreasing)
2947               {
2948                 aU2[i] = aUSurf2f;
2949               }
2950               else
2951               {
2952                 aU2[i] = aUSurf2l;
2953               }
2954             }
2955           }
2956           else
2957           {//aNbPntsWL > 0
2958             if (((aUSurf2l - aUSurf2f) >= aPeriod) &&
2959                 ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
2960                   (Abs(aU2[i] - aUSurf2l) < aTol2D)))
2961             {//end of the line
2962               Standard_Real aU2prev = 0.0, aV2prev = 0.0;
2963               if (isReversed)
2964                 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
2965               else
2966                 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
2967
2968               if (2.0*Abs(aU2prev - aU2[i]) > aPeriod)
2969               {
2970                 if (aU2prev > aU2[i])
2971                   aU2[i] += aPeriod;
2972                 else
2973                   aU2[i] -= aPeriod;
2974               }
2975             }
2976           }
2977
2978           ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
2979                                                               aV1[i], aV2[i]);
2980
2981           if (isFirst)
2982           {
2983             aV1Prev[i] = aV1[i];
2984             aV2Prev[i] = aV2[i];
2985           }
2986         }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2987
2988         isFirst = Standard_False;
2989
2990         //Looking for points into WLine
2991         Standard_Boolean isBroken = Standard_False;
2992         for (Standard_Integer i = 0; i < aNbWLines; i++)
2993         {
2994           if (!isAddingWLEnabled[i])
2995           {
2996             Standard_Boolean isBoundIntersect = Standard_False;
2997             if ((Abs(aV1[i] - aVSurf1f) <= aTol2D) ||
2998                 ((aV1[i] - aVSurf1f)*(aV1Prev[i] - aVSurf1f) < 0.0))
2999             {
3000               isBoundIntersect = Standard_True;
3001             }
3002             else if ((Abs(aV1[i] - aVSurf1l) <= aTol2D) ||
3003                     ((aV1[i] - aVSurf1l)*(aV1Prev[i] - aVSurf1l) < 0.0))
3004             {
3005               isBoundIntersect = Standard_True;
3006             }
3007             else if ((Abs(aV2[i] - aVSurf2f) <= aTol2D) ||
3008                     ((aV2[i] - aVSurf2f)*(aV2Prev[i] - aVSurf2f) < 0.0))
3009             {
3010               isBoundIntersect = Standard_True;
3011             }
3012             else if ((Abs(aV2[i] - aVSurf2l) <= aTol2D) ||
3013                     ((aV2[i] - aVSurf2l)*(aV2Prev[i] - aVSurf2l) < 0.0))
3014             {
3015               isBoundIntersect = Standard_True;
3016             }
3017
3018             if (aWLFindStatus[i] == WLFStatus_Broken)
3019               isBroken = Standard_True;
3020
3021             if (!isBoundIntersect)
3022             {
3023               continue;
3024             }
3025             else
3026             {
3027               anUexpect[i] = anU1;
3028             }
3029           }
3030
3031           // True if the current point already in the domain
3032           const Standard_Boolean isInscribe =
3033               ((aUSurf2f - aU2[i]) <= aTol2D) && ((aU2[i] - aUSurf2l) <= aTol2D) &&
3034               ((aVSurf1f - aV1[i]) <= aTol2D) && ((aV1[i] - aVSurf1l) <= aTol2D) &&
3035               ((aVSurf2f - aV2[i]) <= aTol2D) && ((aV2[i] - aVSurf2l) <= aTol2D);
3036
3037           //isVIntersect == TRUE if intersection line intersects two (!)
3038           //V-bounds of cylinder (1st or 2nd - no matter)
3039           const Standard_Boolean isVIntersect =
3040               (((aVSurf1f - aV1[i])*(aVSurf1f - aV1Prev[i]) < RealSmall()) &&
3041                 ((aVSurf1l - aV1[i])*(aVSurf1l - aV1Prev[i]) < RealSmall())) ||
3042               (((aVSurf2f - aV2[i])*(aVSurf2f - aV2Prev[i]) < RealSmall()) &&
3043                 ((aVSurf2l - aV2[i])*(aVSurf2l - aV2Prev[i]) < RealSmall()));
3044
3045           //isFound1 == TRUE if intersection line intersects V-bounds
3046           //  (First or Last - no matter) of the 1st cylynder
3047           //isFound2 == TRUE if intersection line intersects V-bounds
3048           //  (First or Last - no matter) of the 2nd cylynder
3049           Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
3050           Standard_Boolean isForce = Standard_False;
3051
3052           if (aWLFindStatus[i] == WLFStatus_Absent)
3053           {
3054             if (((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1 - aUSurf1l) < aTol2D))
3055             {
3056               isForce = Standard_True;
3057             }
3058           }
3059
3060           theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
3061                                  aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i, isForce,
3062                                  isFound1, isFound2);
3063
3064           const Standard_Boolean isPrevVBound = !isVIntersect &&
3065                                                 ((Abs(aV1Prev[i] - aVSurf1f) <= aTol2D) ||
3066                                                  (Abs(aV1Prev[i] - aVSurf1l) <= aTol2D) ||
3067                                                  (Abs(aV2Prev[i] - aVSurf2f) <= aTol2D) ||
3068                                                  (Abs(aV2Prev[i] - aVSurf2l) <= aTol2D));
3069
3070           aV1Prev[i] = aV1[i];
3071           aV2Prev[i] = aV2[i];
3072
3073           if ((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
3074           {
3075             aWLFindStatus[i] = WLFStatus_Broken; //start a new line
3076           }
3077           else if (isInscribe)
3078           {
3079             if ((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
3080             {
3081               aWLFindStatus[i] = WLFStatus_Exist;
3082             }
3083
3084             if ((aWLFindStatus[i] != WLFStatus_Broken) ||
3085               (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
3086             {
3087               if (aWLine[i]->NbPnts() > 0)
3088               {
3089                 Standard_Real aU2p = 0.0, aV2p = 0.0;
3090                 if (isReversed)
3091                   aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
3092                 else
3093                   aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
3094
3095                 const Standard_Real aDelta = aU2[i] - aU2p;
3096
3097                 if (2.0 * Abs(aDelta) > aPeriod)
3098                 {
3099                   if (aDelta > 0.0)
3100                   {
3101                     aU2[i] -= aPeriod;
3102                   }
3103                   else
3104                   {
3105                     aU2[i] += aPeriod;
3106                   }
3107                 }
3108               }
3109
3110               if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
3111                                 gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
3112                                 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3113                                 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
3114                                 aWLine[i]->Curve(), i, aTol3D, aTol2D, isForce))
3115               {
3116                 if (aWLFindStatus[i] == WLFStatus_Absent)
3117                 {
3118                   aWLFindStatus[i] = WLFStatus_Exist;
3119                 }
3120               }
3121               else if (!isFound1 && !isFound2)
3122               {//We do not add any point while doing this iteration
3123                 if (aWLFindStatus[i] == WLFStatus_Exist)
3124                 {
3125                   aWLFindStatus[i] = WLFStatus_Broken;
3126                 }
3127               }
3128             }
3129           }
3130           else
3131           {//We do not add any point while doing this iteration
3132             if (aWLFindStatus[i] == WLFStatus_Exist)
3133             {
3134               aWLFindStatus[i] = WLFStatus_Broken;
3135             }
3136           }
3137
3138           if (aWLFindStatus[i] == WLFStatus_Broken)
3139             isBroken = Standard_True;
3140         }//for(Standard_Integer i = 0; i < aNbWLines; i++)
3141
3142         if (isBroken)
3143         {//current lines are filled. Go to the next lines
3144           anUf = anU1;
3145
3146           Standard_Boolean isAdded = Standard_True;
3147
3148           for (Standard_Integer i = 0; i < aNbWLines; i++)
3149           {
3150             if (isAddingWLEnabled[i])
3151             {
3152               continue;
3153             }
3154
3155             isAdded = Standard_False;
3156
3157             Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
3158
3159             theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
3160                                    aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i,
3161                                    Standard_False, isFound1, isFound2);
3162
3163             if (isFound1 || isFound2)
3164             {
3165               isAdded = Standard_True;
3166             }
3167
3168             if (aWLine[i]->NbPnts() > 0)
3169             {
3170               Standard_Real aU2p = 0.0, aV2p = 0.0;
3171               if (isReversed)
3172                 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
3173               else
3174                 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
3175
3176               const Standard_Real aDelta = aU2[i] - aU2p;
3177
3178               if (2 * Abs(aDelta) > aPeriod)
3179               {
3180                 if (aDelta > 0.0)
3181                 {
3182                   aU2[i] -= aPeriod;
3183                 }
3184                 else
3185                 {
3186                   aU2[i] += aPeriod;
3187                 }
3188               }
3189             }
3190
3191             if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
3192                               Standard_True, gp_Pnt2d(anU1, aV1[i]),
3193                               gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3194                               aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3195                               aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
3196                               i, aTol3D, aTol2D, Standard_False))
3197             {
3198               isAdded = Standard_True;
3199             }
3200           }
3201
3202           if (!isAdded)
3203           {
3204             //Before breaking WL, we must complete it correctly
3205             //(e.g. to prolong to the surface boundary).
3206             //Therefore, we take the point last added in some WL
3207             //(have maximal U1-parameter) and try to add it in
3208             //the current WL.
3209             Standard_Real anUmaxAdded = RealFirst();
3210
3211             {
3212               Standard_Boolean isChanged = Standard_False;
3213               for (Standard_Integer i = 0; i < aNbWLines; i++)
3214               {
3215                 if ((aWLFindStatus[i] == WLFStatus_Absent) || (aWLine[i]->NbPnts() == 0))
3216                   continue;
3217
3218                 Standard_Real aU1c = 0.0, aV1c = 0.0;
3219                 if (isReversed)
3220                   aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
3221                 else
3222                   aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
3223
3224                 anUmaxAdded = Max(anUmaxAdded, aU1c);
3225                 isChanged = Standard_True;
3226               }
3227
3228               if (!isChanged)
3229               { //If anUmaxAdded were not changed in previous cycle then
3230                 //we would break existing WLines.
3231                 break;
3232               }
3233             }
3234
3235             for (Standard_Integer i = 0; i < aNbWLines; i++)
3236             {
3237               if (isAddingWLEnabled[i])
3238               {
3239                 continue;
3240               }
3241
3242               ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
3243                 aU2[i], aV1[i], aV2[i]);
3244
3245               AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
3246                              Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
3247                              gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3248                              aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3249                              aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
3250                              i, aTol3D, aTol2D, Standard_False);
3251             }
3252           }
3253
3254           break;
3255         }
3256
3257         //Step computing
3258
3259         {
3260           //Step of aU1-parameter is computed adaptively. The algorithm
3261           //aims to provide given aDeltaV1 and aDeltaV2 values (if it is
3262           //possible because the intersection line can go along V-isoline)
3263           //in every iteration. It allows avoiding "flying" intersection
3264           //points too far each from other (see issue #24915).
3265
3266           const Standard_Real aDeltaV1 = aRangeS1.Delta() / IntToReal(aNbPoints),
3267                               aDeltaV2 = aRangeS2.Delta() / IntToReal(aNbPoints);
3268
3269           math_Matrix aMatr(1, 3, 1, 5);
3270
3271           Standard_Real aMinUexp = RealLast();
3272
3273           for (Standard_Integer i = 0; i < aNbWLines; i++)
3274           {
3275             if (aTol2D < (anUexpect[i] - anU1))
3276             {
3277               continue;
3278             }
3279
3280             if (aWLFindStatus[i] == WLFStatus_Absent)
3281             {
3282               anUexpect[i] += aStepMax;
3283               aMinUexp = Min(aMinUexp, anUexpect[i]);
3284               continue;
3285             }
3286
3287             Standard_Real aStepTmp = aStepMax;
3288
3289             const Standard_Real aSinU1 = sin(anU1),
3290                                 aCosU1 = cos(anU1),
3291                                 aSinU2 = sin(aU2[i]),
3292                                 aCosU2 = cos(aU2[i]);
3293
3294             aMatr.SetCol(1, anEquationCoeffs.mVecC1);
3295             aMatr.SetCol(2, anEquationCoeffs.mVecC2);
3296             aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
3297             aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
3298             aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
3299                             anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
3300                             anEquationCoeffs.mVecD);
3301
3302             //The main idea is in solving of linearized system (2)
3303             //(see description to ComputationMethods class) in order to find new U1-value
3304             //to provide new value V1 or V2, which differs from current one by aDeltaV1 or
3305             //aDeltaV2 respectively.
3306
3307             //While linearizing, following Taylor formulas are used:
3308             //    cos(x0+dx) = cos(x0) - sin(x0)*dx
3309             //    sin(x0+dx) = sin(x0) + cos(x0)*dx
3310
3311             //Consequently cos(U1), cos(U2), sin(U1) and sin(U2) in the system (2)
3312             //must be substituted by corresponding values.
3313
3314             //ATTENTION!!!
3315             //The solution is approximate. More over, all requirements to the
3316             //linearization must be satisfied in order to obtain quality result.
3317
3318             if (!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aStepTmp))
3319             {
3320               //To avoid cycling-up
3321               anUexpect[i] += aStepMax;
3322               aMinUexp = Min(aMinUexp, anUexpect[i]);
3323
3324               continue;
3325             }
3326
3327             if (aStepTmp < aStepMin)
3328               aStepTmp = aStepMin;
3329
3330             if (aStepTmp > aStepMax)
3331               aStepTmp = aStepMax;
3332
3333             anUexpect[i] = anU1 + aStepTmp;
3334             aMinUexp = Min(aMinUexp, anUexpect[i]);
3335           }
3336
3337           anU1 = aMinUexp;
3338         }
3339
3340         if (Precision::PConfusion() >= (anUl - anU1))
3341           anU1 = anUl;
3342
3343         anUf = anU1;
3344
3345         for (Standard_Integer i = 0; i < aNbWLines; i++)
3346         {
3347           if (aWLine[i]->NbPnts() != 1)
3348             isAddedIntoWL[i] = Standard_False;
3349
3350           if (anU1 == anUl)
3351           {//strictly equal. Tolerance is considered above.
3352             anUexpect[i] = anUl;
3353           }
3354         }
3355       }
3356
3357       for (Standard_Integer i = 0; i < aNbWLines; i++)
3358       {
3359         if ((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
3360         {
3361           isTheEmpty = Standard_False;
3362           Standard_Real u1, v1, u2, v2;
3363           aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
3364           IntPatch_Point aP;
3365           aP.SetParameter(u1);
3366           aP.SetParameters(u1, v1, u2, v2);
3367           aP.SetTolerance(aTol3D);
3368           aP.SetValue(aWLine[i]->Point(1).Value());
3369
3370           theSPnt.Append(aP);
3371         }
3372         else if (aWLine[i]->NbPnts() > 1)
3373         {
3374           Standard_Boolean isGood = Standard_True;
3375
3376           if (aWLine[i]->NbPnts() == 2)
3377           {
3378             const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
3379             const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
3380
3381             if (aPf.IsSame(aPl, Precision::Confusion()))
3382               isGood = Standard_False;
3383           }
3384
3385           if (isGood)
3386           {
3387             isTheEmpty = Standard_False;
3388             isAddedIntoWL[i] = Standard_True;
3389             SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
3390                                  anEquationCoeffs, i, aNbPoints, 1,
3391               aWLine[i]->NbPnts(), aTol2D, aPeriod,
3392               isReversed);
3393
3394             aWLine[i]->ComputeVertexParameters(aTol3D);
3395             theSlin.Append(aWLine[i]);
3396           }
3397         }
3398         else
3399         {
3400           isAddedIntoWL[i] = Standard_False;
3401         }
3402
3403 #ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
3404         //aWLine[i]->Dump();
3405 #endif
3406       }
3407     }
3408   }
3409
3410
3411   //Delete the points in theSPnt, which
3412   //lie at least in one of the line in theSlin.
3413   for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
3414   {
3415     for (Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
3416     {
3417       Handle(IntPatch_WLine) aWLine1(Handle(IntPatch_WLine)::
3418         DownCast(theSlin.Value(aNbLin)));
3419
3420       const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
3421       const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
3422
3423       const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
3424       if (aPntCur.IsSame(aPntFWL1, Precision::Confusion()) ||
3425         aPntCur.IsSame(aPntLWL1, Precision::Confusion()))
3426       {
3427         theSPnt.Remove(aNbPnt);
3428         aNbPnt--;
3429         break;
3430       }
3431     }
3432   }
3433
3434   //Try to add new points in the neighborhood of existing point
3435   for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
3436   {
3437     // Standard algorithm (implemented above) could not find any
3438     // continuous curve in neighborhood of aPnt2S (e.g. because
3439     // this curve is too small; see tests\bugs\modalg_5\bug25292_35 and _36).
3440     // Here, we will try to find several new points nearer to aPnt2S.
3441
3442     // The algorithm below tries to find two points in every
3443     // intervals [u1 - aStepMax, u1] and [u1, u1 + aStepMax]
3444     // and every new point will be in maximal distance from
3445     // u1. If these two points exist they will be joined
3446     // by the intersection curve.
3447
3448     const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
3449
3450     Standard_Real u1, v1, u2, v2;
3451     aPnt2S.Parameters(u1, v1, u2, v2);
3452
3453     Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
3454     Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
3455
3456     //Define the index of WLine, which lies the point aPnt2S in.
3457     Standard_Integer anIndex = 0;
3458
3459     Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
3460     if (isReversed)
3461     {
3462       anUf = Max(u2 - aStepMax, aUSurf1f);
3463       anUl = Min(u2 + aStepMax, aUSurf1l);
3464       aCurU2 = u1;
3465     }
3466     else
3467     {
3468       anUf = Max(u1 - aStepMax, aUSurf1f);
3469       anUl = Min(u1 + aStepMax, aUSurf1l);
3470       aCurU2 = u2;
3471     }
3472
3473     const Standard_Real anUinf = anUf, anUsup = anUl, anUmid = 0.5*(anUf + anUl);
3474
3475     {
3476       //Find the value of anIndex variable.
3477       Standard_Real aDelta = RealLast();
3478       for (Standard_Integer i = 0; i < aNbWLines; i++)
3479       {
3480         Standard_Real anU2t = 0.0;
3481         if (!ComputationMethods::CylCylComputeParameters(anUmid, i, anEquationCoeffs, anU2t))
3482           continue;
3483
3484         Standard_Real aDU2 = fmod(Abs(anU2t - aCurU2), aPeriod);
3485         aDU2 = Min(aDU2, Abs(aDU2 - aPeriod));
3486         if (aDU2 < aDelta)
3487         {
3488           aDelta = aDU2;
3489           anIndex = i;
3490         }
3491       }
3492     }
3493
3494     // Bisection method is used in order to find every new point.
3495     // I.e. if we need to find intersection point in the interval [anUinf, anUmid]
3496     // we check the point anUC = 0.5*(anUinf+anUmid). If it is an intersection point
3497     // we try to find another point in the interval [anUinf, anUC] (because we find the point in
3498     // maximal distance from anUmid). If it is not then we try to find another point in the
3499     // interval [anUC, anUmid]. Next iterations will be made analogically.
3500     // When we find intersection point in the interval [anUmid, anUsup] we try to find
3501     // another point in the interval [anUC, anUsup] if anUC is intersection point and
3502     // in the interval [anUmid, anUC], otherwise.
3503
3504     Standard_Real anAddedPar[2] = { anUmid, anUmid };
3505
3506     for (Standard_Integer aParID = 0; aParID < 2; aParID++)
3507     {
3508       if (aParID == 0)
3509       {
3510         anUf = anUinf;
3511         anUl = anUmid;
3512       }
3513       else // if(aParID == 1)
3514       {
3515         anUf = anUmid;
3516         anUl = anUsup;
3517       }
3518
3519       Standard_Real &aPar1 = (aParID == 0) ? anUf : anUl,
3520                     &aPar2 = (aParID == 0) ? anUl : anUf;
3521
3522       while (Abs(aPar2 - aPar1) > aStepMin)
3523       {
3524         Standard_Real anUC = 0.5*(anUf + anUl);
3525         Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
3526         Standard_Boolean isDone = ComputationMethods::
3527                 CylCylComputeParameters(anUC, anIndex, anEquationCoeffs, aU2, aV1, aV2);
3528
3529         if (isDone)
3530         {
3531           if (Abs(aV1 - aVSurf1f) <= aTol2D)
3532             aV1 = aVSurf1f;
3533
3534           if (Abs(aV1 - aVSurf1l) <= aTol2D)
3535             aV1 = aVSurf1l;
3536
3537           if (Abs(aV2 - aVSurf2f) <= aTol2D)
3538             aV2 = aVSurf2f;
3539
3540           if (Abs(aV2 - aVSurf2l) <= aTol2D)
3541             aV2 = aVSurf2l;
3542
3543           isDone = AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
3544                                   Standard_True, gp_Pnt2d(anUC, aV1), gp_Pnt2d(aU2, aV2),
3545                                   aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3546                                   aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
3547                                   aPeriod, aWLine->Curve(), anIndex, aTol3D,
3548                                   aTol2D, Standard_False, Standard_True);
3549         }
3550
3551         if (isDone)
3552         {
3553           anAddedPar[0] = Min(anAddedPar[0], anUC);
3554           anAddedPar[1] = Max(anAddedPar[1], anUC);
3555           aPar2 = anUC;
3556         }
3557         else
3558         {
3559           aPar1 = anUC;
3560         }
3561       }
3562     }
3563
3564     //Fill aWLine by additional points
3565     if (anAddedPar[1] - anAddedPar[0] > aStepMin)
3566     {
3567       for (Standard_Integer aParID = 0; aParID < 2; aParID++)
3568       {
3569         Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
3570         ComputationMethods::CylCylComputeParameters(anAddedPar[aParID], anIndex,
3571                                                   anEquationCoeffs, aU2, aV1, aV2);
3572
3573         AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
3574                         gp_Pnt2d(anAddedPar[aParID], aV1), gp_Pnt2d(aU2, aV2),
3575                         aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3576                         aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod, aWLine->Curve(),
3577                         anIndex, aTol3D, aTol2D, Standard_False, Standard_False);
3578       }
3579
3580       SeekAdditionalPoints(aQuad1, aQuad2, aWLine->Curve(),
3581                             anEquationCoeffs, anIndex, aNbMinPoints,
3582                             1, aWLine->NbPnts(), aTol2D, aPeriod,
3583                             isReversed);
3584
3585       aWLine->ComputeVertexParameters(aTol3D);
3586       theSlin.Append(aWLine);
3587
3588       theSPnt.Remove(aNbPnt);
3589       aNbPnt--;
3590     }
3591   }
3592
3593   return IntPatch_ImpImpIntersection::IntStatus_OK;
3594 }
3595
3596 //=======================================================================
3597 //function : IntCyCy
3598 //purpose  :
3599 //=======================================================================
3600 IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1,
3601                                                const IntSurf_Quadric& theQuad2,
3602                                                const Standard_Real theTol3D,
3603                                                const Standard_Real theTol2D,
3604                                                const Bnd_Box2d& theUVSurf1,
3605                                                const Bnd_Box2d& theUVSurf2,
3606                                                Standard_Boolean& isTheEmpty,
3607                                                Standard_Boolean& isTheSameSurface,
3608                                                Standard_Boolean& isTheMultiplePoint,
3609                                                IntPatch_SequenceOfLine& theSlin,
3610                                                IntPatch_SequenceOfPoint& theSPnt)
3611 {
3612   isTheEmpty = Standard_True;
3613   isTheSameSurface = Standard_False;
3614   isTheMultiplePoint = Standard_False;
3615   theSlin.Clear();
3616   theSPnt.Clear();
3617
3618   const gp_Cylinder aCyl1 = theQuad1.Cylinder(),
3619                     aCyl2 = theQuad2.Cylinder();
3620
3621   IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
3622
3623   if (!anInter.IsDone())
3624   {
3625     return IntPatch_ImpImpIntersection::IntStatus_Fail;
3626   }
3627
3628   if(anInter.TypeInter() != IntAna_NoGeometricSolution)
3629   {
3630     if (CyCyAnalyticalIntersect(theQuad1, theQuad2, anInter,
3631                                 theTol3D, isTheEmpty,
3632                                 isTheSameSurface, isTheMultiplePoint,
3633                                 theSlin, theSPnt))
3634     {
3635       return IntPatch_ImpImpIntersection::IntStatus_OK;
3636     }
3637   }
3638
3639   //Here, intersection line is not an analytical curve(line, circle, ellipsis etc.)
3640
3641   Standard_Real aUSBou[2][2], aVSBou[2][2]; //const
3642
3643   theUVSurf1.Get(aUSBou[0][0], aVSBou[0][0], aUSBou[0][1], aVSBou[0][1]);
3644   theUVSurf2.Get(aUSBou[1][0], aVSBou[1][0], aUSBou[1][1], aVSBou[1][1]);
3645
3646   const Standard_Real aPeriod = 2.0*M_PI;
3647   const Standard_Integer aNbWLines = 2;
3648
3649   const ComputationMethods::stCoeffsValue anEquationCoeffs1(aCyl1, aCyl2);
3650   const ComputationMethods::stCoeffsValue anEquationCoeffs2(aCyl2, aCyl1);
3651
3652   //Boundaries.
3653   //Intersection result can include two non-connected regions
3654   //(see WorkWithBoundaries::BoundariesComputing(...) method).
3655   const Standard_Integer aNbOfBoundaries = 2;
3656   Bnd_Range anURange[2][aNbOfBoundaries];   //const
3657
3658   if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs1, aPeriod, anURange[0]))
3659     return IntPatch_ImpImpIntersection::IntStatus_OK;
3660
3661   if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs2, aPeriod, anURange[1]))
3662     return IntPatch_ImpImpIntersection::IntStatus_OK;
3663
3664   //anURange[*] can be in different periodic regions in
3665   //compare with First-Last surface. E.g. the surface
3666   //is full cylinder [0, 2*PI] but anURange is [5, 7].
3667   //Trivial common range computation returns [5, 2*PI] and
3668   //its summary length is 2*PI-5 == 1.28... only. That is wrong.
3669   //This problem can be solved by the following
3670   //algorithm:
3671   // 1. split anURange[*] by the surface boundary;
3672   // 2. shift every new range in order to inscribe it
3673   //      in [Ufirst, Ulast] of cylinder;
3674   // 3. consider only common ranges between [Ufirst, Ulast]
3675   //    and new ranges.
3676   //
3677   // In above example, we obtain following:
3678   // 1. two ranges: [5, 2*PI] and [2*PI, 7];
3679   // 2. after shifting: [5, 2*PI] and [0, 7-2*PI];
3680   // 3. Common ranges: ([5, 2*PI] and [0, 2*PI]) == [5, 2*PI],
3681   //                   ([0, 7-2*PI] and [0, 2*PI]) == [0, 7-2*PI];
3682   // 4. Their summary length is (2*PI-5)+(7-2*PI-0)==7-5==2 ==> GOOD.
3683
3684   Standard_Real aSumRange[2] = { 0.0, 0.0 };
3685   Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
3686   for (Standard_Integer aCID = 0; aCID < 2; aCID++)
3687   {
3688     anAlloc->Reset();
3689     NCollection_List<Bnd_Range> aListOfRng(anAlloc);
3690
3691     aListOfRng.Append(anURange[aCID][0]);
3692     aListOfRng.Append(anURange[aCID][1]);
3693
3694     const Standard_Real aSplitArr[3] = {aUSBou[aCID][0], aUSBou[aCID][1], 0.0};
3695
3696     NCollection_List<Bnd_Range>::Iterator anITrRng;
3697     for (Standard_Integer aSInd = 0; aSInd < 3; aSInd++)
3698     {
3699       NCollection_List<Bnd_Range> aLstTemp(aListOfRng);
3700       aListOfRng.Clear();
3701       for (anITrRng.Init(aLstTemp); anITrRng.More(); anITrRng.Next())
3702       {
3703         Bnd_Range& aRng = anITrRng.Value();
3704         aRng.Split(aSplitArr[aSInd], aListOfRng, aPeriod);
3705       }
3706     }
3707
3708     anITrRng.Init(aListOfRng);
3709     for (; anITrRng.More(); anITrRng.Next())
3710     {
3711       Bnd_Range& aCurrRange = anITrRng.Value();
3712
3713       Bnd_Range aBoundR;
3714       aBoundR.Add(aUSBou[aCID][0]);
3715       aBoundR.Add(aUSBou[aCID][1]);
3716
3717       if (!InscribeInterval(aUSBou[aCID][0], aUSBou[aCID][1],
3718                                           aCurrRange, theTol2D, aPeriod))
3719       {
3720         //If aCurrRange does not have common block with
3721         //[Ufirst, Ulast] of cylinder then we will try
3722         //to inscribe [Ufirst, Ulast] in the boundaries of aCurrRange.
3723         Standard_Real aF = 1.0, aL = 0.0;
3724         if (!aCurrRange.GetBounds(aF, aL))
3725           continue;
3726
3727         if ((aL < aUSBou[aCID][0]))
3728         {
3729           aCurrRange.Shift(aPeriod);
3730         }
3731         else if (aF > aUSBou[aCID][1])
3732         {
3733           aCurrRange.Shift(-aPeriod);
3734         }
3735       }
3736
3737       aBoundR.Common(aCurrRange);
3738
3739       const Standard_Real aDelta = aBoundR.Delta();
3740
3741       if (aDelta > 0.0)
3742       {
3743         aSumRange[aCID] += aDelta;
3744       }
3745     }
3746   }
3747
3748   //The bigger range the bigger number of points in Walking-line (WLine)
3749   //we will be able to add and consequently we will obtain more
3750   //precise intersection line.
3751   //Every point of WLine is determined as function from U1-parameter,
3752   //where U1 is U-parameter on 1st quadric.
3753   //Therefore, we should use quadric with bigger range as 1st parameter
3754   //in IntCyCy() function.
3755   //On the other hand, there is no point in reversing in case of
3756   //analytical intersection (when result is line, ellipse, point...).
3757   //This result is independent of the arguments order.
3758   const Standard_Boolean isToReverse = (aSumRange[1] > aSumRange[0]);
3759
3760   if (isToReverse)
3761   {
3762     const WorkWithBoundaries aBoundWork(theQuad2, theQuad1, anEquationCoeffs2,
3763                                         theUVSurf2, theUVSurf1, aNbWLines,
3764                                         aPeriod, theTol3D, theTol2D, Standard_True);
3765
3766     return CyCyNoGeometric(aCyl2, aCyl1, aBoundWork, anURange[1], aNbOfBoundaries,
3767                               isTheEmpty, theSlin, theSPnt);
3768   }
3769   else
3770   {
3771     const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs1,
3772                                         theUVSurf1, theUVSurf2, aNbWLines,
3773                                         aPeriod, theTol3D, theTol2D, Standard_False);
3774
3775     return CyCyNoGeometric(aCyl1, aCyl2, aBoundWork, anURange[0], aNbOfBoundaries,
3776                               isTheEmpty, theSlin, theSPnt);
3777   }
3778 }
3779
3780 //=======================================================================
3781 //function : IntCySp
3782 //purpose  :
3783 //=======================================================================
3784 Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
3785                          const IntSurf_Quadric& Quad2,
3786                          const Standard_Real Tol,
3787                          const Standard_Boolean Reversed,
3788                          Standard_Boolean& Empty,
3789                          Standard_Boolean& Multpoint,
3790                          IntPatch_SequenceOfLine& slin,
3791                          IntPatch_SequenceOfPoint& spnt)
3792
3793 {
3794   Standard_Integer i;
3795
3796   IntSurf_TypeTrans trans1,trans2;
3797   IntAna_ResultType typint;
3798   IntPatch_Point ptsol;
3799   gp_Circ cirsol;
3800
3801   gp_Cylinder Cy;
3802   gp_Sphere Sp;
3803
3804   if (!Reversed) {
3805     Cy = Quad1.Cylinder();
3806     Sp = Quad2.Sphere();
3807   }
3808   else {
3809     Cy = Quad2.Cylinder();
3810     Sp = Quad1.Sphere();
3811   }
3812   IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
3813
3814   if (!inter.IsDone()) {return Standard_False;}
3815
3816   typint = inter.TypeInter();
3817   Standard_Integer NbSol = inter.NbSolutions();
3818   Empty = Standard_False;
3819
3820   switch (typint) {
3821
3822   case IntAna_Empty :
3823     {
3824       Empty = Standard_True;
3825     }
3826     break;
3827
3828   case IntAna_Point :
3829     {
3830       gp_Pnt psol(inter.Point(1));
3831       Standard_Real U1,V1,U2,V2;
3832       Quad1.Parameters(psol,U1,V1);
3833       Quad2.Parameters(psol,U2,V2);
3834       ptsol.SetValue(psol,Tol,Standard_True);
3835       ptsol.SetParameters(U1,V1,U2,V2);
3836       spnt.Append(ptsol);
3837     }
3838     break;
3839
3840   case IntAna_Circle:
3841     {
3842       cirsol = inter.Circle(1);
3843       gp_Vec Tgt;
3844       gp_Pnt ptref;
3845       ElCLib::D1(0.,cirsol,ptref,Tgt);
3846
3847       if (NbSol == 1) {
3848         gp_Vec TestCurvature(ptref,Sp.Location());
3849         gp_Vec Normsp,Normcyl;
3850         if (!Reversed) {
3851           Normcyl = Quad1.Normale(ptref);
3852           Normsp  = Quad2.Normale(ptref);
3853         }
3854         else {
3855           Normcyl = Quad2.Normale(ptref);
3856           Normsp  = Quad1.Normale(ptref);
3857         }
3858
3859         IntSurf_Situation situcyl;
3860         IntSurf_Situation situsp;
3861
3862         if (Normcyl.Dot(TestCurvature) > 0.) {
3863           situsp = IntSurf_Outside;
3864           if (Normsp.Dot(Normcyl) > 0.) {
3865             situcyl = IntSurf_Inside;
3866           }
3867           else {
3868             situcyl = IntSurf_Outside;
3869           }
3870         }
3871         else {
3872           situsp = IntSurf_Inside;
3873           if (Normsp.Dot(Normcyl) > 0.) {
3874             situcyl = IntSurf_Outside;
3875           }
3876           else {
3877             situcyl = IntSurf_Inside;
3878           }
3879         }
3880         Handle(IntPatch_GLine) glig;
3881         if (!Reversed) {
3882           glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
3883         }
3884         else {
3885           glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
3886         }
3887         slin.Append(glig);
3888       }
3889       else {
3890         if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
3891           trans1 = IntSurf_Out;
3892           trans2 = IntSurf_In;
3893         }
3894         else {
3895           trans1 = IntSurf_In;
3896           trans2 = IntSurf_Out;
3897         }
3898         Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3899         slin.Append(glig);
3900
3901         cirsol = inter.Circle(2);
3902         ElCLib::D1(0.,cirsol,ptref,Tgt);
3903         Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
3904         if(qwe> 0.0000001) {
3905           trans1 = IntSurf_Out;
3906           trans2 = IntSurf_In;
3907         }
3908         else if(qwe<-0.0000001) {
3909           trans1 = IntSurf_In;
3910           trans2 = IntSurf_Out;
3911         }
3912         else {
3913           trans1=trans2=IntSurf_Undecided;
3914         }
3915         glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3916         slin.Append(glig);
3917       }
3918     }
3919     break;
3920
3921   case IntAna_NoGeometricSolution:
3922     {
3923       gp_Pnt psol;
3924       Standard_Real U1,V1,U2,V2;
3925       IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
3926       if (!anaint.IsDone()) {
3927         return Standard_False;
3928       }
3929
3930       if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
3931         Empty = Standard_True;
3932       }
3933       else {
3934
3935         NbSol = anaint.NbPnt();
3936         for (i = 1; i <= NbSol; i++) {
3937           psol = anaint.Point(i);
3938           Quad1.Parameters(psol,U1,V1);
3939           Quad2.Parameters(psol,U2,V2);
3940           ptsol.SetValue(psol,Tol,Standard_True);
3941           ptsol.SetParameters(U1,V1,U2,V2);
3942           spnt.Append(ptsol);
3943         }
3944
3945         gp_Pnt ptvalid,ptf,ptl;
3946         gp_Vec tgvalid;
3947         Standard_Real first,last,para;
3948         IntAna_Curve curvsol;
3949         Standard_Boolean tgfound;
3950         Standard_Integer kount;
3951
3952         NbSol = anaint.NbCurve();
3953         for (i = 1; i <= NbSol; i++) {
3954           curvsol = anaint.Curve(i);
3955           curvsol.Domain(first,last);
3956           ptf = curvsol.Value(first);
3957           ptl = curvsol.Value(last);
3958
3959           para = last;
3960           kount = 1;
3961           tgfound = Standard_False;
3962
3963           while (!tgfound) {
3964             para = (1.123*first + para)/2.123;
3965             tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3966             if (!tgfound) {
3967               kount ++;
3968               tgfound = kount > 5;
3969             }
3970           }
3971           Handle(IntPatch_ALine) alig;
3972           if (kount <= 5) {
3973             Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
3974                                                  Quad1.Normale(ptvalid));
3975             if(qwe> 0.00000001) {
3976               trans1 = IntSurf_Out;
3977               trans2 = IntSurf_In;
3978             }
3979             else if(qwe<-0.00000001) {
3980               trans1 = IntSurf_In;
3981               trans2 = IntSurf_Out;
3982             }
3983             else {
3984               trans1=trans2=IntSurf_Undecided;
3985             }
3986             alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3987           }
3988           else {
3989             alig = new IntPatch_ALine(curvsol,Standard_False);
3990           }
3991           Standard_Boolean TempFalse1a = Standard_False;
3992           Standard_Boolean TempFalse2a = Standard_False;
3993
3994           //-- ptf et ptl : points debut et fin de alig
3995
3996           ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
3997                         TempFalse2a,ptl,last,Multpoint,Tol);
3998           slin.Append(alig);
3999         } //-- boucle sur les lignes
4000       } //-- solution avec au moins une lihne
4001     }
4002     break;
4003
4004   default:
4005     {
4006       return Standard_False;
4007     }
4008   }
4009   return Standard_True;
4010 }
4011 //=======================================================================
4012 //function : IntCyCo
4013 //purpose  :
4014 //=======================================================================
4015 Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
4016                          const IntSurf_Quadric& Quad2,
4017                          const Standard_Real Tol,
4018                          const Standard_Boolean Reversed,
4019                          Standard_Boolean& Empty,
4020                          Standard_Boolean& Multpoint,
4021                          IntPatch_SequenceOfLine& slin,
4022                          IntPatch_SequenceOfPoint& spnt)
4023
4024 {
4025   IntPatch_Point ptsol;
4026
4027   Standard_Integer i;
4028
4029   IntSurf_TypeTrans trans1,trans2;
4030   IntAna_ResultType typint;
4031   gp_Circ cirsol;
4032
4033   gp_Cylinder Cy;
4034   gp_Cone     Co;
4035
4036   if (!Reversed) {
4037     Cy = Quad1.Cylinder();
4038     Co = Quad2.Cone();
4039   }
4040   else {
4041     Cy = Quad2.Cylinder();
4042     Co = Quad1.Cone();
4043   }
4044   IntAna_QuadQuadGeo inter(Cy,Co,Tol);
4045
4046   if (!inter.IsDone()) {return Standard_False;}
4047
4048   typint = inter.TypeInter();
4049   Standard_Integer NbSol = inter.NbSolutions();
4050   Empty = Standard_False;
4051
4052   switch (typint) {
4053
4054   case IntAna_Empty : {
4055     Empty = Standard_True;
4056   }
4057     break;
4058
4059   case IntAna_Point :{
4060     gp_Pnt psol(inter.Point(1));
4061     Standard_Real U1,V1,U2,V2;
4062     Quad1.Parameters(psol,U1,V1);
4063     Quad1.Parameters(psol,U2,V2);
4064     ptsol.SetValue(psol,Tol,Standard_True);
4065     ptsol.SetParameters(U1,V1,U2,V2);
4066     spnt.Append(ptsol);
4067   }
4068     break;
4069
4070   case IntAna_Circle:  {
4071     gp_Vec Tgt;
4072     gp_Pnt ptref;
4073     Standard_Integer j;
4074     Standard_Real qwe;
4075     //
4076     for(j=1; j<=2; ++j) {
4077       cirsol = inter.Circle(j);
4078       ElCLib::D1(0.,cirsol,ptref,Tgt);
4079       qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
4080       if(qwe> 0.00000001) {
4081         trans1 = IntSurf_Out;
4082         trans2 = IntSurf_In;
4083       }
4084       else if(qwe<-0.00000001) {
4085         trans1 = IntSurf_In;
4086         trans2 = IntSurf_Out;
4087       }
4088       else {
4089         trans1=trans2=IntSurf_Undecided;
4090       }
4091       Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
4092       slin.Append(glig);
4093     }
4094   }
4095     break;
4096
4097   case IntAna_NoGeometricSolution:    {
4098     gp_Pnt psol;
4099     Standard_Real U1,V1,U2,V2;
4100     IntAna_IntQuadQuad anaint(Cy,Co,Tol);
4101     if (!anaint.IsDone()) {
4102       return Standard_False;
4103     }
4104
4105     if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
4106       Empty = Standard_True;
4107     }
4108     else {
4109       NbSol = anaint.NbPnt();
4110       for (i = 1; i <= NbSol; i++) {
4111         psol = anaint.Point(i);
4112         Quad1.Parameters(psol,U1,V1);
4113         Quad2.Parameters(psol,U2,V2);
4114         ptsol.SetValue(psol,Tol,Standard_True);
4115         ptsol.SetParameters(U1,V1,U2,V2);
4116         spnt.Append(ptsol);
4117       }
4118
4119       gp_Pnt ptvalid, ptf, ptl;
4120       gp_Vec tgvalid;
4121       //
4122       Standard_Real first,last,para;
4123       Standard_Boolean tgfound,firstp,lastp,kept;
4124       Standard_Integer kount;
4125       //
4126       //
4127       //IntAna_Curve curvsol;
4128       IntAna_Curve aC;
4129       IntAna_ListOfCurve aLC;
4130       IntAna_ListIteratorOfListOfCurve aIt;
4131
4132       //
4133       NbSol = anaint.NbCurve();
4134       for (i = 1; i <= NbSol; ++i) {
4135         kept = Standard_False;
4136         //curvsol = anaint.Curve(i);
4137         aC=anaint.Curve(i);
4138         aLC.Clear();
4139         ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
4140         //
4141         aIt.Initialize(aLC);
4142         for (; aIt.More(); aIt.Next()) {
4143           IntAna_Curve& curvsol=aIt.Value();
4144           //
4145           curvsol.Domain(first, last);
4146           firstp = !curvsol.IsFirstOpen();
4147           lastp  = !curvsol.IsLastOpen();
4148           if (firstp) {
4149             ptf = curvsol.Value(first);
4150           }
4151           if (lastp) {
4152             ptl = curvsol.Value(last);
4153           }
4154           para = last;
4155           kount = 1;
4156           tgfound = Standard_False;
4157
4158           while (!tgfound) {
4159             para = (1.123*first + para)/2.123;
4160             tgfound = curvsol.D1u(para,ptvalid,tgvalid);
4161             if (!tgfound) {
4162               kount ++;
4163               tgfound = kount > 5;
4164             }
4165           }
4166           Handle(IntPatch_ALine) alig;
4167           if (kount <= 5) {
4168             Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
4169                                                  Quad1.Normale(ptvalid));
4170             if(qwe> 0.00000001) {
4171               trans1 = IntSurf_Out;
4172               trans2 = IntSurf_In;
4173             }
4174             else if(qwe<-0.00000001) {
4175               trans1 = IntSurf_In;
4176               trans2 = IntSurf_Out;
4177             }
4178             else {
4179               trans1=trans2=IntSurf_Undecided;
4180             }
4181             alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
4182             kept = Standard_True;
4183           }
4184           else {
4185             ptvalid = curvsol.Value(para);
4186             alig = new IntPatch_ALine(curvsol,Standard_False);
4187             kept = Standard_True;
4188             //-- std::cout << "Transition indeterminee" << std::endl;
4189           }
4190           if (kept) {
4191             Standard_Boolean Nfirstp = !firstp;
4192             Standard_Boolean Nlastp  = !lastp;
4193             ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
4194                           Nlastp,ptl,last,Multpoint,Tol);
4195             slin.Append(alig);
4196           }
4197         } // for (; aIt.More(); aIt.Next())
4198       } // for (i = 1; i <= NbSol; ++i)
4199     }
4200   }
4201     break;
4202
4203   default:
4204     return Standard_False;
4205
4206   } // switch (typint)
4207
4208   return Standard_True;
4209 }
4210 //=======================================================================
4211 //function : ExploreCurve
4212 //purpose  :
4213 //=======================================================================
4214 Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
4215                               const gp_Cone& aCo,
4216                               IntAna_Curve& aC,
4217                               const Standard_Real aTol,
4218                               IntAna_ListOfCurve& aLC)
4219
4220 {
4221   Standard_Boolean bFind=Standard_False;
4222   Standard_Real aTheta, aT1, aT2, aDst;
4223   gp_Pnt aPapx, aPx;
4224   //
4225   //aC.Dump();
4226   //
4227   aLC.Clear();
4228   aLC.Append(aC);
4229   //
4230   aPapx=aCo.Apex();
4231   //
4232   aC.Domain(aT1, aT2);
4233   //
4234   aPx=aC.Value(aT1);
4235   aDst=aPx.Distance(aPapx);
4236   if (aDst<aTol) {
4237     return bFind;
4238   }
4239   aPx=aC.Value(aT2);
4240   aDst=aPx.Distance(aPapx);
4241   if (aDst<aTol) {
4242     return bFind;
4243   }
4244   //
4245   bFind=aC.FindParameter(aPapx, aTheta);
4246   if (!bFind){
4247     return bFind;
4248   }
4249   //
4250   aPx=aC.Value(aTheta);
4251   aDst=aPx.Distance(aPapx);
4252   if (aDst>aTol) {
4253     return !bFind;
4254   }
4255   //
4256   // need to be splitted at aTheta
4257   IntAna_Curve aC1, aC2;
4258   //
4259   aC1=aC;
4260   aC1.SetDomain(aT1, aTheta);
4261   aC2=aC;
4262   aC2.SetDomain(aTheta, aT2);
4263   //
4264   aLC.Clear();
4265   aLC.Append(aC1);
4266   aLC.Append(aC2);
4267   //
4268   return bFind;
4269 }