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