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 = 0.0;
 754
 755   const Standard_Real aDelta1 = mVecC1.X()*mVecC2.Y()-mVecC1.Y()*mVecC2.X(); //1-2
 756   const Standard_Real aDelta2 = mVecC1.Y()*mVecC2.Z()-mVecC1.Z()*mVecC2.Y(); //2-3
 757   const Standard_Real aDelta3 = mVecC1.X()*mVecC2.Z()-mVecC1.Z()*mVecC2.X(); //1-3
 758   const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
 759   const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
 760   const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
 761
 762   if(anAbsD1 >= anAbsD2)
 763   {
 764     if(anAbsD3 > anAbsD1)
 765     {
 766       aFoundCouple = COE13;
 767       aDetV1V2 = aDelta3;
 768     }
 769     else
 770     {
 771       aFoundCouple = COE12;
 772       aDetV1V2 = aDelta1;
 773     }
 774   }
 775   else
 776   {
 777     if(anAbsD3 > anAbsD2)
 778     {
 779       aFoundCouple = COE13;
 780       aDetV1V2 = aDelta3;
 781     }
 782     else
 783     {
 784       aFoundCouple = COE23;
 785       aDetV1V2 = aDelta2;
 786     }
 787   }
 788
 789   if(Abs(aDetV1V2) < aNulValue)
 790   {
 791     Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
 792   }
 793
 794   switch(aFoundCouple)
 795   {
 796   case COE12:
 797     break;
 798   case COE23:
 799     {
 800       gp_Vec aVTemp = mVecA1;
 801       mVecA1.SetX(aVTemp.Y());
 802       mVecA1.SetY(aVTemp.Z());
 803       mVecA1.SetZ(aVTemp.X());
 804
 805       aVTemp = mVecA2;
 806       mVecA2.SetX(aVTemp.Y());
 807       mVecA2.SetY(aVTemp.Z());
 808       mVecA2.SetZ(aVTemp.X());
 809
 810       aVTemp = mVecB1;
 811       mVecB1.SetX(aVTemp.Y());
 812       mVecB1.SetY(aVTemp.Z());
 813       mVecB1.SetZ(aVTemp.X());
 814
 815       aVTemp = mVecB2;
 816       mVecB2.SetX(aVTemp.Y());
 817       mVecB2.SetY(aVTemp.Z());
 818       mVecB2.SetZ(aVTemp.X());
 819
 820       aVTemp = mVecC1;
 821       mVecC1.SetX(aVTemp.Y());
 822       mVecC1.SetY(aVTemp.Z());
 823       mVecC1.SetZ(aVTemp.X());
 824
 825       aVTemp = mVecC2;
 826       mVecC2.SetX(aVTemp.Y());
 827       mVecC2.SetY(aVTemp.Z());
 828       mVecC2.SetZ(aVTemp.X());
 829
 830       aVTemp = mVecD;
 831       mVecD.SetX(aVTemp.Y());
 832       mVecD.SetY(aVTemp.Z());
 833       mVecD.SetZ(aVTemp.X());
 834
 835       break;
 836     }
 837   case COE13:
 838     {
 839       gp_Vec aVTemp = mVecA1;
 840       mVecA1.SetY(aVTemp.Z());
 841       mVecA1.SetZ(aVTemp.Y());
 842
 843       aVTemp = mVecA2;
 844       mVecA2.SetY(aVTemp.Z());
 845       mVecA2.SetZ(aVTemp.Y());
 846
 847       aVTemp = mVecB1;
 848       mVecB1.SetY(aVTemp.Z());
 849       mVecB1.SetZ(aVTemp.Y());
 850
 851       aVTemp = mVecB2;
 852       mVecB2.SetY(aVTemp.Z());
 853       mVecB2.SetZ(aVTemp.Y());
 854
 855       aVTemp = mVecC1;
 856       mVecC1.SetY(aVTemp.Z());
 857       mVecC1.SetZ(aVTemp.Y());
 858
 859       aVTemp = mVecC2;
 860       mVecC2.SetY(aVTemp.Z());
 861       mVecC2.SetZ(aVTemp.Y());
 862
 863       aVTemp = mVecD;
 864       mVecD.SetY(aVTemp.Z());
 865       mVecD.SetZ(aVTemp.Y());
 866
 867       break;
 868     }
 869   default:
 870     break;
 871   }
 872
 873   //------- For V1 (begin)
 874   //sinU2
 875   mK21 = (mVecC2.Y()*mVecB2.X()-mVecC2.X()*mVecB2.Y())/aDetV1V2;
 876   //sinU1
 877   mK11 = (mVecC2.Y()*mVecB1.X()-mVecC2.X()*mVecB1.Y())/aDetV1V2;
 878   //cosU2
 879   mL21 = (mVecC2.Y()*mVecA2.X()-mVecC2.X()*mVecA2.Y())/aDetV1V2;
 880   //cosU1
 881   mL11 = (mVecC2.Y()*mVecA1.X()-mVecC2.X()*mVecA1.Y())/aDetV1V2;
 882   //Free member
 883   mM1 = (mVecC2.Y()*mVecD.X()-mVecC2.X()*mVecD.Y())/aDetV1V2;
 884   //------- For V1 (end)
 885
 886   //------- For V2 (begin)
 887   //sinU2
 888   mK22 = (mVecC1.X()*mVecB2.Y()-mVecC1.Y()*mVecB2.X())/aDetV1V2;
 889   //sinU1
 890   mK12 = (mVecC1.X()*mVecB1.Y()-mVecC1.Y()*mVecB1.X())/aDetV1V2;
 891   //cosU2
 892   mL22 = (mVecC1.X()*mVecA2.Y()-mVecC1.Y()*mVecA2.X())/aDetV1V2;
 893   //cosU1
 894   mL12 = (mVecC1.X()*mVecA1.Y()-mVecC1.Y()*mVecA1.X())/aDetV1V2;
 895   //Free member
 896   mM2 = (mVecC1.X()*mVecD.Y()-mVecC1.Y()*mVecD.X())/aDetV1V2;
 897   //------- For V1 (end)
 898
 899   ShortCosForm(mL11, mK11, mK1, mFIV1);
 900   ShortCosForm(mL21, mK21, mL1, mPSIV1);
 901   ShortCosForm(mL12, mK12, mK2, mFIV2);
 902   ShortCosForm(mL22, mK22, mL2, mPSIV2);
 903
 904   const Standard_Real aA1=mVecC1.Z()*mK21+mVecC2.Z()*mK22-mVecB2.Z(), //sinU2
 905                       aA2=mVecC1.Z()*mL21+mVecC2.Z()*mL22-mVecA2.Z(), //cosU2
 906                       aB1=mVecB1.Z()-mVecC1.Z()*mK11-mVecC2.Z()*mK12, //sinU1
 907                       aB2=mVecA1.Z()-mVecC1.Z()*mL11-mVecC2.Z()*mL12; //cosU1
 908
 909   mC =mVecD.Z() -mVecC1.Z()*mM1 -mVecC2.Z()*mM2;  //Free
 910
 911   Standard_Real aA = 0.0;
 912
 913   ShortCosForm(aB2,aB1,mB,mFI1);
 914   ShortCosForm(aA2,aA1,aA,mFI2);
 915
 916   mB /= aA;
 917   mC /= aA;
 918 }
 919
 920 //=======================================================================
 921 //function : SearchOnVBounds
 922 //purpose  :
 923 //=======================================================================
 924 static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
 925                                         const stCoeffsValue& theCoeffs,
 926                                         const Standard_Real theVzad,
 927                                         const Standard_Real theInitU2,
 928                                         const Standard_Real theInitMainVar,
 929                                         Standard_Real& theMainVariableValue)
 930 {
 931   const Standard_Real aMaxError = 4.0*M_PI; // two periods
 932
 933   theMainVariableValue = RealLast();
 934   const Standard_Real aTol2 = Precision::PConfusion()*Precision::PConfusion();
 935   Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVzad;
 936
 937   Standard_Real anError = RealLast();
 938   Standard_Integer aNbIter = 0;
 939   do
 940   {
 941     if(++aNbIter > 1000)
 942       return Standard_False;
 943
 944     gp_Vec aVecMainVar = theCoeffs.mVecA1 * sin(aMainVarPrev) - theCoeffs.mVecB1 * cos(aMainVarPrev);
 945     gp_Vec aVecFreeMem =  (theCoeffs.mVecA2 * aU2Prev + theCoeffs.mVecB2) * sin(aU2Prev) -
 946                           (theCoeffs.mVecB2 * aU2Prev - theCoeffs.mVecA2) * cos(aU2Prev) +
 947                           (theCoeffs.mVecA1 * aMainVarPrev + theCoeffs.mVecB1) * sin(aMainVarPrev) -
 948                           (theCoeffs.mVecB1 * aMainVarPrev - theCoeffs.mVecA1) * cos(aMainVarPrev) + theCoeffs.mVecD;
 949
 950     gp_Vec aVecVar1 = theCoeffs.mVecA2 * sin(aU2Prev) - theCoeffs.mVecB2 * cos(aU2Prev);
 951     gp_Vec aVecVar2;
 952
 953     switch(theSBType)
 954     {
 955     case SearchV1:
 956       aVecVar2 = theCoeffs.mVecC2;
 957       aVecFreeMem -= theCoeffs.mVecC1 * theVzad;
 958       break;
 959
 960     case SearchV2:
 961       aVecVar2 = theCoeffs.mVecC1;
 962       aVecFreeMem -= theCoeffs.mVecC2 * theVzad;
 963       break;
 964
 965     default:
 966       return Standard_False;
 967     }
 968
 969     Standard_Real aDetMainSyst =  aVecMainVar.X()*(aVecVar1.Y()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.Y())-
 970                                   aVecMainVar.Y()*(aVecVar1.X()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.X())+
 971                                   aVecMainVar.Z()*(aVecVar1.X()*aVecVar2.Y()-aVecVar1.Y()*aVecVar2.X());
 972
 973     if(IsEqual(aDetMainSyst, 0.0))
 974     {
 975       return Standard_False;
 976     }
 977
 978
 979     Standard_Real aDetMainVar = aVecFreeMem.X()*(aVecVar1.Y()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.Y())-
 980                                 aVecFreeMem.Y()*(aVecVar1.X()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.X())+
 981                                 aVecFreeMem.Z()*(aVecVar1.X()*aVecVar2.Y()-aVecVar1.Y()*aVecVar2.X());
 982
 983     Standard_Real aDetVar1    = aVecMainVar.X()*(aVecFreeMem.Y()*aVecVar2.Z()-aVecFreeMem.Z()*aVecVar2.Y())-
 984                                 aVecMainVar.Y()*(aVecFreeMem.X()*aVecVar2.Z()-aVecFreeMem.Z()*aVecVar2.X())+
 985                                 aVecMainVar.Z()*(aVecFreeMem.X()*aVecVar2.Y()-aVecFreeMem.Y()*aVecVar2.X());
 986
 987     Standard_Real aDetVar2    = aVecMainVar.X()*(aVecVar1.Y()*aVecFreeMem.Z()-aVecVar1.Z()*aVecFreeMem.Y())-
 988                                 aVecMainVar.Y()*(aVecVar1.X()*aVecFreeMem.Z()-aVecVar1.Z()*aVecFreeMem.X())+
 989                                 aVecMainVar.Z()*(aVecVar1.X()*aVecFreeMem.Y()-aVecVar1.Y()*aVecFreeMem.X());
 990
 991     Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
 992
 993     if(Abs(aDelta) > aMaxError)
 994       return Standard_False;
 995
 996     anError = aDelta*aDelta;
 997     aMainVarPrev += aDelta;
 998
 999     ///
1000     aDelta = aDetVar1/aDetMainSyst-aU2Prev;
1001
1002     if(Abs(aDelta) > aMaxError)
1003       return Standard_False;
1004
1005     anError += aDelta*aDelta;
1006     aU2Prev += aDelta;
1007
1008     ///
1009     aDelta = aDetVar2/aDetMainSyst-anOtherVar;
1010     anError += aDelta*aDelta;
1011     anOtherVar += aDelta;
1012   }
1013   while(anError > aTol2);
1014
1015   theMainVariableValue = aMainVarPrev;
1016
1017   return Standard_True;
1018 }
1019
1020 //=======================================================================
1021 //function : InscribePoint
1022 //purpose  : If theFlForce==TRUE theUGiven will be adjasted forceful.
1023 //=======================================================================
1024 static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
1025                                       const Standard_Real theUlTarget,
1026                                       Standard_Real& theUGiven,
1027                                       const Standard_Real theTol2D,
1028                                       const Standard_Real thePeriod,
1029                                       const Standard_Boolean theFlForce)
1030 {
1031   if((theUfTarget - theUGiven <= theTol2D) &&
1032               (theUGiven - theUlTarget <= theTol2D))
1033   {//it has already inscribed
1034
1035     /*
1036              Utf    U      Utl
1037               +     *       +
1038     */
1039
1040     if(theFlForce)
1041     {
1042       Standard_Real anUtemp = theUGiven + thePeriod;
1043       if((theUfTarget - anUtemp <= theTol2D) &&
1044                 (anUtemp - theUlTarget <= theTol2D))
1045       {
1046         theUGiven = anUtemp;
1047         return Standard_True;
1048       }
1049
1050       anUtemp = theUGiven - thePeriod;
1051       if((theUfTarget - anUtemp <= theTol2D) &&
1052                 (anUtemp - theUlTarget <= theTol2D))
1053       {
1054         theUGiven = anUtemp;
1055       }
1056     }
1057
1058     return Standard_True;
1059   }
1060
1061   if(IsEqual(thePeriod, 0.0))
1062   {//it cannot be inscribed
1063     return Standard_False;
1064   }
1065
1066   Standard_Real aDU = theUGiven - theUfTarget;
1067
1068   if(aDU > 0.0)
1069     aDU = -thePeriod;
1070   else
1071     aDU = +thePeriod;
1072
1073   while(((theUGiven - theUfTarget)*aDU < 0.0) && !((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)))
1074   {
1075     theUGiven += aDU;
1076   }
1077
1078   return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D));
1079 }
1080
1081 //=======================================================================
1082 //function : InscribeInterval
1083 //purpose  : Adjusts theUfGiven and after that fits theUlGiven to result
1084 //=======================================================================
1085 static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
1086                                       const Standard_Real theUlTarget,
1087                                       Standard_Real& theUfGiven,
1088                                       Standard_Real& theUlGiven,
1089                                       const Standard_Real theTol2D,
1090                                       const Standard_Real thePeriod)
1091 {
1092   Standard_Real anUpar = theUfGiven;
1093   const Standard_Real aDelta = theUlGiven - theUfGiven;
1094   if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1095                         theTol2D, thePeriod, Standard_False))
1096   {
1097     theUfGiven = anUpar;
1098     theUlGiven = theUfGiven + aDelta;
1099   }
1100   else
1101   {
1102     anUpar = theUlGiven;
1103     if(InscribePoint(theUfTarget, theUlTarget, anUpar,
1104                         theTol2D, thePeriod, Standard_False))
1105     {
1106       theUlGiven = anUpar;
1107       theUfGiven = theUlGiven - aDelta;
1108     }
1109     else
1110     {
1111       return Standard_False;
1112     }
1113   }
1114
1115   return Standard_True;
1116 }
1117
1118 //=======================================================================
1119 //function : InscribeAndSortArray
1120 //purpose  : Sort from Min to Max value
1121 //=======================================================================
1122 static void InscribeAndSortArray( Standard_Real theArr[],
1123                                   const Standard_Integer theNOfMembers,
1124                                   const Standard_Real theUf,
1125                                   const Standard_Real theUl,
1126                                   const Standard_Real theTol2D,
1127                                   const Standard_Real thePeriod)
1128 {
1129   for(Standard_Integer i = 0; i < theNOfMembers; i++)
1130   {
1131     if(Precision::IsInfinite(theArr[i]))
1132     {
1133       if(theArr[i] < 0.0)
1134         theArr[i] = -theArr[i];
1135
1136       continue;
1137     }
1138
1139     InscribePoint(theUf, theUl, theArr[i], theTol2D, thePeriod, Standard_False);
1140
1141     for(Standard_Integer j = i - 1; j >= 0; j--)
1142     {
1143
1144       if(theArr[j+1] < theArr[j])
1145       {
1146         Standard_Real anUtemp = theArr[j+1];
1147         theArr[j+1] = theArr[j];
1148         theArr[j] = anUtemp;
1149       }
1150     }
1151   }
1152 }
1153
1154
1155 //=======================================================================
1156 //function : AddPointIntoWL
1157 //purpose  : Surf1 is a surface, whose U-par is variable.
1158 //=======================================================================
1159 static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
1160                                         const IntSurf_Quadric& theQuad2,
1161                                         const Standard_Boolean isTheReverse,
1162                                         const gp_Pnt2d& thePntOnSurf1,
1163                                         const gp_Pnt2d& thePntOnSurf2,
1164                                         const Standard_Real theUfSurf1,
1165                                         const Standard_Real theUlSurf1,
1166                                         const Standard_Real thePeriodOfSurf1,
1167                                         const Handle(IntSurf_LineOn2S)& theLine,
1168                                         const Standard_Real theTol3D,
1169                                         const Standard_Real theTol2D,
1170                                         const Standard_Boolean theFlForce)
1171 {
1172   gp_Pnt  aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
1173           aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
1174
1175   Standard_Real anUpar = thePntOnSurf1.X();
1176   if(!InscribePoint(theUfSurf1, theUlSurf1, anUpar, theTol2D,
1177                                     thePeriodOfSurf1, theFlForce))
1178     return Standard_False;
1179
1180   IntSurf_PntOn2S aPnt;
1181
1182   if(isTheReverse)
1183   {
1184     aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
1185                         thePntOnSurf2.X(), thePntOnSurf2.Y(),
1186                         anUpar, thePntOnSurf1.Y());
1187   }
1188   else
1189   {
1190     aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
1191                         anUpar, thePntOnSurf1.Y(),
1192                         thePntOnSurf2.X(), thePntOnSurf2.Y());
1193   }
1194
1195   const Standard_Integer aNbPnts = theLine->NbPoints();
1196   if(aNbPnts > 0)
1197   {
1198     Standard_Real aUl = 0.0, aVl = 0.0;
1199     const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts);
1200     if(isTheReverse)
1201       aPlast.ParametersOnS2(aUl, aVl);
1202     else
1203       aPlast.ParametersOnS1(aUl, aVl);
1204
1205     if(anUpar <= aUl)
1206     {//Parameter value will be always increased.
1207       return Standard_False;
1208     }
1209
1210     //theTol2D is minimal step along parameter changed.
1211     //Therefore, if we apply this minimal step two
1212     //neighbour points will be always "same". Consequently,
1213     //we should reduce tolerance for IsSame checking.
1214     const Standard_Real aDTol = 1.0-Epsilon(1.0);
1215     if(aPnt.IsSame(aPlast, theTol3D*aDTol, theTol2D*aDTol))
1216     {
1217       theLine->RemovePoint(aNbPnts);
1218     }
1219   }
1220
1221   theLine->Add(aPnt);
1222   return Standard_True;
1223 }
1224
1225 //=======================================================================
1226 //function : AddBoundaryPoint
1227 //purpose  :
1228 //=======================================================================
1229 static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
1230                                           const IntSurf_Quadric& theQuad2,
1231                                           const Handle(IntPatch_WLine)& theWL,
1232                                           const stCoeffsValue& theCoeffs,
1233                                           const Bnd_Box2d& theUVSurf1,
1234                                           const Bnd_Box2d& theUVSurf2,
1235                                           const Standard_Real theTol3D,
1236                                           const Standard_Real theTol2D,
1237                                           const Standard_Real thePeriod,
1238                                           const Standard_Real theNulValue,
1239                                           const Standard_Real theU1,
1240                                           const Standard_Real theU2,
1241                                           const Standard_Real theV1,
1242                                           const Standard_Real theV1Prev,
1243                                           const Standard_Real theV2,
1244                                           const Standard_Real theV2Prev,
1245                                           const Standard_Boolean isTheReverse,
1246                                           const Standard_Real theArccosFactor,
1247                                           const Standard_Boolean theFlForce,
1248                                           Standard_Boolean& isTheFound1,
1249                                           Standard_Boolean& isTheFound2)
1250 {
1251   Standard_Real aUSurf1f = 0.0, //const
1252                 aUSurf1l = 0.0,
1253                 aVSurf1f = 0.0,
1254                 aVSurf1l = 0.0;
1255   Standard_Real aUSurf2f = 0.0, //const
1256                 aUSurf2l = 0.0,
1257                 aVSurf2f = 0.0,
1258                 aVSurf2l = 0.0;
1259
1260   theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
1261   theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
1262
1263   SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE;
1264   Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f;
1265
1266   Standard_Real anUpar1 = theU1, anUpar2 = theU1;
1267   Standard_Real aVf = theV1, aVl = theV1Prev;
1268   MinMax(aVf, aVl);
1269   if((aVf <= aVSurf1f) && (aVSurf1f <= aVl))
1270   {
1271     aTS1 = SearchV1;
1272     aV1zad = aVSurf1f;
1273     isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theU2, theU1, anUpar1);
1274   }
1275   else if((aVf <= aVSurf1l) && (aVSurf1l <= aVl))
1276   {
1277     aTS1 = SearchV1;
1278     aV1zad = aVSurf1l;
1279     isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theU2, theU1, anUpar1);
1280   }
1281
1282   aVf = theV2;
1283   aVl = theV2Prev;
1284   MinMax(aVf, aVl);
1285
1286   if((aVf <= aVSurf2f) && (aVSurf2f <= aVl))
1287   {
1288     aTS2 = SearchV2;
1289     aV2zad = aVSurf2f;
1290     isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theU2, theU1, anUpar2);
1291   }
1292   else if((aVf <= aVSurf2l) && (aVSurf2l <= aVl))
1293   {
1294     aTS2 = SearchV2;
1295     aV2zad = aVSurf2l;
1296     isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theU2, theU1, anUpar2);
1297   }
1298
1299   if(anUpar1 < anUpar2)
1300   {
1301     if(isTheFound1)
1302     {
1303       Standard_Real anArg = theCoeffs.mB * cos(anUpar1 - theCoeffs.mFI1) + theCoeffs.mC;
1304
1305       if(theNulValue > 1.0 - anArg)
1306         anArg = 1.0;
1307       if(anArg + 1.0 < theNulValue)
1308         anArg = -1.0;
1309
1310       Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
1311
1312       if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False))
1313       {
1314         const Standard_Real aV1 =
1315               (aTS1 == SearchV1) ? aV1zad :
1316                                   theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar1) +
1317                                   theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar1) + theCoeffs.mM1;
1318         const Standard_Real aV2 =
1319               (aTS1 == SearchV2) ? aV2zad :
1320                                   theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar1) +
1321                                   theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar1) + theCoeffs.mM2;
1322
1323         AddPointIntoWL(theQuad1, theQuad2, isTheReverse,
1324                           gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2),
1325                           aUSurf1f, aUSurf1l, thePeriod,
1326                           theWL->Curve(), theTol3D, theTol2D, theFlForce);
1327       }
1328       else
1329       {
1330         isTheFound1 = Standard_False;
1331       }
1332     }
1333
1334     if(isTheFound2)
1335     {
1336       Standard_Real anArg = theCoeffs.mB * cos(anUpar2 - theCoeffs.mFI1) + theCoeffs.mC;
1337
1338       if(theNulValue > 1.0 - anArg)
1339         anArg = 1.0;
1340       if(anArg + 1.0 < theNulValue)
1341         anArg = -1.0;
1342
1343       Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
1344       if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False))
1345       {
1346         const Standard_Real aV1 =
1347               (aTS2 == SearchV1) ? aV1zad :
1348                                   theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar2) +
1349                                   theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar2) + theCoeffs.mM1;
1350         const Standard_Real aV2 =
1351               (aTS2 == SearchV2) ? aV2zad :
1352                                   theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar2) +
1353                                   theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar2) + theCoeffs.mM2;
1354
1355         AddPointIntoWL(theQuad1, theQuad2, isTheReverse,
1356                         gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2),
1357                         aUSurf1f, aUSurf1l, thePeriod,
1358                         theWL->Curve(),theTol3D, theTol2D, theFlForce);
1359       }
1360       else
1361       {
1362         isTheFound2 = Standard_False;
1363       }
1364     }
1365   }
1366   else
1367   {
1368     if(isTheFound2)
1369     {
1370       Standard_Real anArg = theCoeffs.mB * cos(anUpar2 - theCoeffs.mFI1) + theCoeffs.mC;
1371
1372       if(theNulValue > 1.0 - anArg)
1373         anArg = 1.0;
1374       if(anArg + 1.0 < theNulValue)
1375         anArg = -1.0;
1376
1377       Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
1378
1379       if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False))
1380       {
1381         const Standard_Real aV1 = (aTS2 == SearchV1) ? aV1zad :
1382                                   theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar2) +
1383                                   theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar2) + theCoeffs.mM1;
1384         const Standard_Real aV2 = (aTS2 == SearchV2) ? aV2zad :
1385                                   theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar2) +
1386                                   theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar2) + theCoeffs.mM2;
1387
1388         AddPointIntoWL(theQuad1, theQuad2, isTheReverse,
1389                         gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2),
1390                         aUSurf1f, aUSurf1l, thePeriod,
1391                         theWL->Curve(), theTol3D, theTol2D, theFlForce);
1392       }
1393       else
1394       {
1395         isTheFound2 = Standard_False;
1396       }
1397     }
1398
1399     if(isTheFound1)
1400     {
1401       Standard_Real anArg = theCoeffs.mB*cos(anUpar1-theCoeffs.mFI1)+theCoeffs.mC;
1402
1403       if(theNulValue > 1.0 - anArg)
1404         anArg = 1.0;
1405       if(anArg + 1.0 < theNulValue)
1406         anArg = -1.0;
1407
1408       Standard_Real aU2 = theCoeffs.mFI2+theArccosFactor*acos(anArg);
1409       if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False))
1410       {
1411         const Standard_Real aV1 = (aTS1 == SearchV1) ? aV1zad :
1412                                   theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar1) +
1413                                   theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar1) + theCoeffs.mM1;
1414         const Standard_Real aV2 = (aTS1 == SearchV2) ? aV2zad :
1415                                   theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar1) +
1416                                   theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar1) + theCoeffs.mM2;
1417
1418         AddPointIntoWL(theQuad1, theQuad2, isTheReverse,
1419                         gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2),
1420                         aUSurf1f, aUSurf1l, thePeriod,
1421                         theWL->Curve(), theTol3D, theTol2D, theFlForce);
1422       }
1423       else
1424       {
1425         isTheFound1 = Standard_False;
1426       }
1427     }
1428   }
1429
1430   return Standard_True;
1431 }
1432
1433 //=======================================================================
1434 //function : SeekAdditionalPoints
1435 //purpose  :
1436 //=======================================================================
1437 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
1438                                   const IntSurf_Quadric& theQuad2,
1439                                   const Handle(IntSurf_LineOn2S)& theLile,
1440                                   const stCoeffsValue& theCoeffs,
1441                                   const Standard_Integer theMinNbPoints,
1442                                   const Standard_Real theU2f,
1443                                   const Standard_Real theU2l,
1444                                   const Standard_Real theTol2D,
1445                                   const Standard_Real thePeriodOfSurf2,
1446                                   const Standard_Real theArccosFactor,
1447                                   const Standard_Boolean isTheReverse)
1448 {
1449   Standard_Integer aNbPoints = theLile->NbPoints();
1450   if(aNbPoints >= theMinNbPoints)
1451   {
1452     return;
1453   }
1454
1455   Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
1456
1457   Standard_Integer aNbPointsPrev = 0;
1458   while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
1459   {
1460     aNbPointsPrev = aNbPoints;
1461     for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
1462     {
1463       Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
1464       Standard_Real U1l = 0.0, V1l = 0.0; //last  point in 1st suraface
1465
1466       Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface
1467       Standard_Real U2l = 0.0, V2l = 0.0; //last  point in 2nd suraface
1468
1469       lp = fp+1;
1470
1471       if(isTheReverse)
1472       {
1473         theLile->Value(fp).ParametersOnS2(U1f, V1f);
1474         theLile->Value(lp).ParametersOnS2(U1l, V1l);
1475
1476         theLile->Value(fp).ParametersOnS1(U2f, V2f);
1477         theLile->Value(lp).ParametersOnS1(U2l, V2l);
1478       }
1479       else
1480       {
1481         theLile->Value(fp).ParametersOnS1(U1f, V1f);
1482         theLile->Value(lp).ParametersOnS1(U1l, V1l);
1483
1484         theLile->Value(fp).ParametersOnS2(U2f, V2f);
1485         theLile->Value(lp).ParametersOnS2(U2l, V2l);
1486       }
1487
1488       if(Abs(U1l - U1f) <= theTol2D)
1489       {
1490         //Step is minimal. It is not necessary to divide it.
1491         continue;
1492       }
1493
1494       U1prec = 0.5*(U1f+U1l);
1495
1496       Standard_Real anArg = theCoeffs.mB * cos(U1prec - theCoeffs.mFI1) + theCoeffs.mC;
1497       if(anArg > 1.0)
1498         anArg = 1.0;
1499       if(anArg < -1.0)
1500         anArg = -1.0;
1501
1502       U2prec = theCoeffs.mFI2 + theArccosFactor*acos(anArg);
1503       InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False);
1504
1505       gp_Pnt aP1, aP2;
1506
1507       V1prec = theCoeffs.mK21 * sin(U2prec) +
1508                 theCoeffs.mK11 * sin(U1prec) +
1509                 theCoeffs.mL21 * cos(U2prec) +
1510                 theCoeffs.mL11 * cos(U1prec) + theCoeffs.mM1;
1511       V2prec = theCoeffs.mK22 * sin(U2prec) +
1512                 theCoeffs.mK12 * sin(U1prec) +
1513                 theCoeffs.mL22 * cos(U2prec) +
1514                 theCoeffs.mL12 * cos(U1prec) + theCoeffs.mM2;
1515
1516       aP1 = theQuad1.Value(U1prec, V1prec);
1517       aP2 = theQuad2.Value(U2prec, V2prec);
1518
1519       gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
1520
1521 #ifdef OCCT_DEBUG
1522       //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl;
1523 #endif
1524
1525       IntSurf_PntOn2S anIP;
1526       if(isTheReverse)
1527       {
1528         anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec);
1529       }
1530       else
1531       {
1532         anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
1533       }
1534
1535       theLile->InsertBefore(lp, anIP);
1536
1537       aNbPoints = theLile->NbPoints();
1538       if(aNbPoints >= theMinNbPoints)
1539       {
1540         return;
1541       }
1542     }
1543   }
1544 }
1545
1546 //=======================================================================
1547 //function : CriticalPointsComputing
1548 //purpose  :
1549 //=======================================================================
1550 static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
1551                                     const Standard_Real theUSurf1f,
1552                                     const Standard_Real theUSurf1l,
1553                                     const Standard_Real theUSurf2f,
1554                                     const Standard_Real theUSurf2l,
1555                                     const Standard_Real thePeriod,
1556                                     const Standard_Real theTol2D,
1557                                     const Standard_Integer theNbCritPointsMax,
1558                                     Standard_Real theU1crit[])
1559 {
1560   theU1crit[0] = 0.0;
1561   theU1crit[1] = thePeriod;
1562   theU1crit[2] = theUSurf1f;
1563   theU1crit[3] = theUSurf1l;
1564
1565   const Standard_Real aCOS = cos(theCoeffs.mFI2);
1566   const Standard_Real aBSB = Abs(theCoeffs.mB);
1567   if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
1568   {
1569     Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
1570     if(anArg > 1.0)
1571       anArg = 1.0;
1572     if(anArg < -1.0)
1573       anArg = -1.0;
1574
1575     theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
1576     theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
1577   }
1578
1579   Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
1580   Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
1581   MinMax(aSf, aSl);
1582
1583   theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
1584   theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
1585   theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?  acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
1586   theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?  acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
1587
1588   //preparative treatment of array
1589   InscribeAndSortArray(theU1crit, theNbCritPointsMax, 0.0, thePeriod, theTol2D, thePeriod);
1590   for(Standard_Integer i = 1; i < theNbCritPointsMax; i++)
1591   {
1592     Standard_Real &a = theU1crit[i],
1593                   &b = theU1crit[i-1];
1594     const Standard_Real aRemain = fmod(Abs(a - b), thePeriod); // >= 0, because Abs(a - b) >= 0
1595     if((Abs(a - b) < theTol2D) || (aRemain < theTol2D) || (Abs(aRemain - thePeriod) < theTol2D))
1596     {
1597       a = (a + b)/2.0;
1598       b = Precision::Infinite();
1599     }
1600   }
1601 }
1602
1603 //=======================================================================
1604 //function : IntCyCyTrim
1605 //purpose  :
1606 //=======================================================================
1607 Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
1608                               const IntSurf_Quadric& theQuad2,
1609                               const Standard_Real theTol3D,
1610                               const Standard_Real theTol2D,
1611                               const Bnd_Box2d& theUVSurf1,
1612                               const Bnd_Box2d& theUVSurf2,
1613                               const Standard_Boolean isTheReverse,
1614                               Standard_Boolean& isTheEmpty,
1615                               IntPatch_SequenceOfLine& theSlin,
1616                               IntPatch_SequenceOfPoint& theSPnt)
1617 {
1618   Standard_Real aUSurf1f = 0.0, //const
1619                 aUSurf1l = 0.0,
1620                 aVSurf1f = 0.0,
1621                 aVSurf1l = 0.0;
1622   Standard_Real aUSurf2f = 0.0, //const
1623                 aUSurf2l = 0.0,
1624                 aVSurf2f = 0.0,
1625                 aVSurf2l = 0.0;
1626
1627   theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
1628   theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
1629
1630   const Standard_Real aNulValue = 0.01*Precision::PConfusion();
1631
1632   const gp_Cylinder&  aCyl1 = theQuad1.Cylinder(),
1633                       aCyl2 = theQuad2.Cylinder();
1634
1635   IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
1636
1637   if (!anInter.IsDone())
1638   {
1639     return Standard_False;
1640   }
1641
1642   IntAna_ResultType aTypInt = anInter.TypeInter();
1643
1644   if(aTypInt != IntAna_NoGeometricSolution)
1645   { //It is not necessary (because result is an analytic curve) or
1646     //it is impossible to make Walking-line.
1647
1648     return Standard_False;
1649   }
1650
1651   theSlin.Clear();
1652   theSPnt.Clear();
1653   const Standard_Integer aNbMaxPoints = 2000;
1654   const Standard_Integer aNbMinPoints = 200;
1655   const Standard_Integer aNbPoints = Min(Max(aNbMinPoints,
1656                       RealToInt(20.0*aCyl1.Radius())), aNbMaxPoints);
1657   const Standard_Real aPeriod = 2.0*M_PI;
1658   const Standard_Real aStepMin = theTol2D,
1659                       aStepMax = (aUSurf1l-aUSurf1f)/IntToReal(aNbPoints);
1660   const Standard_Integer aNbWLines = 2;
1661
1662   const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
1663
1664   //Boundaries
1665   const Standard_Integer aNbOfBoundaries = 2;
1666   Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
1667   Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
1668
1669   if(anEquationCoeffs.mB > 0.0)
1670   {
1671     if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) < -1.0)
1672     {//There is NOT intersection
1673       return Standard_True;
1674     }
1675     else if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
1676     {//U=[0;2*PI]+aFI1
1677       aU1f[0] = anEquationCoeffs.mFI1;
1678       aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
1679     }
1680     else if((1 + anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
1681             (anEquationCoeffs.mB <= 1 - anEquationCoeffs.mC))
1682     {
1683       Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
1684       if(anArg > 1.0)
1685         anArg = 1.0;
1686       if(anArg < -1.0)
1687         anArg = -1.0;
1688
1689       const Standard_Real aDAngle = acos(anArg);
1690       //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
1691       aU1f[0] = anEquationCoeffs.mFI1;
1692       aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
1693       aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
1694       aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
1695     }
1696     else if((1 - anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
1697             (anEquationCoeffs.mB <= 1 + anEquationCoeffs.mC))
1698     {
1699       Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
1700       if(anArg > 1.0)
1701         anArg = 1.0;
1702       if(anArg < -1.0)
1703         anArg = -1.0;
1704
1705       const Standard_Real aDAngle = acos(anArg);
1706       //U=[aDAngle;2*PI-aDAngle]+aFI1
1707
1708       aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
1709       aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
1710     }
1711     else if(anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
1712     {
1713       Standard_Real anArg1 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB,
1714                     anArg2 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
1715       if(anArg1 > 1.0)
1716         anArg1 = 1.0;
1717       if(anArg1 < -1.0)
1718         anArg1 = -1.0;
1719
1720       if(anArg2 > 1.0)
1721         anArg2 = 1.0;
1722       if(anArg2 < -1.0)
1723         anArg2 = -1.0;
1724
1725       const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
1726       //(U=[aDAngle1;aDAngle2]+aFI1) ||
1727       //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
1728
1729       aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
1730       aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
1731       aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
1732       aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
1733     }
1734     else
1735     {
1736       Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
1737     }
1738   }
1739   else if(anEquationCoeffs.mB < 0.0)
1740   {
1741     if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) > 1.0)
1742     {//There is NOT intersection
1743       return Standard_True;
1744     }
1745     else if(-anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
1746     {//U=[0;2*PI]+aFI1
1747       aU1f[0] = anEquationCoeffs.mFI1;
1748       aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
1749     }
1750     else if((-anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
1751             ( anEquationCoeffs.mB <= anEquationCoeffs.mC - 1))
1752     {
1753       Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
1754       if(anArg > 1.0)
1755         anArg = 1.0;
1756       if(anArg < -1.0)
1757         anArg = -1.0;
1758
1759       const Standard_Real aDAngle = acos(anArg);
1760       //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
1761
1762       aU1f[0] = anEquationCoeffs.mFI1;
1763       aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
1764       aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
1765       aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
1766     }
1767     else if((anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
1768             (anEquationCoeffs.mB <= -anEquationCoeffs.mB - 1))
1769     {
1770       Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
1771       if(anArg > 1.0)
1772         anArg = 1.0;
1773       if(anArg < -1.0)
1774         anArg = -1.0;
1775
1776       const Standard_Real aDAngle = acos(anArg);
1777       //U=[aDAngle;2*PI-aDAngle]+aFI1
1778
1779       aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
1780       aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
1781     }
1782     else if(-anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
1783     {
1784       Standard_Real anArg1 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB,
1785                     anArg2 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
1786       if(anArg1 > 1.0)
1787         anArg1 = 1.0;
1788       if(anArg1 < -1.0)
1789         anArg1 = -1.0;
1790
1791       if(anArg2 > 1.0)
1792         anArg2 = 1.0;
1793       if(anArg2 < -1.0)
1794         anArg2 = -1.0;
1795
1796       const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
1797       //(U=[aDAngle1;aDAngle2]+aFI1) ||
1798       //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
1799
1800       aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
1801       aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
1802       aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
1803       aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
1804     }
1805     else
1806     {
1807       Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
1808     }
1809   }
1810   else
1811   {
1812     Standard_Failure::Raise("Error. Exception. Unhandled case (B-parameter computation)!");
1813   }
1814
1815   for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
1816   {
1817     if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i]))
1818       continue;
1819
1820     InscribeInterval(aUSurf1f, aUSurf1l, aU1f[i], aU1l[i], theTol2D, aPeriod);
1821   }
1822
1823   if( !Precision::IsInfinite(aU1f[0]) && !Precision::IsInfinite(aU1f[1]) &&
1824       !Precision::IsInfinite(aU1l[0]) && !Precision::IsInfinite(aU1l[1]))
1825   {
1826     if( ((aU1f[1] <= aU1l[0]) || (aU1l[1] <= aU1l[0])) &&
1827         ((aU1f[0] <= aU1l[1]) || (aU1l[0] <= aU1l[1])))
1828     {//Join all intervals to one
1829       aU1f[0] = Min(aU1f[0], aU1f[1]);
1830       aU1l[0] = Max(aU1l[0], aU1l[1]);
1831
1832       aU1f[1] = -Precision::Infinite();
1833       aU1l[1] = Precision::Infinite();
1834     }
1835   }
1836
1837   //Critical points
1838   //[0...1] - in these points parameter U1 goes through
1839   //          the seam-edge of the first cylinder.
1840   //[2...3] - First and last U1 parameter.
1841   //[4...5] - in these points parameter U2 goes through
1842   //          the seam-edge of the second cylinder.
1843   //[6...9] - in these points an intersection line goes through
1844   //          U-boundaries of the second surface.
1845   const Standard_Integer aNbCritPointsMax = 10;
1846   Standard_Real anU1crit[aNbCritPointsMax] = {Precision::Infinite(),
1847                                               Precision::Infinite(),
1848                                               Precision::Infinite(),
1849                                               Precision::Infinite(),
1850                                               Precision::Infinite(),
1851                                               Precision::Infinite(),
1852                                               Precision::Infinite(),
1853                                               Precision::Infinite(),
1854                                               Precision::Infinite(),
1855                                               Precision::Infinite()};
1856
1857   CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
1858                                       aPeriod, theTol2D, aNbCritPointsMax, anU1crit);
1859
1860
1861   //Getting Walking-line
1862
1863   enum WLFStatus
1864   {
1865     WLFStatus_Absent = 0,
1866     WLFStatus_Exist  = 1,
1867     WLFStatus_Broken = 2
1868   };
1869
1870   for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++)
1871   {
1872     if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval]))
1873       continue;
1874
1875     Standard_Boolean isAddedIntoWL[aNbWLines];
1876     for(Standard_Integer i = 0; i < aNbWLines; i++)
1877       isAddedIntoWL[i] = Standard_False;
1878
1879     Standard_Real anUf = aU1f[aCurInterval], anUl = aU1l[aCurInterval];
1880     const Standard_Boolean isDeltaPeriod  = IsEqual(anUl-anUf, aPeriod);
1881
1882     //Inscribe and sort critical points
1883     InscribeAndSortArray(anU1crit, aNbCritPointsMax, anUf, anUl, theTol2D, aPeriod);
1884
1885     while(anUf < anUl)
1886     {
1887       Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
1888       WLFStatus aWLFindStatus[aNbWLines];
1889       Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
1890       Standard_Real anArccosFactor[aNbWLines] = {1.0, -1.0};
1891       Standard_Boolean isAddingWLEnabled[aNbWLines];
1892
1893       Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
1894       Handle(IntPatch_WLine) aWLine[aNbWLines];
1895       for(Standard_Integer i = 0; i < aNbWLines; i++)
1896       {
1897         aL2S[i] = new IntSurf_LineOn2S();
1898         aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
1899         aWLFindStatus[i] = WLFStatus_Absent;
1900         isAddingWLEnabled[i] = Standard_True;
1901         aU2[i] = aV1[i] = aV2[i] = 0.0;
1902         aV1Prev[i] = aV2Prev[i] = 0.0;
1903       }
1904
1905       Standard_Real anU1 = anUf;
1906
1907       Standard_Real aCriticalDelta[aNbCritPointsMax];
1908       for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
1909         aCriticalDelta[i] = anU1 - anU1crit[i];
1910
1911       Standard_Boolean isFirst = Standard_True;
1912
1913       while(anU1 <= anUl)
1914       {
1915         if(isDeltaPeriod)
1916         {
1917           if(IsEqual(anU1, anUl))
1918           {
1919             //if isAddedIntoWL* == TRUE WLine contains only one point
1920             //(which was end point of previous WLine). If we will
1921             //add point found on the current step WLine will contain only
1922             //two points. At that both these points will be equal to the
1923             //points found earlier. Therefore, new WLine will repeat
1924             //already existing WLine. Consequently, it is necessary
1925             //to forbid building new line in this case.
1926
1927             for(Standard_Integer i = 0; i < aNbWLines; i++)
1928               isAddingWLEnabled[i] = !isAddedIntoWL[i];
1929           }
1930         }
1931
1932         for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
1933         {
1934           if((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
1935           {
1936             anU1 = anU1crit[i];
1937
1938             for(Standard_Integer i = 0; i < aNbWLines; i++)
1939               aWLFindStatus[i] = WLFStatus_Broken;
1940
1941             break;
1942           }
1943         }
1944
1945         Standard_Real anArg = anEquationCoeffs.mB *
1946               cos(anU1 - anEquationCoeffs.mFI1) + anEquationCoeffs.mC;
1947
1948         if(aNulValue > 1.0 - anArg)
1949           anArg = 1.0;
1950         if(anArg + 1.0 < aNulValue)
1951           anArg = -1.0;
1952
1953         for(Standard_Integer i = 0; i < aNbWLines; i++)
1954         {
1955           const Standard_Integer aNbPntsWL =  aWLine[i].IsNull() ? 0 :
1956                                               aWLine[i]->Curve()->NbPoints();
1957           aU2[i] = anEquationCoeffs.mFI2 + anArccosFactor[i]*acos(anArg);
1958
1959           InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
1960
1961           if(aNbPntsWL == 0)
1962           {//the line has not contained any points yet
1963             if(((aUSurf2l - aUSurf2f) >= aPeriod) &&
1964                 ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
1965                   (Abs(aU2[i]-aUSurf2l) < theTol2D)))
1966             {
1967               const Standard_Real anU1Temp = anU1 + aStepMin;
1968               Standard_Real anArgTemp = anEquationCoeffs.mB *
1969                                         cos(anU1Temp - anEquationCoeffs.mFI1) +
1970                                         anEquationCoeffs.mC;
1971
1972               if(aNulValue > 1.0 - anArgTemp)
1973                 anArgTemp = 1.0;
1974
1975               if(anArgTemp + 1.0 < aNulValue)
1976                 anArgTemp = -1.0;
1977
1978               Standard_Real aU2Temp = anEquationCoeffs.mFI2 +
1979                                       anArccosFactor[i]*acos(anArgTemp);
1980
1981               InscribePoint(aUSurf2f, aUSurf2l, aU2Temp, theTol2D, aPeriod, Standard_False);
1982
1983               if(2.0*Abs(aU2Temp - aU2[i]) > aPeriod)
1984               {
1985                 if(aU2Temp > aU2[i])
1986                   aU2[i] += aPeriod;
1987                 else
1988                   aU2[i] -= aPeriod;
1989               }
1990             }
1991           }
1992
1993           if(aNbPntsWL > 0)
1994           {
1995             if(((aUSurf2l - aUSurf2f) >= aPeriod) &&
1996                 ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
1997                   (Abs(aU2[i]-aUSurf2l) < theTol2D)))
1998             {//end of the line
1999               Standard_Real aU2prev = 0.0, aV2prev = 0.0;
2000               if(isTheReverse)
2001                 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
2002               else
2003                 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
2004
2005               if(2.0*Abs(aU2prev - aU2[i]) > aPeriod)
2006               {
2007                 if(aU2prev > aU2[i])
2008                   aU2[i] += aPeriod;
2009                 else
2010                   aU2[i] -= aPeriod;
2011               }
2012             }
2013           }
2014
2015           if(aWLFindStatus[i] == WLFStatus_Broken)
2016           {
2017             if(Abs(aU2[i]) <= theTol2D)
2018               aU2[i] = 0.0;
2019             else if(Abs(aU2[i] - aPeriod) <= theTol2D)
2020               aU2[i] = aPeriod;
2021             else if(Abs(aU2[i] - aUSurf2f) <= theTol2D)
2022               aU2[i] = aUSurf2f;
2023             else if(Abs(aU2[i] - aUSurf2l) <= theTol2D)
2024               aU2[i] = aUSurf2l;
2025           }
2026
2027           aV1[i] =  anEquationCoeffs.mK21 * sin(aU2[i]) +
2028                     anEquationCoeffs.mK11 * sin(anU1) +
2029                     anEquationCoeffs.mL21 * cos(aU2[i]) +
2030                     anEquationCoeffs.mL11 * cos(anU1) + anEquationCoeffs.mM1;
2031
2032           aV2[i] =  anEquationCoeffs.mK22 * sin(aU2[i]) +
2033                     anEquationCoeffs.mK12 * sin(anU1) +
2034                     anEquationCoeffs.mL22 * cos(aU2[i]) +
2035                     anEquationCoeffs.mL12 * cos(anU1) + anEquationCoeffs.mM2;
2036
2037           if(isFirst)
2038           {
2039             aV1Prev[i] = aV1[i];
2040             aV2Prev[i] = aV2[i];
2041           }
2042         }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2043
2044         isFirst = Standard_False;
2045
2046         //Looking for points into WLine
2047         Standard_Boolean isBroken = Standard_False;
2048         for(Standard_Integer i = 0; i < aNbWLines; i++)
2049         {
2050           if(!isAddingWLEnabled[i])
2051           {
2052             aV1Prev[i] = aV1[i];
2053             aV2Prev[i] = aV2[i];
2054
2055             if(aWLFindStatus[i] == WLFStatus_Broken)
2056               isBroken = Standard_True;
2057
2058             continue;
2059           }
2060
2061           if( ((aUSurf2f-aU2[i]) <= theTol2D) && ((aU2[i]-aUSurf2l) <= theTol2D) &&
2062               ((aVSurf1f - aV1[i]) <= theTol2D) && ((aV1[i] - aVSurf1l) <= theTol2D) &&
2063               ((aVSurf2f - aV2[i]) <= theTol2D) && ((aV2[i] - aVSurf2l) <= theTol2D))
2064           {
2065             Standard_Boolean isForce = Standard_False;
2066             if(aWLFindStatus[i] == WLFStatus_Absent)
2067             {
2068               Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
2069
2070               if(((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1-aUSurf1l) < theTol2D))
2071               {
2072                 isForce = Standard_True;
2073               }
2074
2075               AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs,
2076                                 theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod,
2077                                 aNulValue, anU1, aU2[i], aV1[i], aV1Prev[i],
2078                                 aV2[i], aV2Prev[i], isTheReverse,
2079                                 anArccosFactor[i], isForce, isFound1, isFound2);
2080
2081               if(isFound1 || isFound2)
2082               {
2083                 aWLFindStatus[i] = WLFStatus_Exist;
2084               }
2085             }
2086
2087             if(( aWLFindStatus[i] != WLFStatus_Broken) || (aWLine[i]->NbPnts() >= 1))
2088             {
2089               if(AddPointIntoWL(theQuad1, theQuad2, isTheReverse,
2090                                 gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
2091                                 aUSurf1f, aUSurf1l, aPeriod,
2092                                 aWLine[i]->Curve(), theTol3D, theTol2D, isForce))
2093               {
2094                 if(aWLFindStatus[i] == WLFStatus_Absent)
2095                 {
2096                   aWLFindStatus[i] = WLFStatus_Exist;
2097                 }
2098               }
2099             }
2100           }
2101           else
2102           {
2103             if(aWLFindStatus[i] == WLFStatus_Exist)
2104             {
2105               Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
2106
2107               AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs,
2108                                 theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod,
2109                                 aNulValue, anU1, aU2[i], aV1[i], aV1Prev[i],
2110                                 aV2[i], aV2Prev[i], isTheReverse,
2111                                 anArccosFactor[i], Standard_False, isFound1, isFound2);
2112
2113               if(isFound1 || isFound2)
2114                 aWLFindStatus[i] = WLFStatus_Broken; //start a new line
2115             }
2116           }
2117
2118           aV1Prev[i] = aV1[i];
2119           aV2Prev[i] = aV2[i];
2120
2121           if(aWLFindStatus[i] == WLFStatus_Broken)
2122             isBroken = Standard_True;
2123         }//for(Standard_Integer i = 0; i < aNbWLines; i++)
2124
2125         if(isBroken)
2126         {//current lines are filled. Go to the next lines
2127           anUf = anU1;
2128           break;
2129         }
2130
2131         //Step computing
2132
2133         Standard_Real aStepU1 = aStepMax;
2134
2135         for(Standard_Integer i = 0; i < aNbWLines; i++)
2136         {
2137           Standard_Real aDeltaU1V1 = aStepU1, aDeltaU1V2 = aStepU1;
2138
2139           Standard_Real aFact1 = !IsEqual(sin(aU2[i] - anEquationCoeffs.mFI2), 0.0) ?
2140                                   anEquationCoeffs.mK1 * sin(anU1 - anEquationCoeffs.mFIV1) +
2141                                   anEquationCoeffs.mL1 * anEquationCoeffs.mB * sin(aU2[i] - anEquationCoeffs.mPSIV1) *
2142                                   sin(anU1 - anEquationCoeffs.mFI1)/sin(aU2[i]-anEquationCoeffs.mFI2) : 0.0,
2143                         aFact2 = !IsEqual(sin(aU2[i]-anEquationCoeffs.mFI2), 0.0) ?
2144                                   anEquationCoeffs.mK2 * sin(anU1 - anEquationCoeffs.mFIV2) +
2145                                   anEquationCoeffs.mL2 * anEquationCoeffs.mB * sin(aU2[i] - anEquationCoeffs.mPSIV2) *
2146                                   sin(anU1 - anEquationCoeffs.mFI1)/sin(aU2[i] - anEquationCoeffs.mFI2) : 0.0;
2147
2148           Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints),
2149                         aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints);
2150
2151           if((aV1[i] < aVSurf1f) && (aFact1 < 0.0))
2152           {//Make close to aVSurf1f by increasing anU1
2153             aDeltaV1 = Min(aDeltaV1, Abs(aV1[i]-aVSurf1f));
2154           }
2155
2156           if((aV1[i] > aVSurf1l) && (aFact1 > 0.0))
2157           {//Make close to aVSurf1l by increasing anU1
2158             aDeltaV1 = Min(aDeltaV1, Abs(aV1[i]-aVSurf1l));
2159           }
2160
2161           if((aV2[i] < aVSurf2f) && (aFact2 < 0.0))
2162           {//Make close to aVSurf2f by increasing anU1
2163             aDeltaV2 = Min(aDeltaV2, Abs(aV2[i]-aVSurf2f));
2164           }
2165
2166           if((aV2[i] > aVSurf2l) && (aFact2 > 0.0))
2167           {//Make close to aVSurf2l by increasing anU1
2168             aDeltaV2 = Min(aDeltaV2, Abs(aV2[i]-aVSurf1l));
2169           }
2170
2171           aDeltaU1V1 = !IsEqual(aFact1,0.0)? Abs(aDeltaV1/aFact1) : aStepMax;
2172           aDeltaU1V2 = !IsEqual(aFact2,0.0)? Abs(aDeltaV2/aFact2) : aStepMax;
2173
2174           if(aDeltaU1V1 < aStepU1)
2175             aStepU1 = aDeltaU1V1;
2176
2177           if(aDeltaU1V2 < aStepU1)
2178             aStepU1 = aDeltaU1V2;
2179         }
2180
2181         if(aStepU1 < aStepMin)
2182           aStepU1 = aStepMin;
2183
2184         if(aStepU1 > aStepMax)
2185           aStepU1 = aStepMax;
2186
2187         anU1 += aStepU1;
2188
2189         const Standard_Real aDiff = anU1 - anUl;
2190         if((0.0 < aDiff) && (aDiff < aStepU1-Precision::PConfusion()))
2191           anU1 = anUl;
2192
2193         anUf = anU1;
2194
2195         for(Standard_Integer i = 0; i < aNbWLines; i++)
2196         {
2197           if(aWLine[i]->NbPnts() != 1)
2198             isAddedIntoWL[i] = Standard_False;
2199         }
2200       }
2201
2202       for(Standard_Integer i = 0; i < aNbWLines; i++)
2203       {
2204         if((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
2205         {
2206           isTheEmpty = Standard_False;
2207           Standard_Real u1, v1, u2, v2;
2208           aWLine[i]->Point(1).Parameters(u1, v1, u2, v2);
2209           IntPatch_Point aP;
2210           aP.SetParameter(u1);
2211           aP.SetParameters(u1, v1, u2, v2);
2212           aP.SetTolerance(theTol3D);
2213           aP.SetValue(aWLine[i]->Point(1).Value());
2214
2215           theSPnt.Append(aP);
2216         }
2217         else if(aWLine[i]->NbPnts() > 1)
2218         {
2219           Standard_Boolean isGood = Standard_True;
2220
2221           if(aWLine[i]->NbPnts() == 2)
2222           {
2223             const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
2224             const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
2225
2226             if(aPf.IsSame(aPl, Precision::Confusion()))
2227               isGood = Standard_False;
2228           }
2229
2230           if(isGood)
2231           {
2232             isTheEmpty = Standard_False;
2233             isAddedIntoWL[i] = Standard_True;
2234             SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(),
2235                                   anEquationCoeffs, aNbPoints, aUSurf2f, aUSurf2l,
2236                                   theTol2D, aPeriod, anArccosFactor[i], isTheReverse);
2237
2238             aWLine[i]->ComputeVertexParameters(theTol3D);
2239             theSlin.Append(aWLine[i]);
2240           }
2241         }
2242         else
2243         {
2244           isAddedIntoWL[i] = Standard_False;
2245         }
2246       }
2247     }
2248   }
2249
2250   return Standard_True;
2251 }
2252
2253 //=======================================================================
2254 //function : IntCySp
2255 //purpose  :
2256 //=======================================================================
2257 Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1,
2258                          const IntSurf_Quadric& Quad2,
2259                          const Standard_Real Tol,
2260                          const Standard_Boolean Reversed,
2261                          Standard_Boolean& Empty,
2262                          Standard_Boolean& Multpoint,
2263                          IntPatch_SequenceOfLine& slin,
2264                          IntPatch_SequenceOfPoint& spnt)
2265
2266 {
2267   Standard_Integer i;
2268
2269   IntSurf_TypeTrans trans1,trans2;
2270   IntAna_ResultType typint;
2271   IntPatch_Point ptsol;
2272   gp_Circ cirsol;
2273
2274   gp_Cylinder Cy;
2275   gp_Sphere Sp;
2276
2277   if (!Reversed) {
2278     Cy = Quad1.Cylinder();
2279     Sp = Quad2.Sphere();
2280   }
2281   else {
2282     Cy = Quad2.Cylinder();
2283     Sp = Quad1.Sphere();
2284   }
2285   IntAna_QuadQuadGeo inter(Cy,Sp,Tol);
2286
2287   if (!inter.IsDone()) {return Standard_False;}
2288
2289   typint = inter.TypeInter();
2290   Standard_Integer NbSol = inter.NbSolutions();
2291   Empty = Standard_False;
2292
2293   switch (typint) {
2294
2295   case IntAna_Empty :
2296     {
2297       Empty = Standard_True;
2298     }
2299     break;
2300
2301   case IntAna_Point :
2302     {
2303       gp_Pnt psol(inter.Point(1));
2304       Standard_Real U1,V1,U2,V2;
2305       Quad1.Parameters(psol,U1,V1);
2306       Quad2.Parameters(psol,U2,V2);
2307       ptsol.SetValue(psol,Tol,Standard_True);
2308       ptsol.SetParameters(U1,V1,U2,V2);
2309       spnt.Append(ptsol);
2310     }
2311     break;
2312
2313   case IntAna_Circle:
2314     {
2315       cirsol = inter.Circle(1);
2316       gp_Vec Tgt;
2317       gp_Pnt ptref;
2318       ElCLib::D1(0.,cirsol,ptref,Tgt);
2319
2320       if (NbSol == 1) {
2321         gp_Vec TestCurvature(ptref,Sp.Location());
2322         gp_Vec Normsp,Normcyl;
2323         if (!Reversed) {
2324           Normcyl = Quad1.Normale(ptref);
2325           Normsp  = Quad2.Normale(ptref);
2326         }
2327         else {
2328           Normcyl = Quad2.Normale(ptref);
2329           Normsp  = Quad1.Normale(ptref);
2330         }
2331
2332         IntSurf_Situation situcyl;
2333         IntSurf_Situation situsp;
2334
2335         if (Normcyl.Dot(TestCurvature) > 0.) {
2336           situsp = IntSurf_Outside;
2337           if (Normsp.Dot(Normcyl) > 0.) {
2338             situcyl = IntSurf_Inside;
2339           }
2340           else {
2341             situcyl = IntSurf_Outside;
2342           }
2343         }
2344         else {
2345           situsp = IntSurf_Inside;
2346           if (Normsp.Dot(Normcyl) > 0.) {
2347             situcyl = IntSurf_Outside;
2348           }
2349           else {
2350             situcyl = IntSurf_Inside;
2351           }
2352         }
2353         Handle(IntPatch_GLine) glig;
2354         if (!Reversed) {
2355           glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp);
2356         }
2357         else {
2358           glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl);
2359         }
2360         slin.Append(glig);
2361       }
2362       else {
2363         if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) {
2364           trans1 = IntSurf_Out;
2365           trans2 = IntSurf_In;
2366         }
2367         else {
2368           trans1 = IntSurf_In;
2369           trans2 = IntSurf_Out;
2370         }
2371         Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2372         slin.Append(glig);
2373
2374         cirsol = inter.Circle(2);
2375         ElCLib::D1(0.,cirsol,ptref,Tgt);
2376         Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
2377         if(qwe> 0.0000001) {
2378           trans1 = IntSurf_Out;
2379           trans2 = IntSurf_In;
2380         }
2381         else if(qwe<-0.0000001) {
2382           trans1 = IntSurf_In;
2383           trans2 = IntSurf_Out;
2384         }
2385         else {
2386           trans1=trans2=IntSurf_Undecided;
2387         }
2388         glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2389         slin.Append(glig);
2390       }
2391     }
2392     break;
2393
2394   case IntAna_NoGeometricSolution:
2395     {
2396       gp_Pnt psol;
2397       Standard_Real U1,V1,U2,V2;
2398       IntAna_IntQuadQuad anaint(Cy,Sp,Tol);
2399       if (!anaint.IsDone()) {
2400         return Standard_False;
2401       }
2402
2403       if (anaint.NbPnt()==0 && anaint.NbCurve()==0) {
2404         Empty = Standard_True;
2405       }
2406       else {
2407
2408         NbSol = anaint.NbPnt();
2409         for (i = 1; i <= NbSol; i++) {
2410           psol = anaint.Point(i);
2411           Quad1.Parameters(psol,U1,V1);
2412           Quad2.Parameters(psol,U2,V2);
2413           ptsol.SetValue(psol,Tol,Standard_True);
2414           ptsol.SetParameters(U1,V1,U2,V2);
2415           spnt.Append(ptsol);
2416         }
2417
2418         gp_Pnt ptvalid,ptf,ptl;
2419         gp_Vec tgvalid;
2420         Standard_Real first,last,para;
2421         IntAna_Curve curvsol;
2422         Standard_Boolean tgfound;
2423         Standard_Integer kount;
2424
2425         NbSol = anaint.NbCurve();
2426         for (i = 1; i <= NbSol; i++) {
2427           curvsol = anaint.Curve(i);
2428           curvsol.Domain(first,last);
2429           ptf = curvsol.Value(first);
2430           ptl = curvsol.Value(last);
2431
2432           para = last;
2433           kount = 1;
2434           tgfound = Standard_False;
2435
2436           while (!tgfound) {
2437             para = (1.123*first + para)/2.123;
2438             tgfound = curvsol.D1u(para,ptvalid,tgvalid);
2439             if (!tgfound) {
2440               kount ++;
2441               tgfound = kount > 5;
2442             }
2443           }
2444           Handle(IntPatch_ALine) alig;
2445           if (kount <= 5) {
2446             Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
2447                                                  Quad1.Normale(ptvalid));
2448             if(qwe> 0.00000001) {
2449               trans1 = IntSurf_Out;
2450               trans2 = IntSurf_In;
2451             }
2452             else if(qwe<-0.00000001) {
2453               trans1 = IntSurf_In;
2454               trans2 = IntSurf_Out;
2455             }
2456             else {
2457               trans1=trans2=IntSurf_Undecided;
2458             }
2459             alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
2460           }
2461           else {
2462             alig = new IntPatch_ALine(curvsol,Standard_False);
2463           }
2464           Standard_Boolean TempFalse1a = Standard_False;
2465           Standard_Boolean TempFalse2a = Standard_False;
2466
2467           //-- ptf et ptl : points debut et fin de alig
2468
2469           ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first,
2470                         TempFalse2a,ptl,last,Multpoint,Tol);
2471           slin.Append(alig);
2472         } //-- boucle sur les lignes
2473       } //-- solution avec au moins une lihne
2474     }
2475     break;
2476
2477   default:
2478     {
2479       return Standard_False;
2480     }
2481   }
2482   return Standard_True;
2483 }
2484 //=======================================================================
2485 //function : IntCyCo
2486 //purpose  :
2487 //=======================================================================
2488 Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
2489                          const IntSurf_Quadric& Quad2,
2490                          const Standard_Real Tol,
2491                          const Standard_Boolean Reversed,
2492                          Standard_Boolean& Empty,
2493                          Standard_Boolean& Multpoint,
2494                          IntPatch_SequenceOfLine& slin,
2495                          IntPatch_SequenceOfPoint& spnt)
2496
2497 {
2498   IntPatch_Point ptsol;
2499
2500   Standard_Integer i;
2501
2502   IntSurf_TypeTrans trans1,trans2;
2503   IntAna_ResultType typint;
2504   gp_Circ cirsol;
2505
2506   gp_Cylinder Cy;
2507   gp_Cone     Co;
2508
2509   if (!Reversed) {
2510     Cy = Quad1.Cylinder();
2511     Co = Quad2.Cone();
2512   }
2513   else {
2514     Cy = Quad2.Cylinder();
2515     Co = Quad1.Cone();
2516   }
2517   IntAna_QuadQuadGeo inter(Cy,Co,Tol);
2518
2519   if (!inter.IsDone()) {return Standard_False;}
2520
2521   typint = inter.TypeInter();
2522   Standard_Integer NbSol = inter.NbSolutions();
2523   Empty = Standard_False;
2524
2525   switch (typint) {
2526
2527   case IntAna_Empty : {
2528     Empty = Standard_True;
2529   }
2530     break;
2531
2532   case IntAna_Point :{
2533     gp_Pnt psol(inter.Point(1));
2534     Standard_Real U1,V1,U2,V2;
2535     Quad1.Parameters(psol,U1,V1);
2536     Quad1.Parameters(psol,U2,V2);
2537     ptsol.SetValue(psol,Tol,Standard_True);
2538     ptsol.SetParameters(U1,V1,U2,V2);
2539     spnt.Append(ptsol);
2540   }
2541     break;
2542
2543   case IntAna_Circle:  {
2544     gp_Vec Tgt;
2545     gp_Pnt ptref;
2546     Standard_Integer j;
2547     Standard_Real qwe;
2548     //
2549     for(j=1; j<=2; ++j) {
2550       cirsol = inter.Circle(j);
2551       ElCLib::D1(0.,cirsol,ptref,Tgt);
2552       qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
2553       if(qwe> 0.00000001) {
2554         trans1 = IntSurf_Out;
2555         trans2 = IntSurf_In;
2556       }
2557       else if(qwe<-0.00000001) {
2558         trans1 = IntSurf_In;
2559         trans2 = IntSurf_Out;
2560       }
2561       else {
2562         trans1=trans2=IntSurf_Undecided;
2563       }
2564       Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2);
2565       slin.Append(glig);
2566     }
2567   }
2568     break;
2569
2570   case IntAna_NoGeometricSolution:    {
2571     gp_Pnt psol;
2572     Standard_Real U1,V1,U2,V2;
2573     IntAna_IntQuadQuad anaint(Cy,Co,Tol);
2574     if (!anaint.IsDone()) {
2575       return Standard_False;
2576     }
2577
2578     if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) {
2579       Empty = Standard_True;
2580     }
2581     else {
2582       NbSol = anaint.NbPnt();
2583       for (i = 1; i <= NbSol; i++) {
2584         psol = anaint.Point(i);
2585         Quad1.Parameters(psol,U1,V1);
2586         Quad2.Parameters(psol,U2,V2);
2587         ptsol.SetValue(psol,Tol,Standard_True);
2588         ptsol.SetParameters(U1,V1,U2,V2);
2589         spnt.Append(ptsol);
2590       }
2591
2592       gp_Pnt ptvalid, ptf, ptl;
2593       gp_Vec tgvalid;
2594       //
2595       Standard_Real first,last,para;
2596       Standard_Boolean tgfound,firstp,lastp,kept;
2597       Standard_Integer kount;
2598       //
2599       //
2600       //IntAna_Curve curvsol;
2601       IntAna_Curve aC;
2602       IntAna_ListOfCurve aLC;
2603       IntAna_ListIteratorOfListOfCurve aIt;
2604
2605       //
2606       NbSol = anaint.NbCurve();
2607       for (i = 1; i <= NbSol; ++i) {
2608         kept = Standard_False;
2609         //curvsol = anaint.Curve(i);
2610         aC=anaint.Curve(i);
2611         aLC.Clear();
2612         ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
2613         //
2614         aIt.Initialize(aLC);
2615         for (; aIt.More(); aIt.Next()) {
2616           IntAna_Curve& curvsol=aIt.Value();
2617           //
2618           curvsol.Domain(first, last);
2619           firstp = !curvsol.IsFirstOpen();
2620           lastp  = !curvsol.IsLastOpen();
2621           if (firstp) {
2622             ptf = curvsol.Value(first);
2623           }
2624           if (lastp) {
2625             ptl = curvsol.Value(last);
2626           }
2627           para = last;
2628           kount = 1;
2629           tgfound = Standard_False;
2630
2631           while (!tgfound) {
2632             para = (1.123*first + para)/2.123;
2633             tgfound = curvsol.D1u(para,ptvalid,tgvalid);
2634             if (!tgfound) {
2635               kount ++;
2636               tgfound = kount > 5;
2637             }
2638           }
2639           Handle(IntPatch_ALine) alig;
2640           if (kount <= 5) {
2641             Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid),
2642                                                  Quad1.Normale(ptvalid));
2643             if(qwe> 0.00000001) {
2644               trans1 = IntSurf_Out;
2645               trans2 = IntSurf_In;
2646             }
2647             else if(qwe<-0.00000001) {
2648               trans1 = IntSurf_In;
2649               trans2 = IntSurf_Out;
2650             }
2651             else {
2652               trans1=trans2=IntSurf_Undecided;
2653             }
2654             alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
2655             kept = Standard_True;
2656           }
2657           else {
2658             ptvalid = curvsol.Value(para);
2659             alig = new IntPatch_ALine(curvsol,Standard_False);
2660             kept = Standard_True;
2661             //-- cout << "Transition indeterminee" << endl;
2662           }
2663           if (kept) {
2664             Standard_Boolean Nfirstp = !firstp;
2665             Standard_Boolean Nlastp  = !lastp;
2666             ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first,
2667                           Nlastp,ptl,last,Multpoint,Tol);
2668             slin.Append(alig);
2669           }
2670         } // for (; aIt.More(); aIt.Next())
2671       } // for (i = 1; i <= NbSol; ++i)
2672     }
2673   }
2674     break;
2675
2676   default:
2677     return Standard_False;
2678
2679   } // switch (typint)
2680
2681   return Standard_True;
2682 }
2683 //=======================================================================
2684 //function : ExploreCurve
2685 //purpose  :
2686 //=======================================================================
2687 Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
2688                               const gp_Cone& aCo,
2689                               IntAna_Curve& aC,
2690                               const Standard_Real aTol,
2691                               IntAna_ListOfCurve& aLC)
2692
2693 {
2694   Standard_Boolean bFind=Standard_False;
2695   Standard_Real aTheta, aT1, aT2, aDst;
2696   gp_Pnt aPapx, aPx;
2697   //
2698   //aC.Dump();
2699   //
2700   aLC.Clear();
2701   aLC.Append(aC);
2702   //
2703   aPapx=aCo.Apex();
2704   //
2705   aC.Domain(aT1, aT2);
2706   //
2707   aPx=aC.Value(aT1);
2708   aDst=aPx.Distance(aPapx);
2709   if (aDst<aTol) {
2710     return bFind;
2711   }
2712   aPx=aC.Value(aT2);
2713   aDst=aPx.Distance(aPapx);
2714   if (aDst<aTol) {
2715     return bFind;
2716   }
2717   //
2718   bFind=aC.FindParameter(aPapx, aTheta);
2719   if (!bFind){
2720     return bFind;
2721   }
2722   //
2723   aPx=aC.Value(aTheta);
2724   aDst=aPx.Distance(aPapx);
2725   if (aDst>aTol) {
2726     return !bFind;
2727   }
2728   //
2729   // need to be splitted at aTheta
2730   IntAna_Curve aC1, aC2;
2731   //
2732   aC1=aC;
2733   aC1.SetDomain(aT1, aTheta);
2734   aC2=aC;
2735   aC2.SetDomain(aTheta, aT2);
2736   //
2737   aLC.Clear();
2738   aLC.Append(aC1);
2739   aLC.Append(aC2);
2740   //
2741   return bFind;
2742 }
2743 //=======================================================================
2744 //function : IsToReverse
2745 //purpose  :
2746 //=======================================================================
2747 Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
2748                              const gp_Cylinder& Cy2,
2749                              const Standard_Real Tol)
2750 {
2751   Standard_Boolean bRet;
2752   Standard_Real aR1,aR2, dR, aSc1, aSc2;
2753   //
2754   bRet=Standard_False;
2755   //
2756   aR1=Cy1.Radius();
2757   aR2=Cy2.Radius();
2758   dR=aR1-aR2;
2759   if (dR<0.) {
2760     dR=-dR;
2761   }
2762   if (dR>Tol) {
2763     bRet=aR1>aR2;
2764   }
2765   else {
2766     gp_Dir aDZ(0.,0.,1.);
2767     //
2768     const gp_Dir& aDir1=Cy1.Axis().Direction();
2769     aSc1=aDir1*aDZ;
2770     if (aSc1<0.) {
2771       aSc1=-aSc1;
2772     }
2773     //
2774     const gp_Dir& aDir2=Cy2.Axis().Direction();
2775     aSc2=aDir2*aDZ;
2776     if (aSc2<0.) {
2777       aSc2=-aSc2;
2778     }
2779     //
2780     if (aSc2<aSc1) {
2781       bRet=!bRet;
2782     }
2783   }//if (dR<Tol) {
2784   return bRet;
2785 }