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