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