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