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