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