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