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) &&
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   const ComputationMethods::stCoeffsValue &anEquationCoeffs = theBW.SICoeffs();
2696   const IntSurf_Quadric& aQuad1 = theBW.GetQSurface(1);
2697   const IntSurf_Quadric& aQuad2 = theBW.GetQSurface(2);
2698   const Standard_Boolean isReversed = theBW.IsReversed();
2699   const Standard_Real aTol2D = theBW.Get2dTolerance();
2700   const Standard_Real aTol3D = theBW.Get3dTolerance();
2701   const Standard_Real aPeriod = 2.0*M_PI;
2702   const Standard_Integer aNbMaxPoints = 2000;
2703   const Standard_Integer aNbMinPoints = 200;
2704   const Standard_Integer aNbPoints = Min(Max(aNbMinPoints,
2705                                   RealToInt(20.0*theCyl1.Radius())), aNbMaxPoints);
2706   const Standard_Real aStepMin = aTol2D,
2707                                   aStepMax = (aUSurf1l - aUSurf1f > M_PI / 100.0) ?
2708                                   (aUSurf1l - aUSurf1f) / IntToReal(aNbPoints) :
2709                                   aUSurf1l - aUSurf1f;
2710
2711   //The main idea of the algorithm is to change U1-parameter
2712   //(U-parameter of theCyl1) from aU1f to aU1l with some step
2713   //(step is adaptive) and to obtain set of intersection points.
2714
2715   for (Standard_Integer i = 0; i < theNbOfRanges; i++)
2716   {
2717     if (theRange[i].IsVoid())
2718       continue;
2719
2720     InscribeInterval(aUSurf1f, aUSurf1l, theRange[i], aTol2D, aPeriod);
2721   }
2722
2723   if (theRange[0].Union(theRange[1]))
2724   {
2725     // Works only if (theNbOfRanges == 2).
2726     theRange[1].SetVoid();
2727   }
2728
2729   //Critical points are the value of U1-parameter in the points
2730   //where WL must be decomposed.
2731
2732   //When U1 goes through critical points its value is set up to this
2733   //parameter forcefully and the intersection point is added in the line.
2734   //After that, the WL is broken (next U1 value will be correspond to the new WL).
2735
2736   //See CriticalPointsComputing(...) function to get detail information about this array.
2737   const Standard_Integer aNbCritPointsMax = 12;
2738   Standard_Real anU1crit[aNbCritPointsMax] = { Precision::Infinite(),
2739                                                Precision::Infinite(),
2740                                                Precision::Infinite(),
2741                                                Precision::Infinite(),
2742                                                Precision::Infinite(),
2743                                                Precision::Infinite(),
2744                                                Precision::Infinite(),
2745                                                Precision::Infinite(),
2746                                                Precision::Infinite(),
2747                                                Precision::Infinite(),
2748                                                Precision::Infinite(),
2749                                                Precision::Infinite() };
2750
2751   //This list of critical points is not full because it does not contain any points
2752   //which intersection line goes through V-bounds of cylinders in.
2753   //They are computed by numerical methods on - line (during algorithm working).
2754   //The moment is caught, when intersection line goes through V-bounds of any cylinder.
2755
2756   Standard_Integer aNbCritPoints = aNbCritPointsMax;
2757   CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
2758                           aPeriod, aTol2D, aNbCritPoints, anU1crit);
2759
2760   //Getting Walking-line
2761
2762   enum WLFStatus
2763   {
2764     // No points have been added in WL
2765     WLFStatus_Absent = 0,
2766     // WL contains at least one point
2767     WLFStatus_Exist = 1,
2768     // WL has been finished in some critical point
2769     // We should start new line
2770     WLFStatus_Broken = 2
2771   };
2772
2773   const Standard_Integer aNbWLines = 2;
2774   for (Standard_Integer aCurInterval = 0; aCurInterval < theNbOfRanges; aCurInterval++)
2775   {
2776     //Process every continuous region
2777     Standard_Boolean isAddedIntoWL[aNbWLines];
2778     for (Standard_Integer i = 0; i < aNbWLines; i++)
2779       isAddedIntoWL[i] = Standard_False;
2780
2781     Standard_Real anUf = 1.0, anUl = 0.0;
2782     if (!theRange[aCurInterval].GetBounds(anUf, anUl))
2783       continue;
2784
2785     const Standard_Boolean isDeltaPeriod = IsEqual(anUl - anUf, aPeriod);
2786
2787     //Inscribe and sort critical points
2788     for (Standard_Integer i = 0; i < aNbCritPoints; i++)
2789     {
2790       InscribePoint(anUf, anUl, anU1crit[i], 0.0, aPeriod, Standard_False);
2791     }
2792
2793     std::sort(anU1crit, anU1crit + aNbCritPoints);
2794
2795     while (anUf < anUl)
2796     {
2797       //Change value of U-parameter on the 1st surface from anUf to anUl
2798       //(anUf will be modified in the cycle body).
2799       //Step is computed adaptively (see comments below).
2800
2801       Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
2802       WLFStatus aWLFindStatus[aNbWLines];
2803       Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
2804       Standard_Real anUexpect[aNbWLines];
2805       Standard_Boolean isAddingWLEnabled[aNbWLines];
2806
2807       Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
2808       Handle(IntPatch_WLine) aWLine[aNbWLines];
2809       for (Standard_Integer i = 0; i < aNbWLines; i++)
2810       {
2811         aL2S[i] = new IntSurf_LineOn2S();
2812         aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
2813         aWLFindStatus[i] = WLFStatus_Absent;
2814         isAddingWLEnabled[i] = Standard_True;
2815         aU2[i] = aV1[i] = aV2[i] = 0.0;
2816         aV1Prev[i] = aV2Prev[i] = 0.0;
2817         anUexpect[i] = anUf;
2818       }
2819
2820       Standard_Real aCriticalDelta[aNbCritPointsMax] = { 0 };
2821       for (Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++)
2822       { //We are not interested in elements of aCriticalDelta array
2823         //if their index is greater than or equal to aNbCritPoints
2824
2825         aCriticalDelta[aCritPID] = anUf - anU1crit[aCritPID];
2826       }
2827
2828       Standard_Real anU1 = anUf, aMinCriticalParam = anUf;
2829       Standard_Boolean isFirst = Standard_True;
2830
2831       while (anU1 <= anUl)
2832       {
2833         //Change value of U-parameter on the 1st surface from anUf to anUl
2834         //(anUf will be modified in the cycle body). However, this cycle
2835         //can be broken if WL goes though some critical point.
2836         //Step is computed adaptively (see comments below).
2837
2838         for (Standard_Integer i = 0; i < aNbCritPoints; i++)
2839         {
2840           if ((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
2841           {
2842             //WL has gone through i-th critical point
2843             anU1 = anU1crit[i];
2844
2845             for (Standard_Integer j = 0; j < aNbWLines; j++)
2846             {
2847               aWLFindStatus[j] = WLFStatus_Broken;
2848               anUexpect[j] = anU1;
2849             }
2850
2851             break;
2852           }
2853         }
2854
2855         if (IsEqual(anU1, anUl))
2856         {
2857           for (Standard_Integer i = 0; i < aNbWLines; i++)
2858           {
2859             aWLFindStatus[i] = WLFStatus_Broken;
2860             anUexpect[i] = anU1;
2861
2862             if (isDeltaPeriod)
2863             {
2864               //if isAddedIntoWL[i] == TRUE WLine contains only one point
2865               //(which was end point of previous WLine). If we will
2866               //add point found on the current step WLine will contain only
2867               //two points. At that both these points will be equal to the
2868               //points found earlier. Therefore, new WLine will repeat
2869               //already existing WLine. Consequently, it is necessary
2870               //to forbid building new line in this case.
2871
2872               isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
2873             }
2874             else
2875             {
2876               isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
2877                                       (aWLFindStatus[i] == WLFStatus_Absent));
2878             }
2879           }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2880         }
2881         else
2882         {
2883           for (Standard_Integer i = 0; i < aNbWLines; i++)
2884           {
2885             isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
2886                                     (aWLFindStatus[i] == WLFStatus_Absent));
2887           }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2888         }
2889
2890         for (Standard_Integer i = 0; i < aNbWLines; i++)
2891         {
2892           const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
2893             aWLine[i]->Curve()->NbPoints();
2894
2895           if ((aWLFindStatus[i] == WLFStatus_Broken) ||
2896             (aWLFindStatus[i] == WLFStatus_Absent))
2897           {//Begin and end of WLine must be on boundary point
2898            //or on seam-edge strictly (if it is possible).
2899
2900             Standard_Real aTol = aTol2D;
2901             ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs,
2902                                                         aU2[i], &aTol);
2903             InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
2904
2905             aTol = Max(aTol, aTol2D);
2906
2907             if (Abs(aU2[i]) <= aTol)
2908               aU2[i] = 0.0;
2909             else if (Abs(aU2[i] - aPeriod) <= aTol)
2910               aU2[i] = aPeriod;
2911             else if (Abs(aU2[i] - aUSurf2f) <= aTol)
2912               aU2[i] = aUSurf2f;
2913             else if (Abs(aU2[i] - aUSurf2l) <= aTol)
2914               aU2[i] = aUSurf2l;
2915           }
2916           else
2917           {
2918             ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
2919             InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
2920           }
2921
2922           if (aNbPntsWL == 0)
2923           {//the line has not contained any points yet
2924             if (((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*aTol2D) &&
2925                 ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
2926                   (Abs(aU2[i] - aUSurf2l) < aTol2D)))
2927             {
2928               //In this case aU2[i] can have two values: current aU2[i] or
2929               //aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
2930               //correct value.
2931
2932               Standard_Boolean isIncreasing = Standard_True;
2933               ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs,
2934                                                       aPeriod, isIncreasing);
2935
2936               //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
2937               //then after the next step (when U1 will be increased) U2 will be
2938               //increased too. And we will go out of surface boundary.
2939               //Therefore, If U2(U1) is increasing then U2 must be equal aUSurf2f.
2940               //Analogically, if U2(U1) is decreasing.
2941               if (isIncreasing)
2942               {
2943                 aU2[i] = aUSurf2f;
2944               }
2945               else
2946               {
2947                 aU2[i] = aUSurf2l;
2948               }
2949             }
2950           }
2951           else
2952           {//aNbPntsWL > 0
2953             if (((aUSurf2l - aUSurf2f) >= aPeriod) &&
2954                 ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
2955                   (Abs(aU2[i] - aUSurf2l) < aTol2D)))
2956             {//end of the line
2957               Standard_Real aU2prev = 0.0, aV2prev = 0.0;
2958               if (isReversed)
2959                 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
2960               else
2961                 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
2962
2963               if (2.0*Abs(aU2prev - aU2[i]) > aPeriod)
2964               {
2965                 if (aU2prev > aU2[i])
2966                   aU2[i] += aPeriod;
2967                 else
2968                   aU2[i] -= aPeriod;
2969               }
2970             }
2971           }
2972
2973           ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
2974                                                               aV1[i], aV2[i]);
2975
2976           if (isFirst)
2977           {
2978             aV1Prev[i] = aV1[i];
2979             aV2Prev[i] = aV2[i];
2980           }
2981         }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2982
2983         isFirst = Standard_False;
2984
2985         //Looking for points into WLine
2986         Standard_Boolean isBroken = Standard_False;
2987         for (Standard_Integer i = 0; i < aNbWLines; i++)
2988         {
2989           if (!isAddingWLEnabled[i])
2990           {
2991             Standard_Boolean isBoundIntersect = Standard_False;
2992             if ((Abs(aV1[i] - aVSurf1f) <= aTol2D) ||
2993                 ((aV1[i] - aVSurf1f)*(aV1Prev[i] - aVSurf1f) < 0.0))
2994             {
2995               isBoundIntersect = Standard_True;
2996             }
2997             else if ((Abs(aV1[i] - aVSurf1l) <= aTol2D) ||
2998                     ((aV1[i] - aVSurf1l)*(aV1Prev[i] - aVSurf1l) < 0.0))
2999             {
3000               isBoundIntersect = Standard_True;
3001             }
3002             else if ((Abs(aV2[i] - aVSurf2f) <= aTol2D) ||
3003                     ((aV2[i] - aVSurf2f)*(aV2Prev[i] - aVSurf2f) < 0.0))
3004             {
3005               isBoundIntersect = Standard_True;
3006             }
3007             else if ((Abs(aV2[i] - aVSurf2l) <= aTol2D) ||
3008                     ((aV2[i] - aVSurf2l)*(aV2Prev[i] - aVSurf2l) < 0.0))
3009             {
3010               isBoundIntersect = Standard_True;
3011             }
3012
3013             if (aWLFindStatus[i] == WLFStatus_Broken)
3014               isBroken = Standard_True;
3015
3016             if (!isBoundIntersect)
3017             {
3018               continue;
3019             }
3020             else
3021             {
3022               anUexpect[i] = anU1;
3023             }
3024           }
3025
3026           // True if the current point already in the domain
3027           const Standard_Boolean isInscribe =
3028               ((aUSurf2f - aU2[i]) <= aTol2D) && ((aU2[i] - aUSurf2l) <= aTol2D) &&
3029               ((aVSurf1f - aV1[i]) <= aTol2D) && ((aV1[i] - aVSurf1l) <= aTol2D) &&
3030               ((aVSurf2f - aV2[i]) <= aTol2D) && ((aV2[i] - aVSurf2l) <= aTol2D);
3031
3032           //isVIntersect == TRUE if intersection line intersects two (!)
3033           //V-bounds of cylinder (1st or 2nd - no matter)
3034           const Standard_Boolean isVIntersect =
3035               (((aVSurf1f - aV1[i])*(aVSurf1f - aV1Prev[i]) < RealSmall()) &&
3036                 ((aVSurf1l - aV1[i])*(aVSurf1l - aV1Prev[i]) < RealSmall())) ||
3037               (((aVSurf2f - aV2[i])*(aVSurf2f - aV2Prev[i]) < RealSmall()) &&
3038                 ((aVSurf2l - aV2[i])*(aVSurf2l - aV2Prev[i]) < RealSmall()));
3039
3040           //isFound1 == TRUE if intersection line intersects V-bounds
3041           //  (First or Last - no matter) of the 1st cylynder
3042           //isFound2 == TRUE if intersection line intersects V-bounds
3043           //  (First or Last - no matter) of the 2nd cylynder
3044           Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
3045           Standard_Boolean isForce = Standard_False;
3046
3047           if (aWLFindStatus[i] == WLFStatus_Absent)
3048           {
3049             if (((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1 - aUSurf1l) < aTol2D))
3050             {
3051               isForce = Standard_True;
3052             }
3053           }
3054
3055           theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
3056                                  aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i, isForce,
3057                                  isFound1, isFound2);
3058
3059           const Standard_Boolean isPrevVBound = !isVIntersect &&
3060                                                 ((Abs(aV1Prev[i] - aVSurf1f) <= aTol2D) ||
3061                                                  (Abs(aV1Prev[i] - aVSurf1l) <= aTol2D) ||
3062                                                  (Abs(aV2Prev[i] - aVSurf2f) <= aTol2D) ||
3063                                                  (Abs(aV2Prev[i] - aVSurf2l) <= aTol2D));
3064
3065           aV1Prev[i] = aV1[i];
3066           aV2Prev[i] = aV2[i];
3067
3068           if ((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
3069           {
3070             aWLFindStatus[i] = WLFStatus_Broken; //start a new line
3071           }
3072           else if (isInscribe)
3073           {
3074             if ((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
3075             {
3076               aWLFindStatus[i] = WLFStatus_Exist;
3077             }
3078
3079             if ((aWLFindStatus[i] != WLFStatus_Broken) ||
3080               (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
3081             {
3082               if (aWLine[i]->NbPnts() > 0)
3083               {
3084                 Standard_Real aU2p = 0.0, aV2p = 0.0;
3085                 if (isReversed)
3086                   aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
3087                 else
3088                   aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
3089
3090                 const Standard_Real aDelta = aU2[i] - aU2p;
3091
3092                 if (2.0 * Abs(aDelta) > aPeriod)
3093                 {
3094                   if (aDelta > 0.0)
3095                   {
3096                     aU2[i] -= aPeriod;
3097                   }
3098                   else
3099                   {
3100                     aU2[i] += aPeriod;
3101                   }
3102                 }
3103               }
3104
3105               if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
3106                                 gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
3107                                 aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3108                                 aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
3109                                 aWLine[i]->Curve(), i, aTol3D, aTol2D, isForce))
3110               {
3111                 if (aWLFindStatus[i] == WLFStatus_Absent)
3112                 {
3113                   aWLFindStatus[i] = WLFStatus_Exist;
3114                 }
3115               }
3116               else if (!isFound1 && !isFound2)
3117               {//We do not add any point while doing this iteration
3118                 if (aWLFindStatus[i] == WLFStatus_Exist)
3119                 {
3120                   aWLFindStatus[i] = WLFStatus_Broken;
3121                 }
3122               }
3123             }
3124           }
3125           else
3126           {//We do not add any point while doing this iteration
3127             if (aWLFindStatus[i] == WLFStatus_Exist)
3128             {
3129               aWLFindStatus[i] = WLFStatus_Broken;
3130             }
3131           }
3132
3133           if (aWLFindStatus[i] == WLFStatus_Broken)
3134             isBroken = Standard_True;
3135         }//for(Standard_Integer i = 0; i < aNbWLines; i++)
3136
3137         if (isBroken)
3138         {//current lines are filled. Go to the next lines
3139           anUf = anU1;
3140
3141           Standard_Boolean isAdded = Standard_True;
3142
3143           for (Standard_Integer i = 0; i < aNbWLines; i++)
3144           {
3145             if (isAddingWLEnabled[i])
3146             {
3147               continue;
3148             }
3149
3150             isAdded = Standard_False;
3151
3152             Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
3153
3154             theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
3155                                    aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i,
3156                                    Standard_False, isFound1, isFound2);
3157
3158             if (isFound1 || isFound2)
3159             {
3160               isAdded = Standard_True;
3161             }
3162
3163             if (aWLine[i]->NbPnts() > 0)
3164             {
3165               Standard_Real aU2p = 0.0, aV2p = 0.0;
3166               if (isReversed)
3167                 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
3168               else
3169                 aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
3170
3171               const Standard_Real aDelta = aU2[i] - aU2p;
3172
3173               if (2 * Abs(aDelta) > aPeriod)
3174               {
3175                 if (aDelta > 0.0)
3176                 {
3177                   aU2[i] -= aPeriod;
3178                 }
3179                 else
3180                 {
3181                   aU2[i] += aPeriod;
3182                 }
3183               }
3184             }
3185
3186             if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
3187                               Standard_True, gp_Pnt2d(anU1, aV1[i]),
3188                               gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3189                               aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3190                               aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
3191                               i, aTol3D, aTol2D, Standard_False))
3192             {
3193               isAdded = Standard_True;
3194             }
3195           }
3196
3197           if (!isAdded)
3198           {
3199             //Before breaking WL, we must complete it correctly
3200             //(e.g. to prolong to the surface boundary).
3201             //Therefore, we take the point last added in some WL
3202             //(have maximal U1-parameter) and try to add it in
3203             //the current WL.
3204             Standard_Real anUmaxAdded = RealFirst();
3205
3206             {
3207               Standard_Boolean isChanged = Standard_False;
3208               for (Standard_Integer i = 0; i < aNbWLines; i++)
3209               {
3210                 if ((aWLFindStatus[i] == WLFStatus_Absent) || (aWLine[i]->NbPnts() == 0))
3211                   continue;
3212
3213                 Standard_Real aU1c = 0.0, aV1c = 0.0;
3214                 if (isReversed)
3215                   aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
3216                 else
3217                   aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
3218
3219                 anUmaxAdded = Max(anUmaxAdded, aU1c);
3220                 isChanged = Standard_True;
3221               }
3222
3223               if (!isChanged)
3224               { //If anUmaxAdded were not changed in previous cycle then
3225                 //we would break existing WLines.
3226                 break;
3227               }
3228             }
3229
3230             for (Standard_Integer i = 0; i < aNbWLines; i++)
3231             {
3232               if (isAddingWLEnabled[i])
3233               {
3234                 continue;
3235               }
3236
3237               ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
3238                 aU2[i], aV1[i], aV2[i]);
3239
3240               AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
3241                              Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
3242                              gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
3243                              aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
3244                              aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
3245                              i, aTol3D, aTol2D, Standard_False);
3246             }
3247           }
3248
3249           break;
3250         }
3251
3252         //Step computing
3253
3254         {
3255           //Step of aU1-parameter is computed adaptively. The algorithm
3256           //aims to provide given aDeltaV1 and aDeltaV2 values (if it is
3257           //possible because the intersection line can go along V-isoline)
3258           //in every iteration. It allows avoiding "flying" intersection
3259           //points too far each from other (see issue #24915).
3260
3261           const Standard_Real aDeltaV1 = aRangeS1.Delta() / IntToReal(aNbPoints),
3262                               aDeltaV2 = aRangeS2.Delta() / IntToReal(aNbPoints);
3263
3264           math_Matrix aMatr(1, 3, 1, 5);
3265
3266           Standard_Real aMinUexp = RealLast();
3267
3268           for (Standard_Integer i = 0; i < aNbWLines; i++)
3269           {
3270             if (aTol2D < (anUexpect[i] - anU1))
3271             {
3272               continue;
3273             }
3274
3275             if (aWLFindStatus[i] == WLFStatus_Absent)
3276             {
3277               anUexpect[i] += aStepMax;
3278               aMinUexp = Min(aMinUexp, anUexpect[i]);
3279               continue;
3280             }
3281
3282             Standard_Real aStepTmp = aStepMax;
3283
3284             const Standard_Real aSinU1 = sin(anU1),
3285                                 aCosU1 = cos(anU1),
3286                                 aSinU2 = sin(aU2[i]),
3287                                 aCosU2 = cos(aU2[i]);
3288
3289             aMatr.SetCol(1, anEquationCoeffs.mVecC1);
3290             aMatr.SetCol(2, anEquationCoeffs.mVecC2);
3291             aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
3292             aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
3293             aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
3294                             anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
3295                             anEquationCoeffs.mVecD);
3296
3297             //The main idea is in solving of linearized system (2)
3298             //(see description to ComputationMethods class) in order to find new U1-value
3299             //to provide new value V1 or V2, which differs from current one by aDeltaV1 or
3300             //aDeltaV2 respectively.
3301
3302             //While linearizing, following Taylor formulas are used:
3303             //    cos(x0+dx) = cos(x0) - sin(x0)*dx
3304             //    sin(x0+dx) = sin(x0) + cos(x0)*dx
3305
3306             //Consequently cos(U1), cos(U2), sin(U1) and sin(U2) in the system (2)
3307             //must be substituted by corresponding values.
3308
3309             //ATTENTION!!!
3310             //The solution is approximate. More over, all requirements to the
3311             //linearization must be satisfied in order to obtain quality result.
3312
3313             if (!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aStepTmp))
3314             {
3315               //To avoid cycling-up
3316               anUexpect[i] += aStepMax;
3317               aMinUexp = Min(aMinUexp, anUexpect[i]);
3318
3319               continue;
3320             }
3321
3322             if (aStepTmp < aStepMin)
3323               aStepTmp = aStepMin;
3324
3325             if (aStepTmp > aStepMax)
3326               aStepTmp = aStepMax;
3327
3328             anUexpect[i] = anU1 + aStepTmp;
3329             aMinUexp = Min(aMinUexp, anUexpect[i]);
3330           }
3331
3332           anU1 = aMinUexp;
3333         }
3334
3335         if (Precision::PConfusion() >= (anUl - anU1))
3336           anU1 = anUl;
3337
3338         anUf = anU1;
3339
3340         for (Standard_Integer i = 0; i < aNbWLines; i++)
3341         {
3342           if (aWLine[i]->NbPnts() != 1)
3343             isAddedIntoWL[i] = Standard_False;
3344
3345           if (anU1 == anUl)
3346           {//strictly equal. Tolerance is considered above.
3347             anUexpect[i] = anUl;
3348           }
3349         }
3350       }
3351
3352       for (Standard_Integer i = 0; i < aNbWLines; i++)
3353       {
3354         if ((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
3355         {
3356           isTheEmpty = Standard_False;
3357           Standard_Real u1, v1, u2, v2;
3358           aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
3359           IntPatch_Point aP;
3360           aP.SetParameter(u1);
3361           aP.SetParameters(u1, v1, u2, v2);
3362           aP.SetTolerance(aTol3D);
3363           aP.SetValue(aWLine[i]->Point(1).Value());
3364
3365           theSPnt.Append(aP);
3366         }
3367         else if (aWLine[i]->NbPnts() > 1)
3368         {
3369           Standard_Boolean isGood = Standard_True;
3370
3371           if (aWLine[i]->NbPnts() == 2)
3372           {
3373             const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
3374             const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
3375
3376             if (aPf.IsSame(aPl, Precision::Confusion()))
3377               isGood = Standard_False;
3378           }
3379
3380           if (isGood)
3381           {
3382             isTheEmpty = Standard_False;
3383             isAddedIntoWL[i] = Standard_True;
3384             SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
3385                                  anEquationCoeffs, i, aNbPoints, 1,
3386               aWLine[i]->NbPnts(), aTol2D, aPeriod,
3387               isReversed);
3388
3389             aWLine[i]->ComputeVertexParameters(aTol3D);
3390             theSlin.Append(aWLine[i]);
3391           }
3392         }
3393         else
3394         {
3395           isAddedIntoWL[i] = Standard_False;
3396         }
3397
3398 #ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
3399         //aWLine[i]->Dump();
3400 #endif
3401       }
3402     }
3403   }
3404
3405
3406   //Delete the points in theSPnt, which
3407   //lie at least in one of the line in theSlin.
3408   for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
3409   {
3410     for (Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
3411     {
3412       Handle(IntPatch_WLine) aWLine1(Handle(IntPatch_WLine)::
3413         DownCast(theSlin.Value(aNbLin)));
3414
3415       const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
3416       const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
3417
3418       const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
3419       if (aPntCur.IsSame(aPntFWL1, Precision::Confusion()) ||
3420         aPntCur.IsSame(aPntLWL1, Precision::Confusion()))
3421       {
3422         theSPnt.Remove(aNbPnt);
3423         aNbPnt--;
3424         break;
3425       }
3426     }
3427   }
3428
3429   //Try to add new points in the neighborhood of existing point
3430   for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
3431   {
3432     // Standard algorithm (implemented above) could not find any
3433     // continuous curve in neighborhood of aPnt2S (e.g. because
3434     // this curve is too small; see tests\bugs\modalg_5\bug25292_35 and _36).
3435     // Here, we will try to find several new points nearer to aPnt2S.
3436
3437     // The algorithm below tries to find two points in every
3438     // intervals [u1 - aStepMax, u1] and [u1, u1 + aStepMax]
3439     // and every new point will be in maximal distance from
3440     // u1. If these two points exist they will be joined
3441     // by the intersection curve.
3442
3443     const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
3444
3445     Standard_Real u1, v1, u2, v2;
3446     aPnt2S.Parameters(u1, v1, u2, v2);
3447
3448     Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
3449     Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
3450
3451     //Define the index of WLine, which lies the point aPnt2S in.
3452     Standard_Integer anIndex = 0;
3453
3454     Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
3455     if (isReversed)
3456     {
3457       anUf = Max(u2 - aStepMax, aUSurf1f);
3458       anUl = Min(u2 + aStepMax, aUSurf1l);
3459       aCurU2 = u1;
3460     }
3461     else
3462     {
3463       anUf = Max(u1 - aStepMax, aUSurf1f);
3464       anUl = Min(u1 + aStepMax, aUSurf1l);
3465       aCurU2 = u2;
3466     }
3467
3468     const Standard_Real anUinf = anUf, anUsup = anUl, anUmid = 0.5*(anUf + anUl);
3469
3470     {
3471       //Find the value of anIndex variable.
3472       Standard_Real aDelta = RealLast();
3473       for (Standard_Integer i = 0; i < aNbWLines; i++)
3474       {
3475         Standard_Real anU2t = 0.0;
3476         if (!ComputationMethods::CylCylComputeParameters(anUmid, i, anEquationCoeffs, anU2t))
3477           continue;
3478
3479         Standard_Real aDU2 = fmod(Abs(anU2t - aCurU2), aPeriod);
3480         aDU2 = Min(aDU2, Abs(aDU2 - aPeriod));
3481         if (aDU2 < aDelta)
3482         {
3483           aDelta = aDU2;
3484           anIndex = i;
3485         }
3486       }
3487     }
3488
3489     // Bisection method is used in order to find every new point.
3490     // I.e. if we need to find intersection point in the interval [anUinf, anUmid]
3491     // we check the point anUC = 0.5*(anUinf+anUmid). If it is an intersection point
3492     // we try to find another point in the interval [anUinf, anUC] (because we find the point in
3493     // maximal distance from anUmid). If it is not then we try to find another point in the
3494     // interval [anUC, anUmid]. Next iterations will be made analogically.
3495     // When we find intersection point in the interval [anUmid, anUsup] we try to find
3496     // another point in the interval [anUC, anUsup] if anUC is intersection point and
3497     // in the interval [anUmid, anUC], otherwise.
3498
3499     Standard_Real anAddedPar[2] = { anUmid, anUmid };
3500
3501     for (Standard_Integer aParID = 0; aParID < 2; aParID++)
3502     {
3503       if (aParID == 0)
3504       {
3505         anUf = anUinf;
3506         anUl = anUmid;
3507       }
3508       else // if(aParID == 1)
3509       {
3510         anUf = anUmid;
3511         anUl = anUsup;
3512       }
3513
3514       Standard_Real &aPar1 = (aParID == 0) ? anUf : anUl,
3515                     &aPar2 = (aParID == 0) ? anUl : anUf;
3516
3517       while (Abs(aPar2 - aPar1) > aStepMin)
3518       {
3519         Standard_Real anUC = 0.5*(anUf + anUl);
3520         Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
3521         Standard_Boolean isDone = ComputationMethods::
3522                 CylCylComputeParameters(anUC, anIndex, anEquationCoeffs, aU2, aV1, aV2);
3523
3524         if (isDone)
3525         {
3526           if (Abs(aV1 - aVSurf1f) <= aTol2D)
3527             aV1 = aVSurf1f;
3528
3529           if (Abs(aV1 - aVSurf1l) <= aTol2D)
3530             aV1 = aVSurf1l;
3531
3532           if (Abs(aV2 - aVSurf2f) <= aTol2D)
3533             aV2 = aVSurf2f;
3534
3535           if (Abs(aV2 - aVSurf2l) <= aTol2D)
3536             aV2 = aVSurf2l;
3537
3538           isDone = AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
3539                                   Standard_True, gp_Pnt2d(anUC, aV1), gp_Pnt2d(aU2, aV2),
3540                                   aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3541                                   aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
3542                                   aPeriod, aWLine->Curve(), anIndex, aTol3D,
3543                                   aTol2D, Standard_False, Standard_True);
3544         }
3545
3546         if (isDone)
3547         {
3548           anAddedPar[0] = Min(anAddedPar[0], anUC);
3549           anAddedPar[1] = Max(anAddedPar[1], anUC);
3550           aPar2 = anUC;
3551         }
3552         else
3553         {
3554           aPar1 = anUC;
3555         }
3556       }
3557     }
3558
3559     //Fill aWLine by additional points
3560     if (anAddedPar[1] - anAddedPar[0] > aStepMin)
3561     {
3562       for (Standard_Integer aParID = 0; aParID < 2; aParID++)
3563       {
3564         Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
3565         ComputationMethods::CylCylComputeParameters(anAddedPar[aParID], anIndex,
3566                                                   anEquationCoeffs, aU2, aV1, aV2);
3567
3568         AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
3569                         gp_Pnt2d(anAddedPar[aParID], aV1), gp_Pnt2d(aU2, aV2),
3570                         aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
3571                         aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod, aWLine->Curve(),
3572                         anIndex, aTol3D, aTol2D, Standard_False, Standard_False);
3573       }
3574
3575       SeekAdditionalPoints(aQuad1, aQuad2, aWLine->Curve(),
3576                             anEquationCoeffs, anIndex, aNbMinPoints,
3577                             1, aWLine->NbPnts(), aTol2D, aPeriod,
3578                             isReversed);
3579
3580       aWLine->ComputeVertexParameters(aTol3D);
3581       theSlin.Append(aWLine);
3582
3583       theSPnt.Remove(aNbPnt);
3584       aNbPnt--;
3585     }
3586   }
3587
3588   return IntPatch_ImpImpIntersection::IntStatus_OK;
3589 }
3590
3591 //=======================================================================
3592 //function : IntCyCy
3593 //purpose  :
3594 //=======================================================================
3595 IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1,
3596                                                const IntSurf_Quadric& theQuad2,
3597                                                const Standard_Real theTol3D,
3598                                                const Standard_Real theTol2D,
3599                                                const Bnd_Box2d& theUVSurf1,
3600                                                const Bnd_Box2d& theUVSurf2,
3601                                                Standard_Boolean& isTheEmpty,
3602                                                Standard_Boolean& isTheSameSurface,
3603                                                Standard_Boolean& isTheMultiplePoint,
3604                                                IntPatch_SequenceOfLine& theSlin,
3605                                                IntPatch_SequenceOfPoint& theSPnt)
3606 {
3607   isTheEmpty = Standard_True;
3608   isTheSameSurface = Standard_False;
3609   isTheMultiplePoint = Standard_False;
3610   theSlin.Clear();
3611   theSPnt.Clear();
3612
3613   const gp_Cylinder aCyl1 = theQuad1.Cylinder(),
3614                     aCyl2 = theQuad2.Cylinder();
3615
3616   IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
3617
3618   if (!anInter.IsDone())
3619   {
3620     return IntPatch_ImpImpIntersection::IntStatus_Fail;
3621   }
3622
3623   if(anInter.TypeInter() != IntAna_NoGeometricSolution)
3624   {
3625     if (CyCyAnalyticalIntersect(theQuad1, theQuad2, anInter,
3626                                 theTol3D, isTheEmpty,
3627                                 isTheSameSurface, isTheMultiplePoint,
3628                                 theSlin, theSPnt))
3629     {
3630       return IntPatch_ImpImpIntersection::IntStatus_OK;
3631     }
3632   }
3633
3634   //Here, intersection line is not an analytical curve(line, circle, ellipsis etc.)
3635
3636   Standard_Real aUSBou[2][2], aVSBou[2][2]; //const
3637
3638   theUVSurf1.Get(aUSBou[0][0], aVSBou[0][0], aUSBou[0][1], aVSBou[0][1]);
3639   theUVSurf2.Get(aUSBou[1][0], aVSBou[1][0], aUSBou[1][1], aVSBou[1][1]);
3640
3641   const Standard_Real aPeriod = 2.0*M_PI;
3642   const Standard_Integer aNbWLines = 2;
3643
3644   const ComputationMethods::stCoeffsValue anEquationCoeffs1(aCyl1, aCyl2);
3645   const ComputationMethods::stCoeffsValue anEquationCoeffs2(aCyl2, aCyl1);
3646
3647   //Boundaries.
3648   //Intersection result can include two non-connected regions
3649   //(see WorkWithBoundaries::BoundariesComputing(...) method).
3650   const Standard_Integer aNbOfBoundaries = 2;
3651   Bnd_Range anURange[2][aNbOfBoundaries];   //const
3652
3653   if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs1, aPeriod, anURange[0]))
3654     return IntPatch_ImpImpIntersection::IntStatus_OK;
3655
3656   if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs2, aPeriod, anURange[1]))
3657     return IntPatch_ImpImpIntersection::IntStatus_OK;
3658
3659   //anURange[*] can be in different periodic regions in
3660   //compare with First-Last surface. E.g. the surface
3661   //is full cylinder [0, 2*PI] but anURange is [5, 7].
3662   //Trivial common range computation returns [5, 2*PI] and
3663   //its summary length is 2*PI-5 == 1.28... only. That is wrong.
3664   //This problem can be solved by the following
3665   //algorithm:
3666   // 1. split anURange[*] by the surface boundary;
3667   // 2. shift every new range in order to inscribe it
3668   //      in [Ufirst, Ulast] of cylinder;
3669   // 3. consider only common ranges between [Ufirst, Ulast]
3670   //    and new ranges.
3671   //
3672   // In above example, we obtain following:
3673   // 1. two ranges: [5, 2*PI] and [2*PI, 7];
3674   // 2. after shifting: [5, 2*PI] and [0, 7-2*PI];
3675   // 3. Common ranges: ([5, 2*PI] and [0, 2*PI]) == [5, 2*PI],
3676   //                   ([0, 7-2*PI] and [0, 2*PI]) == [0, 7-2*PI];
3677   // 4. Their summary length is (2*PI-5)+(7-2*PI-0)==7-5==2 ==> GOOD.
3678
3679   Standard_Real aSumRange[2] = { 0.0, 0.0 };
3680   Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
3681   for (Standard_Integer aCID = 0; aCID < 2; aCID++)
3682   {
3683     anAlloc->Reset();
3684     NCollection_List<Bnd_Range> aListOfRng(anAlloc);
3685
3686     aListOfRng.Append(anURange[aCID][0]);
3687     aListOfRng.Append(anURange[aCID][1]);
3688
3689     const Standard_Real aSplitArr[3] = {aUSBou[aCID][0], aUSBou[aCID][1], 0.0};
3690
3691     NCollection_List<Bnd_Range>::Iterator anITrRng;
3692     for (Standard_Integer aSInd = 0; aSInd < 3; aSInd++)
3693     {
3694       NCollection_List<Bnd_Range> aLstTemp(aListOfRng);
3695       aListOfRng.Clear();
3696       for (anITrRng.Init(aLstTemp); anITrRng.More(); anITrRng.Next())
3697       {
3698         Bnd_Range& aRng = anITrRng.Value();
3699         aRng.Split(aSplitArr[aSInd], aListOfRng, aPeriod);
3700       }
3701     }
3702
3703     anITrRng.Init(aListOfRng);
3704     for (; anITrRng.More(); anITrRng.Next())
3705     {
3706       Bnd_Range& aCurrRange = anITrRng.Value();
3707
3708       Bnd_Range aBoundR;
3709       aBoundR.Add(aUSBou[aCID][0]);
3710       aBoundR.Add(aUSBou[aCID][1]);
3711
3712       if (!InscribeInterval(aUSBou[aCID][0], aUSBou[aCID][1],
3713                                           aCurrRange, theTol2D, aPeriod))
3714       {
3715         //If aCurrRange does not have common block with
3716         //[Ufirst, Ulast] of cylinder then we will try
3717         //to inscribe [Ufirst, Ulast] in the boundaries of aCurrRange.
3718         Standard_Real aF = 1.0, aL = 0.0;
3719         if (!aCurrRange.GetBounds(aF, aL))
3720           continue;
3721
3722         if ((aL < aUSBou[aCID][0]))
3723         {
3724           aCurrRange.Shift(aPeriod);
3725         }
3726         else if (aF > aUSBou[aCID][1])
3727         {
3728           aCurrRange.Shift(-aPeriod);
3729         }
3730       }
3731
3732       aBoundR.Common(aCurrRange);
3733
3734       const Standard_Real aDelta = aBoundR.Delta();
3735
3736       if (aDelta > 0.0)
3737       {
3738         aSumRange[aCID] += aDelta;
3739       }
3740     }
3741   }
3742
3743   //The bigger range the bigger number of points in Walking-line (WLine)
3744   //we will be able to add and consequently we will obtain more
3745   //precise intersection line.
3746   //Every point of WLine is determined as function from U1-parameter,
3747   //where U1 is U-parameter on 1st quadric.
3748   //Therefore, we should use quadric with bigger range as 1st parameter
3749   //in IntCyCy() function.
3750   //On the other hand, there is no point in reversing in case of
3751   //analytical intersection (when result is line, ellipse, point...).
3752   //This result is independent of the arguments order.
3753   const Standard_Boolean isToReverse = (aSumRange[1] > aSumRange[0]);
3754
3755   if (isToReverse)
3756   {
3757     const WorkWithBoundaries aBoundWork(theQuad2, theQuad1, anEquationCoeffs2,
3758                                         theUVSurf2, theUVSurf1, aNbWLines,
3759                                         aPeriod, theTol3D, theTol2D, Standard_True);
3760
3761     return CyCyNoGeometric(aCyl2, aCyl1, aBoundWork, anURange[1], aNbOfBoundaries,
3762                               isTheEmpty, theSlin, theSPnt);
3763   }
3764   else
3765   {
3766     const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs1,
3767                                         theUVSurf1, theUVSurf2, aNbWLines,
3768                                         aPeriod, theTol3D, theTol2D, Standard_False);
3769
3770     return CyCyNoGeometric(aCyl1, aCyl2, aBoundWork, anURange[0], aNbOfBoundaries,
3771                               isTheEmpty, theSlin, theSPnt);
3772   }
3773 }
3774
3775 //=======================================================================
3776 //function : IntCySp
3777 //purpose  :
3778 //=======================================================================
3779 Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
3780                          const IntSurf_Quadric& Quad2,
3781                          const Standard_Real Tol,
3782                          const Standard_Boolean Reversed,
3783                          Standard_Boolean& Empty,
3784                          Standard_Boolean& Multpoint,
3785                          IntPatch_SequenceOfLine& slin,
3786                          IntPatch_SequenceOfPoint& spnt)
3787
3788 {
3789   Standard_Integer i;
3790
3791   IntSurf_TypeTrans trans1,trans2;
3792   IntAna_ResultType typint;
3793   IntPatch_Point ptsol;
3794   gp_Circ cirsol;
3795
3796   gp_Cylinder Cy;
3797   gp_Sphere Sp;
3798
3799   if (!Reversed) {
3800     Cy = Quad1.Cylinder();
3801     Sp = Quad2.Sphere();
3802   }
3803   else {
3804     Cy = Quad2.Cylinder();
3805     Sp = Quad1.Sphere();
3806   }
3807   IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
3808
3809   if (!inter.IsDone()) {return Standard_False;}
3810
3811   typint = inter.TypeInter();
3812   Standard_Integer NbSol = inter.NbSolutions();
3813   Empty = Standard_False;
3814
3815   switch (typint) {
3816
3817   case IntAna_Empty :
3818     {
3819       Empty = Standard_True;
3820     }
3821     break;
3822
3823   case IntAna_Point :
3824     {
3825       gp_Pnt psol(inter.Point(1));
3826       Standard_Real U1,V1,U2,V2;
3827       Quad1.Parameters(psol,U1,V1);
3828       Quad2.Parameters(psol,U2,V2);
3829       ptsol.SetValue(psol,Tol,Standard_True);
3830       ptsol.SetParameters(U1,V1,U2,V2);
3831       spnt.Append(ptsol);
3832     }
3833     break;
3834
3835   case IntAna_Circle:
3836     {
3837       cirsol = inter.Circle(1);
3838       gp_Vec Tgt;
3839       gp_Pnt ptref;
3840       ElCLib::D1(0.,cirsol,ptref,Tgt);
3841
3842       if (NbSol == 1) {
3843         gp_Vec TestCurvature(ptref,Sp.Location());
3844         gp_Vec Normsp,Normcyl;
3845         if (!Reversed) {
3846           Normcyl = Quad1.Normale(ptref);
3847           Normsp  = Quad2.Normale(ptref);
3848         }
3849         else {
3850           Normcyl = Quad2.Normale(ptref);
3851           Normsp  = Quad1.Normale(ptref);
3852         }
3853
3854         IntSurf_Situation situcyl;
3855         IntSurf_Situation situsp;
3856
3857         if (Normcyl.Dot(TestCurvature) > 0.) {
3858           situsp = IntSurf_Outside;
3859           if (Normsp.Dot(Normcyl) > 0.) {
3860             situcyl = IntSurf_Inside;
3861           }
3862           else {
3863             situcyl = IntSurf_Outside;
3864           }
3865         }
3866         else {
3867           situsp = IntSurf_Inside;
3868           if (Normsp.Dot(Normcyl) > 0.) {
3869             situcyl = IntSurf_Outside;
3870           }
3871           else {
3872             situcyl = IntSurf_Inside;
3873           }
3874         }
3875         Handle(IntPatch_GLine) glig;
3876         if (!Reversed) {
3877           glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
3878         }
3879         else {
3880           glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
3881         }
3882         slin.Append(glig);
3883       }
3884       else {
3885         if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
3886           trans1 = IntSurf_Out;
3887           trans2 = IntSurf_In;
3888         }
3889         else {
3890           trans1 = IntSurf_In;
3891           trans2 = IntSurf_Out;
3892         }
3893         Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3894         slin.Append(glig);
3895
3896         cirsol = inter.Circle(2);
3897         ElCLib::D1(0.,cirsol,ptref,Tgt);
3898         Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
3899         if(qwe> 0.0000001) {
3900           trans1 = IntSurf_Out;
3901           trans2 = IntSurf_In;
3902         }
3903         else if(qwe<-0.0000001) {
3904           trans1 = IntSurf_In;
3905           trans2 = IntSurf_Out;
3906         }
3907         else {
3908           trans1=trans2=IntSurf_Undecided;
3909         }
3910         glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
3911         slin.Append(glig);
3912       }
3913     }
3914     break;
3915
3916   case IntAna_NoGeometricSolution:
3917     {
3918       gp_Pnt psol;
3919       Standard_Real U1,V1,U2,V2;
3920       IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
3921       if (!anaint.IsDone()) {
3922         return Standard_False;
3923       }
3924
3925       if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
3926         Empty = Standard_True;
3927       }
3928       else {
3929
3930         NbSol = anaint.NbPnt();
3931         for (i = 1; i <= NbSol; i++) {
3932           psol = anaint.Point(i);
3933           Quad1.Parameters(psol,U1,V1);
3934           Quad2.Parameters(psol,U2,V2);
3935           ptsol.SetValue(psol,Tol,Standard_True);
3936           ptsol.SetParameters(U1,V1,U2,V2);
3937           spnt.Append(ptsol);
3938         }
3939
3940         gp_Pnt ptvalid,ptf,ptl;
3941         gp_Vec tgvalid;
3942         Standard_Real first,last,para;
3943         IntAna_Curve curvsol;
3944         Standard_Boolean tgfound;
3945         Standard_Integer kount;
3946
3947         NbSol = anaint.NbCurve();
3948         for (i = 1; i <= NbSol; i++) {
3949           curvsol = anaint.Curve(i);
3950           curvsol.Domain(first,last);
3951           ptf = curvsol.Value(first);
3952           ptl = curvsol.Value(last);
3953
3954           para = last;
3955           kount = 1;
3956           tgfound = Standard_False;
3957
3958           while (!tgfound) {
3959             para = (1.123*first + para)/2.123;
3960             tgfound = curvsol.D1u(para,ptvalid,tgvalid);
3961             if (!tgfound) {
3962               kount ++;
3963               tgfound = kount > 5;
3964             }
3965           }
3966           Handle(IntPatch_ALine) alig;
3967           if (kount <= 5) {
3968             Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
3969                                                  Quad1.Normale(ptvalid));
3970             if(qwe> 0.00000001) {
3971               trans1 = IntSurf_Out;
3972               trans2 = IntSurf_In;
3973             }
3974             else if(qwe<-0.00000001) {
3975               trans1 = IntSurf_In;
3976               trans2 = IntSurf_Out;
3977             }
3978             else {
3979               trans1=trans2=IntSurf_Undecided;
3980             }
3981             alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
3982           }
3983           else {
3984             alig = new IntPatch_ALine(curvsol,Standard_False);
3985           }
3986           Standard_Boolean TempFalse1a = Standard_False;
3987           Standard_Boolean TempFalse2a = Standard_False;
3988
3989           //-- ptf et ptl : points debut et fin de alig
3990
3991           ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
3992                         TempFalse2a,ptl,last,Multpoint,Tol);
3993           slin.Append(alig);
3994         } //-- boucle sur les lignes
3995       } //-- solution avec au moins une lihne
3996     }
3997     break;
3998
3999   default:
4000     {
4001       return Standard_False;
4002     }
4003   }
4004   return Standard_True;
4005 }
4006 //=======================================================================
4007 //function : IntCyCo
4008 //purpose  :
4009 //=======================================================================
4010 Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
4011                          const IntSurf_Quadric& Quad2,
4012                          const Standard_Real Tol,
4013                          const Standard_Boolean Reversed,
4014                          Standard_Boolean& Empty,
4015                          Standard_Boolean& Multpoint,
4016                          IntPatch_SequenceOfLine& slin,
4017                          IntPatch_SequenceOfPoint& spnt)
4018
4019 {
4020   IntPatch_Point ptsol;
4021
4022   Standard_Integer i;
4023
4024   IntSurf_TypeTrans trans1,trans2;
4025   IntAna_ResultType typint;
4026   gp_Circ cirsol;
4027
4028   gp_Cylinder Cy;
4029   gp_Cone     Co;
4030
4031   if (!Reversed) {
4032     Cy = Quad1.Cylinder();
4033     Co = Quad2.Cone();
4034   }
4035   else {
4036     Cy = Quad2.Cylinder();
4037     Co = Quad1.Cone();
4038   }
4039   IntAna_QuadQuadGeo inter(Cy,Co,Tol);
4040
4041   if (!inter.IsDone()) {return Standard_False;}
4042
4043   typint = inter.TypeInter();
4044   Standard_Integer NbSol = inter.NbSolutions();
4045   Empty = Standard_False;
4046
4047   switch (typint) {
4048
4049   case IntAna_Empty : {
4050     Empty = Standard_True;
4051   }
4052     break;
4053
4054   case IntAna_Point :{
4055     gp_Pnt psol(inter.Point(1));
4056     Standard_Real U1,V1,U2,V2;
4057     Quad1.Parameters(psol,U1,V1);
4058     Quad1.Parameters(psol,U2,V2);
4059     ptsol.SetValue(psol,Tol,Standard_True);
4060     ptsol.SetParameters(U1,V1,U2,V2);
4061     spnt.Append(ptsol);
4062   }
4063     break;
4064
4065   case IntAna_Circle:  {
4066     gp_Vec Tgt;
4067     gp_Pnt ptref;
4068     Standard_Integer j;
4069     Standard_Real qwe;
4070     //
4071     for(j=1; j<=2; ++j) {
4072       cirsol = inter.Circle(j);
4073       ElCLib::D1(0.,cirsol,ptref,Tgt);
4074       qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
4075       if(qwe> 0.00000001) {
4076         trans1 = IntSurf_Out;
4077         trans2 = IntSurf_In;
4078       }
4079       else if(qwe<-0.00000001) {
4080         trans1 = IntSurf_In;
4081         trans2 = IntSurf_Out;
4082       }
4083       else {
4084         trans1=trans2=IntSurf_Undecided;
4085       }
4086       Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
4087       slin.Append(glig);
4088     }
4089   }
4090     break;
4091
4092   case IntAna_NoGeometricSolution:    {
4093     gp_Pnt psol;
4094     Standard_Real U1,V1,U2,V2;
4095     IntAna_IntQuadQuad anaint(Cy,Co,Tol);
4096     if (!anaint.IsDone()) {
4097       return Standard_False;
4098     }
4099
4100     if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
4101       Empty = Standard_True;
4102     }
4103     else {
4104       NbSol = anaint.NbPnt();
4105       for (i = 1; i <= NbSol; i++) {
4106         psol = anaint.Point(i);
4107         Quad1.Parameters(psol,U1,V1);
4108         Quad2.Parameters(psol,U2,V2);
4109         ptsol.SetValue(psol,Tol,Standard_True);
4110         ptsol.SetParameters(U1,V1,U2,V2);
4111         spnt.Append(ptsol);
4112       }
4113
4114       gp_Pnt ptvalid, ptf, ptl;
4115       gp_Vec tgvalid;
4116       //
4117       Standard_Real first,last,para;
4118       Standard_Boolean tgfound,firstp,lastp,kept;
4119       Standard_Integer kount;
4120       //
4121       //
4122       //IntAna_Curve curvsol;
4123       IntAna_Curve aC;
4124       IntAna_ListOfCurve aLC;
4125       IntAna_ListIteratorOfListOfCurve aIt;
4126
4127       //
4128       NbSol = anaint.NbCurve();
4129       for (i = 1; i <= NbSol; ++i) {
4130         kept = Standard_False;
4131         //curvsol = anaint.Curve(i);
4132         aC=anaint.Curve(i);
4133         aLC.Clear();
4134         ExploreCurve(Co, aC, 10.*Tol, aLC);
4135         //
4136         aIt.Initialize(aLC);
4137         for (; aIt.More(); aIt.Next()) {
4138           IntAna_Curve& curvsol=aIt.Value();
4139           //
4140           curvsol.Domain(first, last);
4141           firstp = !curvsol.IsFirstOpen();
4142           lastp  = !curvsol.IsLastOpen();
4143           if (firstp) {
4144             ptf = curvsol.Value(first);
4145           }
4146           if (lastp) {
4147             ptl = curvsol.Value(last);
4148           }
4149           para = last;
4150           kount = 1;
4151           tgfound = Standard_False;
4152
4153           while (!tgfound) {
4154             para = (1.123*first + para)/2.123;
4155             tgfound = curvsol.D1u(para,ptvalid,tgvalid);
4156             if (!tgfound) {
4157               kount ++;
4158               tgfound = kount > 5;
4159             }
4160           }
4161           Handle(IntPatch_ALine) alig;
4162           if (kount <= 5) {
4163             Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
4164                                                  Quad1.Normale(ptvalid));
4165             if(qwe> 0.00000001) {
4166               trans1 = IntSurf_Out;
4167               trans2 = IntSurf_In;
4168             }
4169             else if(qwe<-0.00000001) {
4170               trans1 = IntSurf_In;
4171               trans2 = IntSurf_Out;
4172             }
4173             else {
4174               trans1=trans2=IntSurf_Undecided;
4175             }
4176             alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
4177             kept = Standard_True;
4178           }
4179           else {
4180             ptvalid = curvsol.Value(para);
4181             alig = new IntPatch_ALine(curvsol,Standard_False);
4182             kept = Standard_True;
4183             //-- std::cout << "Transition indeterminee" << std::endl;
4184           }
4185           if (kept) {
4186             Standard_Boolean Nfirstp = !firstp;
4187             Standard_Boolean Nlastp  = !lastp;
4188             ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
4189                           Nlastp,ptl,last,Multpoint,Tol);
4190             slin.Append(alig);
4191           }
4192         } // for (; aIt.More(); aIt.Next())
4193       } // for (i = 1; i <= NbSol; ++i)
4194     }
4195   }
4196     break;
4197
4198   default:
4199     return Standard_False;
4200
4201   } // switch (typint)
4202
4203   return Standard_True;
4204 }
4205 //=======================================================================
4206 //function : ExploreCurve
4207 //purpose  : Splits aC on several curves in the cone apex points.
4208 //=======================================================================
4209 Standard_Boolean ExploreCurve(const gp_Cone& theCo,
4210                               IntAna_Curve& theCrv,
4211                               const Standard_Real theTol,
4212                               IntAna_ListOfCurve& theLC)
4213 {
4214   const Standard_Real aSqTol = theTol*theTol;
4215   const gp_Pnt aPapx(theCo.Apex());
4216
4217   Standard_Real aT1, aT2;
4218   theCrv.Domain(aT1, aT2);
4219
4220   theLC.Clear();
4221   //
4222   TColStd_ListOfReal aLParams;
4223   theCrv.FindParameter(aPapx, aLParams);
4224   if (aLParams.IsEmpty())
4225   {
4226     theLC.Append(theCrv);
4227     return Standard_False;
4228   }
4229
4230   for (TColStd_ListIteratorOfListOfReal anItr(aLParams); anItr.More(); anItr.Next())
4231   {
4232     Standard_Real aPrm = anItr.Value();
4233
4234     if ((aPrm - aT1) < Precision::PConfusion())
4235       continue;
4236
4237     Standard_Boolean isLast = Standard_False;
4238     if ((aT2 - aPrm) < Precision::PConfusion())
4239     {
4240       aPrm = aT2;
4241       isLast = Standard_True;
4242     }
4243
4244     const gp_Pnt aP = theCrv.Value(aPrm);
4245     const Standard_Real aSqD = aP.SquareDistance(aPapx);
4246     if (aSqD < aSqTol)
4247     {
4248       IntAna_Curve aC1 = theCrv;
4249       aC1.SetDomain(aT1, aPrm);
4250       aT1 = aPrm;
4251       theLC.Append(aC1);
4252     }
4253
4254     if (isLast)
4255       break;
4256   }
4257
4258   if (theLC.IsEmpty())
4259   {
4260     theLC.Append(theCrv);
4261     return Standard_False;
4262   }
4263
4264   if ((aT2 - aT1) > Precision::PConfusion())
4265   {
4266     IntAna_Curve aC1 = theCrv;
4267     aC1.SetDomain(aT1, aT2);
4268     theLC.Append(aC1);
4269   }
4270
4271   return Standard_True;
4272 }