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