src/IntPatch/IntPatch_ImpImpIntersection_4.gxx

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