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